Transcriptome Genetics with HapMap and RNA-Seq

April 22, 2010 by Dan Koboldt

Two papers in Nature this month leverage the power of second-generation sequencing technologies to investigate gene expression variation in human cell lines. By performing RNA-Seq in HapMap cell lines, the authors generated the most extensive gene expression data to date for these samples, and were able to use publicly available HapMap genotypes to associate expression differences with genetic variation. This strategy was applied to the HapMap samples two years ago using expression microarrays. Using RNA-Seq instead of microarrays, however, offers a few key advantages:

More accurate quantification of highly abundant transcripts, where microarrays reach saturation
Access to rare transcripts below the sensitivity threshold for microarrays
Detection of novel gene structure from alternative splicing and unannotated exons
Identification of allele-specific expression

The first study, from Jonathan Pritchard’s lab at the University of Chicago, sequenced RNA from 69 Yoruban (African) individuals on the Illumina GAII platform. They generated at least two lanes per individual, for a total of 1.2 billion reads, of which 964 million (80%) mapped uniquely to the genome or to exon-exon boundaries. The second study, from Emmanouil Dermitzakis’s group at the Sanger center, sequenced RNA from 60 CEU (CEPH Europeans from Utah) individuals, also on the Illumina GAII platform. They generated one lane of paired-end data per individual, for a total of about 1.0 billion reads. Since neither study provided a table summarizing their data (which I’d have liked), I put one together:

Study	Pickrell et al.	Montgomery et al.
Samples	69, African descent (YRI)	60, European descent (CEU)
Sequencing	Illumina 1X35 or 1X46	Illumina 2X37
Reads/Sample	17.4 million	16.9 million
SNP Dataset	HapMap II/III	HapMap III
Total SNPs	3.8 million	1.2 million

Up to this point, the two studies sounded nearly identical. For the data analysis, however, each group went in a different (and interesting) direction.

Pooled Data for Discovery of Novel Gene Structures

Novel Exons. The Pritchard group pooled all data to examine the completeness of current gene annotations. Some 86% of uniquely mapped reads corresponded to known exons. Using conservation data from alignments of 28 vertebrate exomes, the authors identified 4,031 regions that are evolutionarily conserved and show evidence of transcription. About one-quarter of these appear to be part of spliced transcripts, but most appeared to be novel untranslated regions (UTRs). Some 115 regions, however, had sequences consistent with protein-coding exons. To investigate the possibility that their novel exons are real, the authors used RNA-Seq data from several human tissues and chimpanzee cell lines. The evidence suggests that their regions do represent novel exons, but ones that are expressed in a more tissue-specific fashion than annotated exons.

Novel Poly-A Sites. The authors next screened the ~70 million unmapped sequence reads for long runs of A or T nucleotides, which might indicate novel poly-adenylation sites. Of the ~8,000 novel sites that they identified, some 45% fell within 10 bp of a known cleavage site. To further validate their findings, they screened their poly-A regions for the binding site of the CPSF polyadenylation factor, and found a 32-fold enrichment for the CPSF target hexamer. The net result was a high confidence set of 3,481 cleavage sites that show evidence of poly-A (from RNA-Seq data) and CPSF binding.

RNA-Seq: 10 million Reads Is All You Need

The Dermitzakis study generated 16.9 million (+/- 5.9 million) reads per individual, which were mapped to the NCBI 36 reference sequence using Maq with a maximum insert size of 2 megabases). The resulting alignments were filtered to remove alignments with low mapping quality or to the X, Y, or MT chromosomes. Discordant read pairs (by distance or orientation) were also removed. To quantify the expression of known exons/transcripts/genes, the authors scaled read counts for each individual to a theoretical yield of 10 million reads, and only considered exons with data in >90% of individuals. This resulted in data for 90,064 exons from 10,777 genes, of which 95% had at least 10 reads (on average) per individual. While the normalization seems to reduce the dataset to less than half of known genes, it nevertheless provided an extensive view of gene expression across these 60 individuals.

Cis-Regulatory Effects on Gene Expression

Using HapMap genotypes for 1.2 million SNPs, the Dermitzakis group identified 836 genes associated with cis-regulatory variants (compared to 539 genes identified in microarray studies of the same individuals). Even when normalized for the number of genes tested, the increased resolution of RNA-Seq over microarrays yielded a larger number of genetic regulatory effects. The RNA-Seq exon eQTLs (expression quantitative trait loci) were enriched for abundant transcripts, suggesting that saturation of highly expressed exons reduces the sensitivity for microarrays to detect some cis-regulatory effects.

The Pritchard group searched for cis-regulatory variation with an even larger dataset – RNA-Seq for 69 individuals and 3.8 million HapMap SNPs. They identified 929 genes with local eQTLs (4.6% of annotated genes); consistent with previous findings, virtually all SNPs associated with expression level were near the corresponding gene. They also reported the overlap with the CEU study results: the top 500 associations reported in CEU samples were enriched 10 to 40-fold for significant eQTLs in YRI samples. Given the marked genetic differences between these two populations, this result suggests that these studies are identifying replicable cis-regulatory events.

Mechanism of Cis-Regulatory Effects

An important feature of RNA-Seq data is that it can be used not only to detect cis-regulatory variation, but to assess the mechanism by which these variants act. The Pritchard group looked at 222 of their 929 eQTLs for which the associated SNPs fell within the gene exons. They classified the RNA-Seq reads as originating from the high-expression haplotype or the low-expression haplotype, and found that for 195 of the genes (88%), more than 50% of the expressed transcripts carried the allele associated with high expression. Therefore, the modulation of gene expression is a direct result of the associated variation (probably by activating nearby cis-regulatory elements). In other words: the eQTL tells us that variants near the gene are associated with its expression. That means something nearby is regulating it. The fact that the haplotype associated with increased expression is the haplotype that predominates tells us that the high-expression allele is what drives the expression of its nearby gene. As opposed to, say, driving expression of the gene from both chromosomes.

Allelic Effects on Splicing

Finally, both groups looked at the actual content of expressed transcripts, to find SNPs associated with alternative splicing. The Pritchard group calls these splicing quantitative trait loci (sQTLs), and found 187 genes with significant associations. Binding sites for known splice factors (U1 snRNP and U2AF) were enriched for sQTLs, as were SNPs within 2 bp of a canonical splice site. The Dermitzakis group found 110 genes with significant associations, and stratified splicing-associated variants according to their position in the gene structure. When tested against the exons upstream and downstream of where they resided, splice donor variants were enriched 3.17-fold with the upstream (5′) exon, while splice acceptor variants were enriched 7.02 fold with the downstream (3′) exon. Thus, these SNPs affect the inclusion/exclusion of their exons in the mature transcript.

Dermitzakis’s group visually examined their most significant associations to characterize the mechanism of splicing regulation. Of the 110 significant sQTLs identified in CEU samples:

41% were single exon skipping events
17% created an alternate acceptor
13% were double or triple exon skipping events
6% created an alternate donor
5% were mutually exclusive exons
5% were retained introns.

In summary, these studies establish the feasibility of transcriptome sequencing to assess gene expression and characterize regulatory variation. Indeed, as the title of one study suggests, RNA sequencing is a powerful tool for studying the mechanisms underlying human gene expression variation, and will undoubtedly yield better understanding of the complex relationships between genotype and phenotype.

References
Pickrell JK, Marioni JC, Pai AA, Degner JF, Engelhardt BE, Nkadori E, Veyrieras JB, Stephens M, Gilad Y, & Pritchard JK (2010). Understanding mechanisms underlying human gene expression variation with RNA sequencing. Nature, 464 (7289), 768-72 PMID: 20220758

Montgomery SB, Sammeth M, Gutierrez-Arcelus M, Lach RP, Ingle C, Nisbett J, Guigo R, & Dermitzakis ET (2010). Transcriptome genetics using second generation sequencing in a Caucasian population. Nature, 464 (7289), 773-7 PMID: 20220756

Next-Gen Sequencing in 2010

March 9, 2010 by Dan Koboldt

On the shuttle from Marco Island to the airport last week, I happened to sit next to a very nice gentleman from Illumina. We got to talking, of course, and I asked him if they saw a threat from any of the new sequencing platforms presented at AGBT. I’m aware that Illumina currently enjoys a greater-than-50% share of the next-gen sequencing market, so I was curious about his impressions.

“We definitely see a segmentation of the market,” he admitted.

Something had been bothering me about the sequencing-company presentations this year, and I finally realized what it was. During AGBT 2009, every player was gunning to take over the world. This year it seems like every sequencing platform has a niche in mind.

General Sequencing: Illumina vs. Life Technologies

Illumina’s HiSeq2000 and Life Tech’s SOLiD 4 are after the general sequencing market – whole genome, transcriptome, and targeted (capture) sequencing. It’s a constant game of one-upmanship in throughput and claimed accuracy. In February this year, Illumina launched the HiSeq2000 with expected throughput of 200 GB per run. Life Technologies launched SOLiD 4 with 100 GB per run, but promised 300GB per run later this year. On the read length front, Illumina remains the clear winner – 2×100 is in production at many genome centers, and even longer reads have been promised. Life Tech, to their credit, is pushing the SOLiD 4 platform pretty hard.

When Length Matters: 454

Roche/454 has wisely backed away from large-scale sequencing, and instead seems to be targeting applications where longer (450 bp) reads are a requirement. At AGBT, Henry Erlich (Roche) gave an interesting talk about genotyping and haplotyping human HLA regions to improve donor matching for organ transplants. Here’s a key challenge of modern medicine where sequencing can offer tangible benefits. Here at the genome center, we use 454 runs for validation and for small-scale targeted sequencing. There are many applications where relatively inexpensive long-read sequencing runs are idea; full-length cDNA sequencing, for example, comes to mind.

Complete Genomics: Sequencing as a Service

The business model of Complete Genomics seems a bit of a gamble to me. They aim to be the provider of relatively inexpensive, start-to-finish sequencing services. No technology or reagent sales for these guys. Instead, they want to take your samples and give you back the SNPs. In the coming years, they hope to build as many as 10 facilities throughout the world that provide these services. I’m a bit leery of Complete Genomics, not only because their proprietary technology lags behind others (currently it’s at 2X35 bp), but because they’ll need to do something like 10,000 genomes a year just to stay in business. I don’t think we’re ready for that.

Sequencing for the Masses: IonTorrent

Many of us were impressed by IonTorrent this year at AGBT. The incredibly low cost of their instrument ($50K) and sequencing runs ($300-500) mean that nearly any lab could write a grant around this technology. The sample prep, accuracy, and throughput are still a grey area, but if they prove to be good enough, high-throughput sequencing will suddenly be available to just about everyone.

Single Molecule Applications: Pac Bio and Oxford Nanopore

The true single-molecule sequencing platforms that are close to market are certainly getting everyone excited. In the next few years, however, it’s unlikely that Pacific Biosciences, Oxford Nanopore, mystery-Chinese-platform, or other companies will displace massively parallel sequencing. No, I think Illumina and SOLiD will remain the “work horses” for discovery, certainly at major genome centers. Where SMS technologies can excel, however, is ultra-long reads – think about PacBio’s strobe sequencing to resolve structural variation or finish assemblies – and lots of molecule-kinetics stuff that I don’t understand.

I think that 2010 will be an exciting and telling time for all of these platforms. In a year’s time, we should have results in hand from HiSeq, SOLiD4, PacBio, and even IonTorrent, and be able to distinguish between marketing claims and sequencing reality.

AGBT: PacBio Somewhat Unveiled

February 27, 2010 by Dan Koboldt

Yesterday the Pacific Biosciences commercial instrument (photo) was at last unveiled to a packed room of conference attendees. The road to this third generation sequencer’s release has been paved with nearly $300 million of investment capital since leaving a basement at Cornell University. PacBio, in addition to becoming something of a media darling, has quietly swelled to a several-hundred-employee company.

Since last year, PacBio claims to have achieved read lengths of up to 10.3 kbp, although I haven’t spoken to anyone outside the company who has seen reads that long. Even so, a few vignettes presented in the workshop told of how PacBio has been applied to influenza strain identification and detection of stuctural variants (SVs).

Strobe Sequencing in Real Time

Of particular interest is the “strobe sequencing” mode of the instrument, in which the detection laser is turned off for precise amounts of time to generate mate-pair-like reads spanning large fragments. This feature relies on the real time sequencing, which occurs at a very consistent per-base rate. In fact, it’s possible to infer sequence insertions and deletions as spikes or dips (respectively) in the time required to sequence a template of known size.

Kinetic Variation Applications

The kinetics of real-time sequencing offer an informative new dimension of information from the PacBio data. In a talk today, Eric Schadt of PacBio showed that the kinetics of sequencing vary significantly for “modified” bases, i.e. methylated residues. In a collaboration with Carrie Harwood (UW), PacBio is sequencing the genomes and transcriptomes of 132 isolates of a hydrogen-producing species of Rhodopseudomonas. It turned out that kinetic variation exists at many bases as a “mixture” of sequencing times; by mining these, they identified thousands of methylated bases that caused up to 12-fold variation in sequencing kinetics.

Burning Questions Unanswered

Personally, I was not entirely satisfied with the PacBio workshop. When it opened for questions, I asked the first: whether PacBio had improved any upon the “dark bases” that go by undetected in single molecule sequencing. The presenter — Stephen Turner of PacBio — first gave me a nice 2-minute lecture on why there are no such thing as “dark bases” on PacBio’s sequencing platform due to its inherent awesomeness (sarcasm mine). There is still a problem with “missed bases” but Turner was almost comically evasive (as Daniel MacArthur put it) in stating how often they occur. The next question concerned read lengths, a second topic on which Turner refused to provide concrete information.

Thus, I find myself cautious in my excitement about this new platform, and will reserve judgment until later this year, when the first of the golden-ticket early access partners begin generating data on their own PacBio SMRT sequencers.

AGBT: Focus on Cancer Genomics

February 26, 2010 by Dan Koboldt

As usual, the quality of the scientific presentations at this meeting has been outstanding. The weather, too, has improved at last:

p_00014

There are too many to cover (or even attend) completely, but one area of interest with a strong focus this year is cancer genomics. Yesterday during plenary sessions, Stacey Gabriel of the Broad Institute of MIT and Harvard presented sequencing of multiple myeloma, a liquid tumor affecting 50,000 people in the US. Around 5,200 gigabases of sequence was generated across 26 tumor samples and matched controls, yielding ~30x average depth per genome. Their mutation detection pipeline achieved an admirable validation rate for somatic SNVs (95%). Short indels were more challenging (~50% validated), and candidate rearrangements even more so (30-50% validated). However, their study validated ~40 somatic mutations per tumor, implicating known MM genes (NRAS, KRAS, TP53) as well as novel ones (DIS3, FAM46C).

Elliott Margulies on Melanoma

Last night, there was a concurrent session devoted to cancer genomics. Eliott Margulies (NIH/NHGRI) led the lineup with his work sequencing the tumor genome and matched normal of a melanoma patient. Using the Illumina platform (2×100 bp), his group achieved 36x and 43x haploid coverage for tumor and normal, respectively, with ~99% of the genome covered by at least one read. Much of the talk was devoted to their analysis pipeline, summarized as:

Initial alignment of Illumina reads with ELAND
Partitioning the reads into “genome” bins of several kilobases
Local realignment with cross_match in highly parallelized fashion
SNV calling with their “Most Probable Genotype” (MPG) method
Removal of variants with any evidence in the Germline, or ones in dbSNP

The 175,768 novel tumor-specific SNVs were classified as coding (807) or noncoding (174,961). Some 513 of 807 coding variants were nonsynonymous. Of these, 101 were selected for validation; 84 got validation results and 75 somatic coding mutations (89%) were confirmed. Unsurprisingly, Dr. Margulies used his group’s expertise in comparative genomics to closely examine the noncoding variants as well. His group recently annotated “Chai” regions of the human genome, which bear evidence of evolutionary constraint that suggest functional relevance. Some 10,285 of the 174,961 fell within Chai regions, and among them were ~2,000 variants predicted to dramatically alter the local structure of DNA (suggesting regulatory changes).

Sequencing Pre- and Post-Treatment Lung Cancer

Ian Bosdet of BC Cancer Agency presented some very interesting work on mutational profiling of pre- and post-treatment lung cancer tumors. His group had the opportunity to participate in a clinical trial at BCCA in which carefully-selected, treatment-naive NSCLC patients underwent a standard therapeutic program. First, each patient underwent a pre-treatment evaluation and biopsy. Next, they received erlotinib (an EGFR inhibitor) until the disease inevitably progressed. Then, another biopsy that was sent for pathology review, as well as DNA/RNA extraction for sequencing. Transcriptome sequencing yielded some interesting findings. For example, the expression of one gene (IER5L or IER5C, it’s hard to read my own handwriting) was highly expressed in smokers that did not respond to treatment. A screen of unmapped transcript reads against viral genomes revealed the presence of Epstein-Barr Virus transcripts in one tumor that was later re-classified as EBV-positive lymphadenocarcinoma (?).

Mutational profiling for three patients was obtained via exome capture (Agilent) and sequencing of normal, pre-treatment tumor, and post-treatment tumor samples. Somatic mutations in PHACTR2 were seen only in pre-treatment samples. Mutations in a few genes (PRMT10, RanBP2) were found at both times, but a few (YY1AP1, SNX9) were only present after treatment, suggesting a role for these genes in progressive disease.

« Previous Page