Processed pseudogenes acquired somatically during cancer development

Cancer evolves by mutation, with somatic reactivation of retrotransposons being one such mutational process. Germline retrotransposition can cause processed pseudogenes, but whether this occurs somatically has not been evaluated. Here we screen sequencing data from 660 cancer samples for somatically acquired pseudogenes. We find 42 events in 17 samples, especially non-small cell lung cancer (5/27) and colorectal cancer (2/11). Genomic features mirror those of germline LINE element retrotranspositions, with frequent target-site duplications (67%), consensus TTTTAA sites at insertion points, inverted rearrangements (21%), 5′ truncation (74%) and polyA tails (88%). Transcriptional consequences include expression of pseudogenes from UTRs or introns of target genes. In addition, a somatic pseudogene that integrated into the promoter and first exon of the tumour suppressor gene, MGA, abrogated expression from that allele. Thus, formation of processed pseudogenes represents a new class of mutation occurring during cancer development, with potentially diverse functional consequences depending on genomic context.

M utation underpins both the evolution of species and the development of cancer 1 . In the germline, repetitive DNA, such as long interspersed elements (LINE) and Alu repeats, form a sizable proportion of the human genome. LINE elements continue to propagate in the genome through their retrotransposition machinery. Essentially, a functional LINE element, when transcribed and translated by the host machinery, encodes two proteins that co-ordinate reverse transcription of the mRNA template and its integration back into the genome at a distant site from the original element. Insertion of transposable elements has considerably reshaped the human genome over evolutionary time. A role for retrotransposition as a mutational force in somatic cells has been increasingly recognized in the last few years, documented to occur during both normal neurogenesis [2][3][4][5] and cancer development [6][7][8][9][10] .
Germline pseudogenes, representing cDNA copies of mRNA transcripts, are a by-product of LINE-mediated retrotransposition 11 . This copying of DNA sequence through an mRNA intermediate creates several distinctive genomic features of pseudogenes, including the presence of polyA tails, absence of intronic sequence and target-site duplications. In the germline, processed pseudogenes influence evolution through gene duplication, novel exons, gene fusions 12 , sequestration of miRNAs 13 and antisense transcript production 14 .
Formation of processed pseudogenes has not been systematically studied in cancer. From a screen of 660 cancer samples, we find that somatically acquired pseudogenes are present in 2.6% of cancers, especially lung and colorectal cancers. We find a diverse array of transcriptional consequences, including expression from UTR and intronic insertions, as well as inactivation of transcription from exonic insertions.

Results
We developed bioinformatic methods to detect somatically acquired processed pseudogenes in massively parallel sequencing data from both targeted exome and genome-wide studies in cancer. Paired-end sequencing reads were aligned to the genome and transcriptome with a view to identifying reads that split exactly across canonical splice sites, were mapped to exons separated by more than the library insert size, or mapped between a pseudogene and its insertion site ( Supplementary Fig. 1). To define a putative pseudogene, we required that at least three exons from a single gene were represented in the tumour DNA, with at least two canonical splice junctions directly observed from either split reads or confirmatory capillary sequencing. To establish that pseudogenes were not germline, we analysed sequencing data from the matched normal DNA for the given patient and screened for identical events in other patients. We performed PCR on tumour and matched normal DNA for all predicted pseudogenes with mapped insertion sites, and excluded variants with a positive germline PCR band from this analysis ( Supplementary Fig. 2).
We identified 42 somatically acquired pseudogenes in 14 out of 629 primary cases and 3 out of 31 cell lines sequenced (Table 1,  Supplementary Data 1 and Supplementary Table 1). As a typical example, in a lung cancer we identified an insertion of all five exons of the gene FOPNL, including a portion of the 5 0 UTR, the full coding sequence and the full 3 0 UTR, into the eleventh intron of SND1 in the opposite orientation (Fig. 1a). All four canonical exon-exon junctions of FOPNL were crossed by sequencing reads in the tumour DNA, and a polyA tail of at least 50 bp was present at the 3 0 end of the pseudogene. No such evidence was seen in the matched germline DNA from this patient or any other, and PCR confirmed the 5 0 and 3 0 insertion points as somatic. The two insertion points were mapped to base-pair resolution and revealed a target-site duplication of 10 bp that included the motif TTTT at either end of the pseudogene.
We observed variations on this theme (Figs 1b and 2a, Supplementary Fig. 3). PolyA tails were usually, but not universally, present, seen in 88% (21/24) pseudogenes with insertion sites fully mapped. Target-site duplications were frequent (67% with insertion sites fully mapped; 16/24), although in five cases target-site deletions occurred. The majority (74%; 31/42) of somatic pseudogenes did not contain the full coding sequence of the gene, as the 5 0 end was frequently truncated during reverse transcription and insertion. Two somatic pseudogenes were novel isoforms of their template genes not recorded in the Ensembl database. We did not observe any base substitutions in the somatic pseudogenes, implying that the reverse transcriptase has reasonable fidelity and the template transcripts had not undergone RNA editing.
On the basis of these characteristics, somatic pseudogenes, like their germline counterparts, are probably a product of the reverse transcription and transposition capability of endogenous LINE elements acting on mature mRNA transcripts. Most compelling is the presence of target-site duplications of 8-20 bp, with the point of insertion of the polyA tail almost universally occurring within a TTTTAA or very similar motif (Fig. 2b), which is the classic signature of LINE retrotransposition 15,16 . Where we observed target-site deletions at insertion points the sequence showed less consensus, although there was still predilection for AT-rich sequences (Fig. 2b). We also found a somatic pseudogene inserted at the breakpoint of a somatically acquired genomic rearrangement ( Supplementary Fig. 4), a process that can occur in the germline 17 .
Interestingly, 21% (9/42) of somatic pseudogenes showed a single inverted rearrangement within the cDNA occurring away from the canonical splice sites of the gene (Fig. 2a,c, Supplementary Fig. 5). Such internal inverted rearrangements are observed in B8% of L1 LINE elements in the genome 18 and are thought to occur by a 'twin priming' mechanism. Here, not only does the TTTT overhang prime reverse transcription from the polyA tail, but the opposite strand, nicked 10-20 bp downstream, can fold back onto the mRNA transcript internally and prime another cDNA copy. These two partial copies of the mRNA are then resolved by non-homologous endjoining, leading to an inverted rearrangement 19 . Our data are consistent with this mechanism. The internal rearrangements lead to a shuffling of the exon order and direction in the pseudogene, but without duplication of sequence (Fig. 2c). The insertion points of the 5 0 end of the inverted pseudogenes showed 1-3 bp of microhomology, consistent with a second priming event. Microhomology of between 1 bp and 4 bp was also usually present at the internal rearrangement breakpoint, indicative of non-homologous end-joining repair.
The 629 primary cases included data from 18 different tumour types ( Supplementary Fig. 6, Supplementary Table 2). Somatic pseudogenes were most frequent in non-small cell lung cancer (19%; 5/27 patients) and colorectal cancer (18%; 2/11). The fact that such events occur in colorectal and lung cancers is consistent with the observation of high rates of somatic retrotransposition of LINE elements in these tumour types [7][8][9] .
For four patients with somatic pseudogenes, we sequenced more than one tumour sample (Fig. 2d). In one patient, we analysed four clonally related samples from two bronchial lesions, both of which evolved from carcinoma in situ to invasive cancer. We found somatic pseudogenes common to all four lesions and others unique to a single lesion (Fig. 2d). In another case of lung cancer we found four somatic pseudogenes, all of which were present in both the carcinoma in situ and the invasive tumour. Similarly, for two patients with colorectal cancer, we sequenced the primary tumour and a liver metastasis. In both cases, the somatic pseudogene was present in both primary and metastasis. Taken together, these data indicate that somatic pseudogenes can form relatively early in cancer development, before the tumour becomes invasive or metastatic, but can also occur in more advanced stages of disease.
Highly expressed transcripts were especially likely to be templates for somatic pseudogenes (Fig. 3a). Overall, 63% of genes acting as the template for somatic pseudogenes were among the top quartile of expressed genes for that tumour type 20 , a statistically significant enrichment (Po0.0001, Wilcoxon test). A similar property has been noted for germline pseudogenes 21 . This may explain why we see several examples in which a gene was recurrently retrotransposed as a somatic pseudogene. For example, we found two somatic pseudogene copies of HDAC1, one in a colorectal cancer and one in a lung cancer ( Supplementary Fig. 7a). Similarly, we see somatic pseudogenes involving two aldo-keto reductase genes, AKR1C1 and AKR1C3 ( Supplementary Fig. 7b). The aldo-keto reductase proteins activate polycyclic aromatic hydrocarbons in tobacco smoke, a key step in inducing the genotoxicity of these critical carcinogens 22,23 . Overexpression of these genes has been documented in lung cancer 24 , which may explain why they recurrently serve as templates for somatic pseudogenes in this disease.
Like any mutational process, the majority of somatic pseudogenes are likely to be passenger mutations, but a few will have functional consequences that may be oncogenic. Somatic pseudogenes could exert functional consequences through many different mechanisms 12 , including fusion gene formation, increased expression of the pseudogene, disruption of a gene at the insertion site, sequestration of miRNAs from the template gene 13 and production of antisense transcripts 14 . Among the 31 somatic pseudogenes with insertion sites identified, nine were inserted into introns and three into 3 0 UTRs. Of these, five were in the same orientation as the host gene, three in the opposite orientation and four had internal inverted rearrangements.
To assess transcriptional consequences of somatic pseudogene insertion, we performed RNA-sequencing on two primary cancers and three cell lines, as well as analysed 20 non-small lung cancers sequenced by TCGA 25 (Supplementary Table 3). Across these samples, there were 16 somatic pseudogenes, of ARTICLE which 10 were inserted in intergenic regions, three in introns, two in 3 0 UTRs and one in the first exon of the MGA gene. We found no evidence of expression from somatic pseudogenes inserted in intergenic regions, whereas one of the three intronic insertions was expressed. In this case, a partial KTN1 pseudogene inserted into the last intron of PSD3 in a primary squamous cell lung cancer. We saw a clear peak of expression arising from the PSD3 intron immediately adjacent to the insertion, with aberrantly mapping read pairs aligning to the KTN1 UTR on one end and the PSD3 intron on the other end ( Supplementary Fig. 8). Of the two insertions into 3 0 UTRs, both were expressed. In one, a KRT6A pseudogene was inserted into the 3 0 UTR of MLL, the latter being a well-known fusion gene in leukaemias ( Supplementary Fig. 9). The RNA-sequencing data show that the last 1.2 kb of the MLL 3 0 UTR is lost from the mature transcript and replaced by the somatic pseudogene, since we found paired-end reads spanning the 5 0 insertion point but not the 3 0 insertion point (Fig. 3b). The 3 0 UTR of MLL has been shown to regulate transcript levels 26 , a feedback loop that could be disrupted by such a change, although expression of an aberrant transcript does not in itself imply oncogenicity. Similarly, a KIF18A pseudogene inserted into the 3 0 UTRs of two overlapping genes on opposite strands, KIAA1967 and BIN3. Reads reporting both KIAA1967-KIF18A and BIN3-KIF18A fusion junctions were found in the RNA-sequencing data ( Supplementary Fig. 10).
We also observed that somatic pseudogene insertion could abrogate expression of a target gene at the insertion site. In lung adenocarcinoma cell line NCI-H2009, a PTPN12 pseudogene caused an 8 kb target-site deletion that removed the promoter and first exon of MGA (Fig. 3c)  expression is derived from the intact MGA allele, and loss of the promoter and exon 1 has eliminated expression from the disrupted allele. MGA encodes a MAX-interacting protein 27 , with focal deletions and inactivating mutations in lymphoid malignancies 28,29 . From a compendium of 7,651 exome sequences 30 , we previously found that MGA is a likely tumour suppressor gene on the basis of a statistically significant excess of nonsense mutations (qo10 À 6 ) especially in lung adenocarcinoma 31 (Supplementary Fig. 9b).
The diversity, complexity and iniquity of mutational processes operative during the development of cancer have been laid bare by whole-genome sequencing, and here we describe another novel mechanism of somatic mutation. There has been much recent interest in how retrotransposition of repeat elements in somatic cells reshapes the genome during normal brain development and during cancer development 4,8 . The formation of pseudogenes in somatic cells represents a companion mutational process, with considerable flexibility in potential mechanisms to alter a cell's transcriptional activity.   Pseudogene detection. Data were aligned to both the reference genome (GRCh37) and the reference Ensembl transcriptome using BWA 39 and the alignment coordinates from the transcriptome mapping were converted back to genome space. Owing to the different characteristics of exome versus genome data, the analyses of these alignments were optimized to deal with the different data types (Supplementary Fig. 1). It is important to note that, regardless of the analysis method, exome versus genome data are likely to have differential sensitivities for detecting pseudogenes. Targeted exome data includes less than 2% (50 Mb) of the human genome and provides high depth over exons but does not contain the intronic and intergenic regions where the majority of genomic structural rearrangements lie. Analysis of genomic alignments through our standard structural variant algorithm 35,40 was therefore sufficient for pseudogene detection in these data, as it provided high confidence calls for exon-exon junctions with few calls representing other types of genomic rearrangement. Candidate pseudogenes were required to have at least two apparent deletion events in the same gene, that is, involvement of at least three exons, with breakpoints separated by a distance of at least 500 bp but less than 50 kb, a size range that includes the majority of introns of the human genome 41   intervening intron(s) aligned with a deletion annotated in the Database of Genomic Variants, as this indicates a potential germline pseudogene. Transcriptome alignments were used to identify reads that split across exon-exon boundaries by visual inspection in IGV (Integrative Genomics Viewer).
Due to the large size of the genome, the presence of repetitive elements, which tend to be more prevalent in non-coding sequences, and the presence of multiple genomic structural rearrangements in many tumour types, genome data required a custom pipeline to be developed (Perl scripts available on request). After alignment to the transcriptome and conversion of the read mapping coordinates back to genome space, read pairs that are derived from pseudogenes are characterized by a large insert size, due to either reads mapping to adjacent exons, or a split in the CIGAR string where two halves of a read have aligned across a splice junction in the transcriptome. Read pairs were therefore filtered on these characteristics. Retained pairs were required to have an insert size of between 300 bp and 50 kb, two or fewer splits in the CIGAR string with the first and last match being at least 10 bp, and a mapping quality of at least 10. Read pairs passing these criteria were annotated against gene positions at both the 5 0 and 3 0 end of both read 1 and read 2. Read pairs where all four coordinates mapped within the same protein coding gene in Ensembl were retained. The number of supporting read pairs per gene for the tumour, matched normal and an unmatched normal panel of 23 low-depth genomes were calculated and genes with no supporting reads in the normal/normal panel but more than three supporting reads in the tumour were validated by visual inspection in IGV.
Insertion site mapping. Somatic pseudogenes were best distinguished from contamination of the genomic DNA library by either RNA-sequencing libraries or plasmids containing cDNA expression constructs by mapping the insertion sites of the pseudogene to base-pair resolution. Mappings to the reference genome were used to identify insertion sites using reads mapping between the pseudogene and elsewhere in the genome and/or reads that aligned to one side of the insertion point with soft clipping. Soft-clipped reads were realigned using BLAT, to obtain exact breakpoint coordinates. However, the feasibility of mapping insertion points depended on the data type and the position of the insertion.
In whole-genome shotgun sequencing, the insertion sites were identified for 73% (16/22) of pseudogenes whereas 86% of insertion sites could be mapped in low-depth genomes (Supplementary Data 1). The lower rate of mapping in spite of the higher coverage in whole-genome shotgun sequencing may reflect pseudogene insertions into repetitive sequences in PD7354, which would be unmappable given the short read lengths characteristic of Illumina sequencing.
Alignments to either the genome or transcriptome rarely provided insertion sites for exome data. The majority of somatic pseudogenes consist of the 3 0 UTR and varying extents of 5 0 truncation. As UTRs are not included in the bait design for exome pull-down, insertion junction sequences are commonly not represented in the bam files of these data. The 5 0 end insertion sites were mapped in 29% (2/7) of somatic pseudogenes, with 40% of the remaining pseudogenes being full-length and therefore including a 5 0 UTR. The 3 0 insertion site could only be mapped in 14% (1/7) of pseudogenes identified in exome data. To assess the relative sensitivity of genome versus exome data for the detection of pseudogenes and capacity to map the insertion site, we sequenced both the exome and whole genome of three cell lines containing somatic pseudogenes to typical depth (Supplementary Table 4). The number of read pairs supporting the presence of a pseudogene through exon-to-exon mappings and the number of discordant read pairs reporting the insertion site for matched exome and genome data were compared. Overall, exome data shows greater sensitivity for splice junctions, whereas genome data are more likely to include the insertion junctions.
Statistical analysis of recurrent point mutations. To assess whether any of the insertions disrupted tumour suppressor genes at the insertion site, we analysed point mutation calls from 7,651 exomes available from previous publications, TCGA, ICGC and in-house for evidence of a higher than expected rate of inactivating mutation, using an adaptation of the method described in Greenman et al., 42 . This analysis has been described in detail elsewhere 31 .
PCR validation and capillary sequencing. The somatic status of pseudogenes was confirmed by designing primer pairs between the pseudogene and insertion site (where known) and performing PCR on both tumour and matched normal genomic DNA. DNA for both tumour and matched normal samples was available for 12/16 somatic pseudogene insertion sites predicted from whole-genome (430 Â ) data, all 10 insertion site predictions from low-depth data, and all 3 mapped exome insertions. A total of 92% (11/12) whole-genome predictions validated, with the remaining 8% (1/12) not producing product in either tumour or normal, consistent with PCR/primer failure. All (3/3) exome predictions had confirmed somatic insertion points. Unsurprisingly, low-depth genomes, with no or low coverage sequencing of matching normals, gave a high rate of germline pseudogenes. Overall, 40% (4/10) PCRs showed product in both tumour and normal; 50% (5/10) were somatic and the remaining 10% (1/10) failed to produce product in either tumour or normal. The germline pseudogenes have not been included in this manuscript.
RNA-sequencing. RNA-sequencing was performed on the three cell lines in which somatically/in vitro acquired pseudogenes were identified and two of the primary samples in which somatic pseudogenes were identified.
Total RNA was extracted using Trizol and sequencing libraries prepared by standard Illumina RNA-sequencing of polyA-selected RNA. A 75-bp paired-end sequencing was used and between 24 and 31 Gbp of data were generated per sample. TopHat (v 1.3.3) was used to map reads to the genome. The quality of the RNA-Seq data was assessed using a number of metrics, including absence of 3 0 bias and low amounts of ribosomal RNA.
We confirmed the positions, insertion sites and expression of somatic pseudogenes in RNA-Seq using two algorithms (Shlien A., unpublished). First, we looked for discordantly mapped read pairs where one end mapped to the pseudogene and the other near the integration site. We then evaluated grouping of pairs by the consistency of the orientation in which the reads mapped, their position and their overlap with regions of high homology (multi-mapping). Second, we looked for reads that could be partially mapped to the pseudogene and partially mapped to the insertion site (split reads). To do so, we mapped all the RNA-seq data to a transcriptome database containing all known exon-exon junctions. We then shattered and remapped all of the remaining reads into k-mers (13 bp) to an indexed version of the human genome. Mapped fragments were extended one base at a time so long as they maintained a single mapping position. In this way we resolved the breakpoints of the expressed somatic pseudogenes. Note that in the absence of fusion transcript formation, re-expression of the pseudogene from its insertion site was indistinguishable from expression from the template copy.