Detection and benchmarking of somatic mutations in cancer genomes using RNA-seq data

Methods overview, sample preparation, data origin and databases used

We developed a new pipeline to detect somatic mutations in RNA-seq data, combining RNA-seq alignment using a STAR 2-pass procedure with somatic mutation detection using MuTect2 for variant calling (Cibulskis et al., 2013). Variants from RNA-seq and WES were compared, first, on a pair of RNA-seq/WES from a GBM tumor that had already been analyzed in our laboratory (Hall et al., 2018) and then on a set of nine pairs of RNA-seq and WES data from GBM tumors analyzed by the Cancer Genome Atlas (TCGA) (Brennan et al., 2013). We compared and evaluated RNA-seq and WES mutations in four steps. First, we estimated the proportion of germline or somatic mutations by comparison of identified variants to the dbSNP database (Kitts et al., 2013) which catalogs known germline variants, and the Catalogue Of Somatic Mutations in Cancer or COSMIC database (Forbes et al., 2015), respectively. The use of these databases allowed us to evaluate whether a variant was a germline (included in dbSNP but not in COSMIC) or a somatic mutation (included in COSMIC but not in dbSNP). Second, somatic mutations detected in RNA-seq-only data were consolidated to highlight mutations present in multiple tumor samples. Third, their functional impact on proteins was evaluated by using two scoring systems: SIFT and functional analysis through hidden markov models (FATHMM) with cancer-weights (FATHMMcw) (Ng & Henikoff, 2003; Shihab et al., 2013b). Fourth, we focused on mutations affecting a set of 29 genes already shown to be implicated in GBM by a previous TCGA study (Cancer Genome Atlas Research Network, 2008). Mutations falling into coding regions of these 29 genes and showing high likelihood of altered protein function were assumed to be the best GBM-related mutations and potential cancer-drivers. Finally, we repeated this analysis on an independent validation dataset consisting of 15 pairs of RNA-seq and WES data from TCGA.

We generated paired RNA-seq and WES data from one GBM tumor (SD01) collected at St. David’s Medical Center (Austin, TX, USA) after written informed consent, in a study approved by the Institutional Review Boards of St. David’s Medical Center and of the University of Texas at Austin (approval numbers AMIRB 10-5-03 and 2012-01-0040). For WES and RNA-seq, we used the exome capture kit NimbleGen SeqCap EZ (Roche, Pleasanton, CA, USA) and the NEBNext small RNA kit (NEB, Ipswich, MA, USA), respectively. Sequencing was carried out at the NGS Core Facility of the MD Anderson Cancer Center Science Park on an Illumina HiSeq 2500. Data is available in dbGaP (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs001389.v1.p1). For GBM data from TCGA, BAM files resulting from alignment were downloaded from the Genomic Data Commons data portal and used directly in the subsequent analysis pipeline since they were already aligned with STAR. To evaluate variants, two databases were used: dbSNP (Kitts et al., 2013) with the b147 build on the GRCh38 reference (37 × 10⁶ variants), and the COSMIC database v78 (Forbes et al., 2015), which contains 3.3 × 10⁶ known somatic variants. We carried out all analyses using the GRCh38 primary assembly reference acquired from GENCODE (Harrow et al., 2012). ANNOVAR (v.2016Feb01) (Wang, Li & Hakonarson, 2010) was used to annotate variants relative to RefSeq annotations (release 73) (O’Leary et al., 2016).

A pipeline to detect variants from RNA-seq data with STAR 2-pass and GATK MuTect2 and distinguish GBM-related mutations

The general pipeline used is shown in Fig. 1, with slight differences between samples (SD01 and TCGA) or techniques (RNA-seq and WES) as depicted in Fig. 1B. The workflow was adapted from GATK best practices for variant calling (Van der Auwera, 2014; Van der Auwera et al., 2013) but using MuTect2 for variant calling. The process first involved trimming the adapters with cutadapt (v1.10) (Martin, 2011) from fastq files, removing sequences that were shorter than 36 bases after trimming, and removing rRNA and tRNA sequences by aligning with BWA (v0.7.12-r1039) (Li & Durbin, 2009) to a reference built with known rRNA/tRNA. Filtered reads were then aligned with STAR aligner (v2.4.2a) using a 2-pass procedure (Dobin & Gingeras, 2015). Before variant calling, aligned reads in BAM format were sorted, duplicate reads were flagged (MarkDuplicates, Picard v2.5.0), the base scores recalibrated (BaseRecalibrator, GATK v3.6) and RNA-seq reads were split into exons (SplitNCigarReads, GATK v3.6). Variant calling was done with MuTect2 in tumor versus normal mode as described below. Variants recovered in VCF files were then separated into RNA-seq-only, Intersection and WES-only. ANNOVAR (v.2016Feb01) (Wang, Li & Hakonarson, 2010) was used to annotate variants relative to RefSeq annotations (release 73) (O’Leary et al., 2016). SIFT score/prediction (v2.3) (Ng & Henikoff, 2003), and FATHMM score/prediction with cancer weights (v2.3) (Shihab et al., 2013a, 2013b) were used to evaluate the functional impact of non-synonymous SNVs and frameshift indels. Finally, a set of 29 genes known to be related to GBM (Cancer Genome Atlas Research Network, 2008) was used to evaluate GBM-related mutations in specific pathways.

Variant calling using MuTect2 from genome analysis toolkit

MuTect2 infers genotypes with two log-odd ratios (Cibulskis et al., 2013) which score the confidence that a mutation is present in the tumor sample (TLOD score) and is absent from the matched-normal sample (NLOD score). The thresholds used by MuTect2 to consider a variant as being real and somatic (leading to the annotation “PASS”) are by default TLOD > 6.3 and NLOD > 2.2. For dbSNP variants, a higher NLOD threshold of 5.5 is used, except if the variant is also present in the COSMIC database.

Building a panel of normals for variant calling with MuTect2

The creation of a Panel of Normals (PoN) is an optional step that improves variant calling by filtering out method-specific artefacts, by doing variant calling (MuTect2) on a set of normal samples (Fig. 1A). The samples for the PoN should ideally be obtained through protocols and data processing steps closely matched to the tumor sample. For this reason, two PoN were built, one with RNA-seq data from normal samples and another with WES data from normal samples. Then, variants identified by MuTect2 in at least two normal samples were compiled together into one PoN VCF file. Although using 30 normal samples is recommended by GATK, we used only 12 normal samples as they were matched to the 12 GBM tumor samples from TCGA.

MuTect2 filters

Based on the TLOD score, MuTect2 will reject a variant when a specific TLOD > 6.3 threshold is not reached, suggesting insufficient evidence of its presence in the tumor sample (t_lod_fstar filter). homologous_mapping_event is a filter that detects homologous sequences and filters out variants falling into sequences that have three or more events observed in the tumor. clustered_events is a filter for clustered artifacts. str_contraction filters out variants from short tandem repeat regions. alt_allele_in_normal filters out variants if enough evidence is shown of its presence in the normal sample (NLOD threshold > 2.0). multi_event_alt_allele_in_normal filters out a variant when multiple events are detected at the same position in the matched-normal sample. germline_risk filters out variants that show sufficient evidence of being germline based on dbSNP, COSMIC and the matched-normal sample (NLOD value). panel_of_normals filters out variants present in at least two samples of the panel of normals.

RefSeq annotations with ANNOVAR

ANNOVAR (v.2016Feb01) (Wang, Li & Hakonarson, 2010) was used to annotate the variants in the VCF file with RefSeq Genes annotations (release 73 with reference GRCh38) (O’Leary et al., 2016) and SIFT scores/predictions (v2.3) (Ng & Henikoff, 2003). RefSeq gives the closest gene name, or the two closest genes whenever a variant falls within intergenic regions. RefSeq also gives information about the type of mutation and the eventual amino acid change, whenever a variant falls in a coding region. For effects on alternative splicing, RefSeq gives a list of all possible transcripts.

Scoring non-synonymous SNVs and indels with SIFT score

One way to assess the functional impact of an amino-acid (AA) change is to use SIFT (Ng & Henikoff, 2003), which uses homologous sequence comparison. SIFT (v2.3) gives a score based on the frequency at which an AA appears at a specific location in functionally related protein sequences. The AA change is given a predicted score: Tolerated (p > 0.05) or Deleterious (p < 0.05). Low scores typically occur in highly conserved regions that tend to be intolerant to most substitutions. On the contrary, unconserved regions tend to be more tolerant to AA changes. SIFT indel has been developed for scoring frameshifting indels (Hu & Ng, 2013), which relies on a different algorithm based on a machine learning model. It gives a prediction of damaging or neutral along with a confidence score.

Scoring non-synonymous SNVs and indels with FATHMM cancer-weighted scores

Functional analysis through hidden markov models (FATHMM v2.3) also uses homologous protein sequences to find the probability of an amino acid substitution at a given position. The algorithm relies on Hidden Markov models to compute probabilities, its final scores being a ratio between the probability of the wild-type and the mutant AA. The version used here (Shihab et al., 2013b) also incorporates cancer weights (FATHMMcw), the frequency of cancer-associated variants from the CanProVar database and wild type weights, the frequency of neutral polymorphisms from UniRef database falling in the same protein region as the variants. The final score is an indication whether an AA substitution is deleterious and associated with cancer (prediction CANCER given for score < −0.75) or neutral (prediction PASSENGER given for score > −0.75). FATHMM for indels (Shihab et al., 2015) works on indels shorter than 20 bp and emits a prediction (pathogenic or neutral) together with a confidence score (expressed in %).

Criterion to build a set of 29 genes previously shown to be altered in GBM

A set of 29 genes that were shown to be the most frequently mutated genes in GBM by a TCGA study on 91 GBM samples (Cancer Genome Atlas Research Network, 2008) was used to look for somatic mutations in GBM-related pathways. Genes selected to be part of the set were ARF, BRCA2, CBL, CDK4, CDKN2B, CDKN2C, EGFR, EP300, ERBB2, ERBB3, FGFR2, IRS1, MDM2, MDM4, MET, MSH6, NF1, P16, PDGFRB, PIK3C2B, PIK3C2G, PIK3CA, PIK3R1, PRKCZ, PTEN, RB1, SPRY2, TP53 and TSC2. These genes were shown to bear mutations in at least 2% of samples, the most altered being ARF (49%), EGFR (45%), PTEN (36%) and TP53 (35%). The “Best GBM-related mutation” (Table 1) is indicated when a mutation was included in this set of 29 genes, part of COSMIC database but not in dbSNP, resulted in an AA change and retained based on both SIFT and FATHMM scores as being functionally deleterious for protein function.

Read counts and variant features highlight differences between RNA-seq and WES variants in TCGA samples

In the majority of samples, RNA-seq showed fewer uniquely mapped reads than WES (Fig. 2A; Fig. S1). Secondary alignments and unmapped reads were generally higher in the RNA-seq data, which could be due in part to unmapped splice junction reads and mismatches in RNA-seq due to RNA editing. Adenosine to inosine is the most common form of RNA editing in humans, leading mainly to A > G and T > C base substitutions (Picardi et al., 2015), which were clearly enriched in RNA-seq compared to WES data (Fig. 2B). RNA editing site databases like DARNED (Kiran et al., 2013), RADAR (Ramaswami & Li, 2014) or Inosinome Atlas (Picardi et al., 2015) could potentially be used to filter out such variants (Piskol, Ramaswami & Li, 2013).

The proportion of variants filtered by the different MuTect2 filters are shown in Fig. 2C. MuTect2 generates two log-odd ratios, TLOD and NLOD, which can be used to infer the somatic origin of a variant (Materials and Methods). RNA-seq variants showed lower TLOD scores and slightly higher NLOD scores than WES variants. Low read counts or poor base qualities supporting the altered allele in tumor can lead to low TLOD values. Fewer RNA-seq variants met the TLOD threshold (Fig. 2C, TLODfstar). Interestingly, TLOD scores of COSMIC variants were higher than non-COSMIC variants (Fig. S2), suggesting that TLOD reflects the higher true positive rate. On the other hand, variants that also occur in the matched-normal samples could be filtered by the AltAlleleInNormal MuTect2 filter based on NLOD values. RNA-seq data from TCGA samples showed particularly low numbers of variants excluded by this filter (Fig. 2C), which could be due to coverage differences between tumor RNA-seq and matched-normal (the latter being WES data). Thus, distinct variant features given as an output by MuTect2 could be used to build a variant filtering model (Ding et al., 2012).

Variants accepted as true and somatic (PASS) by MuTect2 were higher in RNA-seq than WES for all TCGA samples (Fig. 2D). The overlap between RNA-seq and WES was small in all samples, but interestingly, the overlap increased with increasing significance of the variants. An average of only 6.60% of WES variants retained by MuTect2 (PASS) were also present in RNA-seq, while 15.9% of WES variants from coding regions and 17.2% of functional mutations were common to RNA-seq (Fig. 2D). Coverage differences between RNA-seq and WES could partially explain the phenomenon. A previous study indeed found that ∼71% of RNA-seq variants fell outside the WES capture boundaries (O’Brien et al., 2015). Moreover, they showed that a high proportion of RNA-seq-only variants were missed by WES because of their low allele fraction (AF).

As expected, the RNA-seq/WES intersection was enriched in variants from coding regions (89.7% of coding variants), since both RNA-seq and WES query exons. RNA-seq data also showed an unexpected level of intronic/intergenic variants. Intronic mRNA reads could partly come from unspliced RNA (pre-mRNA). A previous study has indeed detected many intronic mRNA variants, which could come from inefficient splicing in cancer (Sowalsky et al., 2015). On the other hand, intergenic RNA-seq variants could come from unannotated genes, non-coding RNA, retrotransposons, splicing errors (Pickrell et al., 2010) and sequencing/mapping errors.

Allele fraction and coverage are useful features to further classify variants

In theory, heterozygous mutations would show an AF around 0.5. However, somatic mutations from cancer cells are expected to appear at lower frequencies, as tumor samples are heterogeneous and not pure clones. Moreover, CNVs can lead to gain/loss of chromosomes and/or duplications of genes (Yin et al., 2009). RNA-seq-only variants showed a surprising AF distribution in that 38.2% showed an AF > 0.95 (Figs. 3A and 3B) versus only 0.50% of WES-only variants. These high AF RNA-seq-only variants showed low coverage, and the majority of them occurred in intronic/intergenic regions (85.6% of RNA-seq-only variants with AF > 0.95). Conversely, we also found a high number of RNA-seq-only variants showing AF < 0.05 (36.3% of RNA-seq variants representing 4,267 variants in nine TCGA samples). In comparison, WES-only data showed only 502 variants (22.8%) with AF < 0.05. Low AF RNA-seq-only variants mainly originated from coding regions (81.1% of RNA-seq-only variants with AF < 0.05) and often showed high coverage, which distinguished them from WES-only variants (Fig. 3A). We examined the coverage data for RNA-seq-only variants with AF < 0.05 and coverage > 500, and found only one variant (out of 2,192) that was also present in WES data, and showing only one altered read. This region of high coverage/low AF is of particular interest as it is likely to contain true somatic mutations that are missed in WES data.

COSMIC/dbSNP overlap can be used as an indicator of the somatic/germline content in TCGA samples

For each of the three classes of variants—WES-only, Intersection and RNA-seq-only—we examined the proportion in different genomic regions (Fig. 3C), potential for affecting protein function (Fig. 3D) and representation in dbSNP and COSMIC databases (Fig. 3E). The proportion of variants included in the dbSNP database is potentially an indicator of germline content among identified variants, while the overlap with the COSMIC database can serve as an indicator of somatic mutations (Fig. 3E). It must be noted that with the increasing coverage in dbSNP of variants from ever-increasing numbers of human genomes, inclusion in dbSNP cannot always rule out a somatic variant (Nadarajah et al., 2016). Nevertheless, their overlap is small, at least in the versions of the databases we used, with only 0.15% of dbSNP variants included in COSMIC (Fig. S3). Coding variants identified by both RNA-seq & WES (Intersection) showed a particularly high proportion (87.7%) included in COSMIC but not in dbSNP (COSMIC-only), which may be considered the most likely candidates for somatic mutations. A high proportion of WES-only coding variants (39.5%) and a low proportion of RNA-seq-only coding variants (3.0%) were likewise found in COSMIC-only but although the proportions were very different, both WES-only and RNA-seq-only variants contained the same order-of-magnitude COSMIC-only variants (Fig. 3E). Thus, RNA-seq-only identified 138 COSMIC-only variants from coding regions that were therefore missed by WES-only. Because COSMIC contains variants discovered mainly by WES, it is possible that many of the RNA-seq-only variants unknown to both COSMIC and dbSNP, representing 96.4% of RNA-seq-only variants from coding regions (4,402 variants in nine samples), could include many bonafide cancer somatic mutations. We therefore explored this possibility further.

Genes showing somatic mutations in multiple TCGA samples only in RNA-seq data

There were 63 genes with RNA-seq-only variants that were mutated in five or more tumors, and many genes from this group have been implicated in cancer (Fig. 4). For example, a set of three complement related genes—complement C3, α-2 macroglobulin and the complement lysis inhibitor SP-40/clusterin (CLU)—that have been implicated in various cancers including gliomas (Reis et al., 2018; Saratsis et al., 2014; Shinoura et al., 1994; Suman et al., 2016) were present in this group, and interestingly, these three proteins have been recently shown to form a network of related biomarkers in B-ALL (Cavalcante et al., 2016). One tumor contained a cluster of highly mutated genes (Fig. 4, bottom left), including SPARC and FLNA, which are associated with cell-matrix interactions and cell motility (Neuzillet et al., 2013; Xu et al., 2010), and thus possibly involved in metastasis. On the other hand, MAGED1 was linked with cell-death mechanisms (Mouri et al., 2013), which are often disrupted in cancer. One frameshift insertion was detected in the ARF1 gene located at the exact same position (G14fs) in all nine samples. This was a COSMIC-only variant with plausible AF and coverage. Despite high coverages in WES at the variant position, the insertion was never present in WES data, and since indels have been shown to be more prone to artefacts (Kroigard et al., 2016), it was not retained in Tables 1 and 2 (see below). Note that the mutational landscape presented here is distinct from the one obtained by a TCGA study on WES data (Brennan et al., 2013), which is not surprising as RNA-seq-only data is likely interrogating other regions of the genome relative to WES.

Analysis of somatic mutations found by RNA-seq without a corresponding matched normal sample

The SD01 GBM tumor sample had no corresponding matched normal to enable reliable distinction of somatic mutations from germline variants, so it presented unusual challenges. However, it is worthwhile to consider such samples because often, RNA-seq data may be available from a tumor without a corresponding matched normal sample. The total number of variant called in SD01 was much higher than the average TCGA sample (by 10.5-fold for RNA-seq and 17.7-fold for WES). SD01 had a similar number of aligned reads as the TCGA samples for both RNA-seq and WES, so the higher number of somatic variants could be in part due to the absence of matched-normal, the small panel of normals used and/or by a higher underlying mutation rate in this specific tumor. MuTect2 variant calling was carried out in tumor-only mode and only relied on TLOD values without distinction between somatic and germline variants. Many dbSNP variants were indeed observed (Fig. S4). The distribution of SD01 variants by chromosome showed a remarkably high number of variants on Chromosome 7 (Fig. S5), which could reflect amplification of Chromosome 7, a common feature in GBM (Cancer Genome Atlas Research Network, 2008). SD01 also showed a higher density of transition variants (T > C, C > T, A > G and G > A), which tend to be less deleterious, as expected for germline variants (Campbell & Eichler, 2013). Nevertheless, SD01 RNA-seq-only variants included several interesting candidate somatic mutations. One of these RNA-seq-only mutations was EGFR-A702S, found in COSMIC but not in dbSNP, and retained by both SIFT and FATHMM scores (see below). Two other frameshift insertions were also found by RNA-seq-only data in EGFR (S229fs and W477fs), with COSMIC variants found at the same AA coordinates (Table 2). Moreover, the intersection between RNA-seq and WES data in SD01 showed other interesting candidates, such as a point mutation in EGFR (A289V—retained by both SIFT/FATHMM, and present in COSMIC but not in dbSNP).

Analyzing the functional impact of somatic mutations on protein function in relation to cancer and GBM pathways

We used the algorithms FATHMM and SIFT to evaluate the potential impact of somatic variants on protein function in cancer pathways (Materials and Methods). The FATHMM and SIFT score distributions showed a significant difference only for FATHMM scores between the Intersection and WES-only (Fig. 5). Many RNA-seq-only variants scored below both FATHMM and SIFT thresholds, indicating they could be potential functional mutations. The overall proportion of variants retained by FATHMM and SIFT was higher for Intersection variants (11.8%, Fig. 5D), and slightly higher in RNA-seq-only than WES-only. Mutations present among a set of 29 hand-curated GBM-related genes were designated as the “best GBM-related mutations” (Table 1), and comprised 11 mutations. RNA-seq-only detected three of these 11 mutations, while WES-only found two and the Intersection between RNA-seq and WES found six of the 11 GBM mutations. These three RNA-seq-only mutations (EGFR-A702S, TP53-I254S and TSC2-V296fs) are thus cancer-driver candidates found only by RNA-seq and should therefore motivate the use of RNA-seq as they were missed by WES. Taken together, our results suggest that the intersection between RNA-seq data and WES yielded the highest quality GBM-related mutations in TCGA samples, for three reasons. First, variants from RNA-seq/WES intersection showed 90.5% of COSMIC-only variants (Fig. 3E), an average much higher than WES-only or RNA-seq-only data. Second, coding variants from the intersection also showed more evidence of functional alteration through their SIFT and FATHMM scores (Fig. 5). Third, 6/11 of the “best GBM-related mutations” were identified in the intersection (Table 1), even though it was the smallest group in term of variant number. Thus combining RNA-seq and WES greatly improves the confidence in certain variants discovered by WES, particularly in highly expressed genes.

New somatic/GBM-related mutations evaluation from unknown variants

Another group of findings are shown in Table 2 as being potentially undiscovered variants, as they were neither in COSMIC nor in dbSNP, but affected one of the 29 GBM-related genes and retained either by SIFT or FATHMM scores. These mutations are therefore the best candidates for being new discoveries as they implicate known GBM-related pathways. RNA-seq-only data allowed the discovery of 8/9 potentially new mutations, against only one new variant in WES-only data, which suggests that variant calling from RNA-seq has considerable potential to generate new discoveries, including in already well known pathways. For example, an RNA-seq-only variant, EGFR-L718R, showed 22 variant reads out of a total of 4,792 (AF 4.5 × 10⁻⁸). WES showed 101 reads at the same position, giving a probability of only 0.39 of at least one variant read occurring in the WES data (based on binomial probability). Interestingly, COSMIC has cataloged a different variant, L718M at the same position (Table 2).

In order to confirm the overall findings from the preceding analysis, we repeated the entire pipeline on an independent validation dataset comprising 15 GBM tumors downloaded from TCGA. The overall characteristics of the variants in this validation dataset, as well as the genomic locations, the nature of variants found in the WES, RNA-seq-only and Intersection sets, and the identities of genes showing significant RNA-seq-only variants matched well with our previous analysis (Fig. S6).