HITS-CLIP Unravels microRNA-mRNA Interactions

September 3, 2009 by Dan Koboldt

Micro-RNAs (miRNAs) are short (18-26 nt) sequences that act as post-transcriptional repressors of gene expression. Over 700 miRNAs have been reported in the human genome; each is believed to bind directly to many mRNAs to regulate their translation or stability. Thus, miRNAs represent a key regulatory mechanism affecting numerous cellular activities, and are of particular interest in cancer research. Understanding the complex relationships between miRNAs and mRNAs remains challenging, however, and computational approaches alone have been largely unsuccessful.

HITS-CLIP: Isolation and Sequencing of Argonaute-miRNA-mRNA Complexes

Enter HITS-CLIP, a new approach that applies high throughput sequencing of RNAs isolated by crosslinking immunoprecipitation. Essentially, it’s a method by which radition is used to cross-link protein-RNA complexes and stringently purify them. Then, massively parallel sequencing yields all of the RNA “tags” bound by the protein of interest.

Ago-miRNA-mRNA Complex (Image Credit: Nature 460: 479-486, 2009)

In a recent Nature paper, Chi et al used HITS-CLIP to isolate RNA bound by the Argonaute protein (Ago), which mediates miRNA-mRNA interaction (see figure). The purified complexes showed two different modal sizes (110 kDa and 130kDa), suggesting that Ago (97 kDa) was crosslinked to two different RNA species – hopefully, miRNAs (small) and the mRNAs that they were targeting (large).

The authors applied Illumina high-throughput sequencing to characterize Ago-bound miRNAs and the mRNA “tags” to which they were linked. With relatively straightforward bio-informatics approaches, it was possible to cross-reference expressed miRNAs with complementary sequences of mRNA tags. The resulting “ternary map” of miRNA-mRNA interaction sites yields a wealth of information about this post-transcriptional regulatory mechanism.

Decoding miRNA-mRNA Interaction

The authors identified 454 unique miRNAs crosslinked to Ago in the mouse brain; mir-30e was the most abundant species, representing 14% of all miRNA tags. In silico clustering and normalization of the messenger RNA tags yielded 1,463 robust clusters from 829 different brain transcripts.

Locations of Ago-bound mRNA tags (Image Credit: Nature 460: 479-486, 2009)

When these tags were overlaid with gene annotations, several patterns emerged. As expected, a substantial portion (40%) of Ago-bound tags were in 3′ UTRs where miRNA activity is known to have high efficacy. Some 8% (one-fifth of the 40%) were actually outside of the UTR but <10kb downstream, regions likely to harbor unannotated 3′ UTRs.

Unsurprisingly, very few Ago-bound tags were in 5′ UTRs. However, a substantial fraction of tags fell in coding sequences (25%), introns (12%), and non-coding RNAs (4%), suggesting that miRNA activity occurs in these regions as well. Another 6% of tags were in intergenic regions, possibly in as-yet-unannotated transcripts. These unexpected locations of miRNA binding may offer additional insights into the mechanisms of miRNA regulation.

Next, the authors sought to define the Ago-mRNA “footprint” in which the majority of tags were contained. The distribution of tags in a defined cluster, at least in their figure, looks like a bell curve, with a sharp peak in the middle. About 95% of the time, Ago bound within 45-62 nucleotides of this peak, so the authors defined this region as the average Ago-miRNA footprint. Linear regression analysis of all 6-8 base motifs in clusters yielded numerous “enriched” seed sequences; the most prevalent corresponded to the binding site of miR-124, a well known brain-specific miRNA. Indeed, Ago-mRNA footprints were rich in miRNA binding sites, suggesting that this approach may predict active sites with far better specificity than other methods.

HITS-CLIP Implications

By reducing the search space for miRNA binding sites to a 45-60-nucleotide Ago footprint, HITS-CLIP offers a powerful complementary approach to bioinformatic methods for miRNA binding site prediction. Computational approaches alone are known to have high false positive rates, whereas the authors estimate FP rates of just 15% for HITS-CLIP. The new method offers dramatic improvement for transcripts with highly conserved 3′ UTRs, which often have many “predicted” miRNA binding sites because so many computational methods rely on conservation. Analysis of the HITS-CLIP ternary map revealed that real miRNA-mRNA binding events are very specific, with an average of just 2.6 Ago-mRNA clusters per regulated transcript. Despite the thousands of predicted binding sites, each miRNA bound an average of 655 targets. These results suggest that miRNA selectivity is much higher than previously believed. Yet Ago-mRNA clusters seemed to show no apparent sequence preference (data not shown), so it’s likely that other RNA-binding proteins are involved.

Thus, this study sets the stage for large scale genome-wide RNA-protein maps that include other proteins, tissues, and species, which should yield an unprecedented new level of understanding of this complex regulatory process.

References
Chi SW, Zang JB, Mele A, & Darnell RB (2009). Argonaute HITS-CLIP decodes microRNA-mRNA interaction maps. Nature, 460 (7254), 479-86 PMID: 19536157

WUCGI: WashU Cancer Genomics Initiative

August 27, 2009 by Dan Koboldt

WUCGI

Yesterday afternoon was the kickoff party launching WashU’s Cancer Genomics Initiative (CGI), better known as our goal to sequence 150 cancer genomes in the coming year.

Cancer Sequencing Ramps Up

Under the leadership of Wilson, Ley, and Elaine Mardis, our group sequenced the first cancer genome, from a woman who had died of AML (M1), and published the results in Nature last fall. Three weeks ago came the sequel. In the New England Journal of Medicine, we published the complete genome of another M1 leukemia, this time from a man who’s been treated and remains in full remission. In less than a year, the number of Illumina runs required to sequence a tumor genome dropped by over 80%, from 98 runs in AML1 to just 16.5 runs in AML2.

It’s not just the sequencing throughput that makes WUCGI a realistic effort. Many groups have Illumina sequencers, some even more than we do. Some of the most critical advances have taken place behind the scenes – for example, the variant detection pipelines developed by David Larson, Ken Chen, Chris Harris, and others. Sequencing on this scale would not be possible without the IT and informatics infrastructure, built under the leadership of David Dooling and Gary Stiehr, that gives us the computational firepower to run whole-genome analyses.

Two Genomes Down, 150 To Go

With two genomes published, the center leadership has set an ambitious goal: To sequence 150 cancer genomes in the coming year. Obviously, these will include more AML samples, hopefully some with therapy-related changes or abnormal cytogenetics. In collaboration with Matt Ellis and others at the Siteman Cancer Center, we’ll be tackling breast cancer as well. No doubt we’ll be revisiting lung cancer, for which we sequenced candidate genes as part of the Tumor Sequencing Project (TSP) consortium. As part of the Cancer Genome Atlas (TCGA) consortium, we’re working on glioblastoma multiforme (brain cancer) and ovarian cancer. Also, intriguingly, I hear rumors that there will be some sequencing of less common, largely unexplored cancers like multiple myeloma.

As Tim Ley said yesterday, it’s thrilling to be a part of this. We truly are entering the golden age of cancer genomics.

Second Cancer Genome in New England Journal

August 6, 2009 by Dan Koboldt

Today our group published the second cancer genome, AML2, in the New England Journal of Medicine. In this study, we sequenced the complete genomes of tumor cells and matched normal (skin) cells from a patient with cytogenetically normal de novo FAB M1 AML. This is an exciting publication for many reasons, the foremost of which may be the venue: with an impact factor of 52.59, the NEJM is almost certainly the most widely read biomedical journal in the world.

Diagnosed with Leukemia: It Could Happen to You

The story begins three years ago, with a previously healthy 38-year-old man of European ancestry who went to his doctor complaining of fatigue and a persistent cough. After an elevated white blood cell count, his physician ordered a bone marrow biopsy, which revealed 90% cellularity and 86% blasts. Diagnosis: Leukemia.

The patient underwent ten days of chemotherapy with cytarabine (7 days) followed by daunorubicin (3 days). Five weeks later he’d obtained morphologically complete remission and recovered counts. Now, three years later, he remains in complete remission. According to my conversations with an oncologist, this kind of happy ending is not very common with leukemia. Most leukemia patients are diagnosed at an advanced age, and don’t do as well.

Acute myelogenous leukemia cells. Credit: Univ. of Virginia

Moving Beyond Cytogenetics

At the time of his diagnosis, routine cytogenetic analysis of the patient’s tumor cells showed a normal 46XY karyotype. Bone marrow and skin samples were banked with informed consent for whole genome sequencing in accordance with our IRB. There was no family history of leukemia, though the patient’s mother had developed breast cancer and later non-Hodgkins lymphoma. Her half sister had also developed breast cancer. The field for discovery of mutations underlying this AML was wide open.

Whole Genome Sequencing with Illumina

We sequenced the genomes of tumor cells and matched normal (skin) cells to high depth (23.3x and 21.3x, respectively) on the Illumina/Solexa platform. The tumor sample required just 16.5 runs (most of which were 2×75 PE) to reach 98% diploid coverage. That’s a dramatic improvement over our first cancer genome, AML1, which took 98 runs (36 bp SE) to achieve 91% diploid coverage. At current rates, we really can sequence a genome a week. As any bioinformatician knows, however, the analysis usually takes a bit longer.

Dave Larson in my group really deserves the credit for the whole genome variant detection pipeline applied to AML2. With direction from Elaine Mardis, Rick Wilson, and Tim Ley, and others, Dave created a pipeline for automated variant calling, somatic scoring, and tiered classification of variants for cancer genomes (see Figure 1 of the paper). We identified 3.87 million single nucleotide variants (SNVs) in the tumor genome, of which 97.5% were in the skin genome and another 1.7% were previously described (i.e. dbSNP). That left 20,256 putative somatic variants which we classified as follows:

Tier 1 variants were coding variants that alter amino acid sequences, like nonsynonymous, nonstop, and splice-site mutations.
Tier 2 variants were variants in evolutionarily conserved or regulatory-potential sequences of the genome.
Tier 3 were the remaining variants that were in non-repetitive regions of the genome.
Tier 4 were the remaining variants that were in repetitive regions of the genome.

Validation and Deep 454 Read Counts

We used 3730 sequencing to validate somatic variants in Tiers 1 and 2. Some 62 mutations were validated, of which 10 were tier 1 (amino acid-altering) mutations. Additionally, we validated two somatic indels, one of which (NPM1) was previously described; the other was an insertion in the CEP170 gene predicted to add a leucine residue to the encoded protein.

In the absence of true functional validation, there are at least two approaches to evaluating whether or not a somatic mutation is a driver – a mutation that confers some advantage to drive tumor development – or a passenger – a background mutation that’s just along for the ride. First, driver mutations should be present most tumor cells, since the dominant clone will be the most “fit” in the tumor population. To assess mutation frequencies in our patient’s tumor cells, we applied 454 sequencing of mutation-containing amplicons in the tumor DNA, tumor cDNA, and skin DNA. Deep read counts for somatic events on the X and Y chromosomes showed allele frequencies of around 98%, consistent with the fact that nearly all cells in the bone marrow sample were part of the malignant clone. For the rest of the somatic mutations, variant frequencies hovered near the 50% mark (as expected) with a few exceptions. The CEP170 indel had a reduced (~35%) frequency in tumor DNA, suggesting that perhaps it’s not a driver mutation.

Recurrence of Mutations in Other AMLs

The other measure of importance of a somatic mutation is recurrence in other tumors of the same type. Thus, we screened for the presence of validated somatic mutations in a panel of 187 additional leukemia patients to see if any were recurrent. Most, unfortunately, were not. However, two variants were found in other samples, suggesting an important role in the development of AML. One was a noncoding conserved mutation (tier 2) on chromosome 10 which was detected in one other sample. Recurrence in just one other sample might not seem impressive, but by our estimation, the odds of such an event happening by chance are 1.1 x 10^-9. Thus, we may have uncovered a noncoding functional mutation that contributes to carcinogenesis via an as-yet-unknown mechanism.

The other was a nonsynonymous (tier 1) mutation in IDH1 at residue 132. Sixteen of 187 other leukemia samples carried mutations at the same residue in IDH1, suggesting an important role for this gene in the development of AML. Somatic mutations in IDH1 were recently characterized in glioblastoma (GBM) by our friends at Johns Hopkins, but this is the first time that IDH1 mutations were described in AML.

Conclusions: Lots of Passengers, Not Many Drivers

After sequencing the complete tumor genomes of two AML patients, we have estimated that these cancers carry an estimated 750 somatic events. Most such events will be background passenger mutations, acquired in the progenitor tumor cell before it became cancerous. Admittedly, that means there’s much more work to do to fully characterize the sequence changes underlying development of AML and other cancers. Our group is eager for the challenge. With the ever-growing throughput of the Illumina platform and our automated pipelines for whole-cancer-genome analysis, we hope to sequence at least a hundred more cancers in the coming year.

References

Mardis, E., Ding, L., Dooling, D., Larson, D., McLellan, M., Chen, K., Koboldt, D., Fulton, R., Delehaunty, K., McGrath, S., Fulton, L., Locke, D., Magrini, V., Abbott, R., Vickery, T., Reed, J., Robinson, J., Wylie, T., Smith, S., Carmichael, L., Eldred, J., Harris, C., Walker, J., Peck, J., Du, F., Dukes, A., Sanderson, G., Brummett, A., Clark, E., McMichael, J., Meyer, R., Schindler, J., Pohl, C., Wallis, J., Shi, X., Lin, L., Schmidt, H., Tang, Y., Haipek, C., Wiechert, M., Ivy, J., Kalicki, J., Elliott, G., Ries, R., Payton, J., Westervelt, P., Tomasson, M., Watson, M., Baty, J., Heath, S., Shannon, W., Nagarajan, R., Link, D., Walter, M., Graubert, T., DiPersio, J., Wilson, R., & Ley, T. (2009). Recurring Mutations Found by Sequencing an Acute Myeloid Leukemia Genome New England Journal of Medicine DOI: 10.1056/NEJMoa0903840

Maq, BWA, and Bowtie Compared

July 30, 2009 by Dan Koboldt

Until recently, Maq has provided the central alignment/assembly/variant-detection functionality for our Illumina pipeline. As technologies and algorithms evolve, however, we continue to investigate possible improvements. Heng Li’s sequel to Maq, called BWA, utilizes the incredibly fast Burrows-Wheeler indexing algorithm to speed up alignment time by orders of magnitude. Also, BWA generates alignments in SAM/BAM format by default, which is convenient for our large-scale sequencing projects where BAM files are becoming the standard format.

These features, along with our impression that Heng Li and company do not plan future updates to Maq, lead me to infer that BWA is the heir-apparent for our Illumina pipeline. Before the transition, however, we must compare Maq results with BWA results on the same dataset, to identify any differences that may affect downstream analysis. Also, we are continuing to evaluate other aligners, especially Bowtie, which offer comparable or even better speed at short read alignment.

Test Data: WGS and Targeted Sequencing of a Single Sample

We have a sample in-house for which we performed whole genome sequencing (WGS) and subsequently validated numerous novel variants. We also performed capture-based targeted resequencing (Illumina 2x75bp PE) of 6,000 genes in the same sample. To compare the performance of BWA, Maq, and Bowtie, we aligned the capture data with each tool separately, and looked at about a dozen sites where we’d validated novel variants from WGS.

Sensitivity – Total Reads Mapped

Here’s a histogram of the read depth at each of the 12 variant sites by aligner:

bwa-maq-bowtie-coverage

These results surprised me. Based on previous experience, I’d guessed that Maq would yield the highest depth, followed by BWA, and then Bowtie. Instead, with one exception, it was the other way around – Bowtie was more sensitive than BWA, which in turn was more sensitive than Maq. Yet these differences were relatively minor; overall, the coverage seems very comparable across all three aligners. I think that’s good news.

Variant Frequency by Read Count

Next, we looked at the observed variant frequencies, calculated as the relative fraction of reads supporting reference or variant alleles.

bwa-maq-bowtie-varfreq

When it comes to variant frequency, Maq and BWA yield almost identical results (despite slight coverage disparities). Bowtie yielded slighly higher frequencies in some cases, slightly lower frequencies in others. Again, these were very minor differences from three very different alignment algorithms, suggesting that each of them yields fairly robust results.

Farewell to Maq

Unfortunately, the results of my analysis do not bode well for Maq, only because Maq took a few days to align data that BWA and Bowtie processed in a matter of hours. So which Burrows-Wheeler aligner will prevail? It’s difficult to say. As far as SNP detection goes, BWA and Bowtie seem comparable.

« Previous Page