ABSTRACT
Elucidation of tumor-DNA virus associations in many cancer types has enhanced our knowledge of fundamental oncogenesis mechanisms and provided a basis for cancer prevention initiatives. RNA-Seq is a novel tool to comprehensively assess such associations. We interrogated RNA-Seq data from 3,775 malignant neoplasms in The Cancer Genome Atlas database for the presence of viral sequences. Viral integration sites were also detected in expressed transcripts using a novel approach. The detection capacity of RNA-Seq was compared to available clinical laboratory data. Human papillomavirus (HPV) transcripts were detected using RNA-Seq analysis in head-and-neck squamous cell carcinoma, uterine endometrioid carcinoma, and squamous cell carcinoma of the lung. Detection of HPV by RNA-Seq correlated with detection by in situ hybridization and immunohistochemistry in squamous cell carcinoma tumors of the head and neck. Hepatitis B virus and Epstein-Barr virus (EBV) were detected using RNA-Seq in hepatocellular carcinoma and gastric carcinoma tumors, respectively. Integration sites of viral genes and oncogenes were detected in cancers harboring HPV or hepatitis B virus but not in EBV-positive gastric carcinoma. Integration sites of expressed viral transcripts frequently involved known coding areas of the host genome. No DNA virus transcripts were detected in acute myeloid leukemia, cutaneous melanoma, low- and high-grade gliomas of the brain, and adenocarcinomas of the breast, colon and rectum, lung, prostate, ovary, kidney, and thyroid. In conclusion, this study provides a large-scale overview of the landscape of DNA viruses in human malignant cancers. While further validation is necessary for specific cancer types, our findings highlight the utility of RNA-Seq in detecting tumor-associated DNA viruses and identifying viral integration sites that may unravel novel mechanisms of cancer pathogenesis.
INTRODUCTION
The association between infection with DNA viruses and neoplasia is well established in a variety of cancer types (1). The oncogenic potential of DNA viruses is variable, and their role in cancer pathogenesis is mediated through various mechanisms, including, for example, mutagenic integration into the host genome and expression of oncogenic viral proteins (2, 3). The elucidation of such mechanisms has played a key role in enhancing our understanding of cancer pathogenesis even as novel aspects of DNA virus biology continue to be unraveled (4). Furthermore, since eradication of viral infections through prevention and vaccination initiatives could have a drastic epidemiologic impact on the incidence of virus-associated cancers, there has been a vigorous search for tumor-associated viruses in neoplasms where such an association has remained elusive.
One of the best understood causal relationships is between human papillomavirus (HPV) infection and squamous neoplasia of the anogenital and head-and-neck regions. HPV is a small, 50- to 55-nm-diameter, nonenveloped, double-stranded DNA virus that carries out its life cycle in either mucosal or cutaneous stratified squamous epithelia (5). The viral genome (8 kb in size) is amplified initially as extrachromosomal circular elements (episomes) but may eventually integrate into the host genome. Over 120 types of HPV have been identified, of which those capable of infecting humans are designated high risk or low risk on the basis of their association with human neoplasms and oncogenic potential. The oncoproteins E5, E6, and E7 are the primary agents responsible for initiation and progression of HPV-associated cancers, and they operate primarily by abrogating negative growth regulators and inducing genomic instability. The integration of HPV DNA into the host cell genome is considered an important step in malignant progression and is commonly identified in noninvasive and invasive carcinomas associated with high-risk types HPV16 and HPV18 (6–8). HPV integration sites, with a predilection for sites of known genomic fragility, have been found to be distributed randomly over the whole genome in one study (9), and the majority of integrated HPV genomes appear to be actively transcribed (10, 11).
Hepatocellular carcinoma (HCC) is one of the leading causes of cancer-related mortality in the world and is strongly associated with chronic hepatitis B virus (HBV) infection. The HBV virion consists of partially double-stranded DNA packaged with a core protein (HBcAg) and DNA polymerase within envelope proteins (HBsAg) (12). Integration of viral DNA into the genome of HCC cells has been demonstrated in several studies, and insertional mutagenesis has been identified as a critical step in HBV-mediated HCC pathogenesis (13, 14). Integration sites were initially thought to be distributed randomly throughout the host genome, but data supporting a more deliberate process that preferentially involves transcribed regions of critical genes have been reported.(15–19).
The Epstein-Barr virus (EBV) (also known as human herpesvirus 5; HHV5) is a DNA virus that infects over 90% of the world's population before adolescence. It has been associated with a wide variety of human malignancies of epithelial, hematolymphoid, and mesenchymal derivation (2, 20, 21). Gastric carcinoma associated with EBV appears to comprise a distinct entity that is predominant in younger male individuals (22, 23). This subset of gastric carcinoma, 8 to 10% of cases, is more prevalent in Caucasian and Hispanic patients than Asians, and it shows no association with Helicobacter pylori infection (22). In these cases, EBV appears to play a direct oncogenic role through genome-wide alteration of promoter methylation (24), microRNA (miRNA) expression, and expression of genes involved in cell motility and transformation pathways (25). The integration status of EBV in gastric carcinoma remains poorly understood.
The discovery of novel cancer-associated viruses and the genomic effects of known viruses on the cancer cell genome remains incomplete, technically complex (26), and an arena in which significant resources have been and continue to be consumed. To provide an overview of the landscape of DNA virus-tumor associations in malignant cancers, we developed a method to conduct a comprehensive survey of The Cancer Genome Atlas (TCGA) RNA-Seq database for viral transcripts and, where applicable, their integration sites within the host genome.
MATERIALS AND METHODS
The Cancer Genome Atlas data.Annotated RNA-Seq data were downloaded from the Cancer Genomics Hub (CGHub; https://cghub.ucsc.edu) and the TCGA Data Portal (https://tcga-data.nci.nih.gov/tcga/) at various time points between 18 May and 3 August 2012. Patient enrollment and utilization of data were conducted in accordance with TCGA human subjects protection and data access policies (http://cancergenome.nih.gov/PublishedContent/Files/pdfs/6.3.1_TCGA_Human_Subjects_and_Data_Access_policies_FINAL_011211.pdf).
RNA sequencing.Total RNA for each sample was converted into a template molecule library for sequencing on an Illumina HiSeq 2000 according to the Illumina TruSeq sample preparation kit protocol. Briefly, poly(A) mRNA was purified from total RNA (1 μg) using poly(T) oligonucleotide-attached magnetic beads. The mRNA then was fragmented, and the first strand of cDNA was synthesized from the cleaved RNA fragments using reverse transcriptase and random primers. Following the synthesis of the second strand of cDNA, end repair was performed on overhangs using T4 DNA polymerase and Klenow DNA polymerase, followed by adenylation of the 3′ ends and ligation of sequencing adapters to the ends of the cDNA fragments. The purified cDNA templates were enriched by 15 cycles of PCR amplification and validated using BioAnalyzer (Agilent, Santa Clara, CA) to assess size, purity, and concentration of the purified cDNA libraries.
The cDNA libraries were placed on an Illumina c-Bot for paired-end (PE) cluster generation according to the protocol outlined in the Illumina HiSeq analysis user guide (part 11251649, RevA). The template cDNA libraries (1.5 μg) were hybridized to a flow cell, amplified, linearized, and denatured to create a flow cell with single-stranded DNA (ssDNA) ready for sequencing. Each flow cell was sequenced on an Illumina HiSeq Genome Analyzer. Each sample underwent one lane of PE sequencing according to the protocol outlined in the Illumina HiSeq user guide (part 11251649, RevA). After completion of the 50-cycle PE sequencing run, bases and quality values (FASTQ files) were generated for each read with the current Illumina pipeline.
Mapping/alignment.We first performed quality checks on sequencing data using the HTSeq package (http://www-huber.embl.de/users/anders/HTSeq/doc/count.html). The raw paired-end (PE) reads in FASTQ format were then aligned to the human reference genome, GRCh37/hg19 (hg19), using MOSAIK open-source alignment software. MOSAIK works with PE reads and uses both a hashing scheme and the Smith-Waterman algorithm to produce gapped optimal alignments and to map exon junction-spanning reads with a local alignment option for RNA-Seq data.
Virus detection from RNA-Seq.Our approach (39) started by computationally subtracting human sequences, followed by generating a set of nonhuman sequences (e.g., viruses) on RNA-Seq. Once raw PE reads from RNA-Seq were aligned to the human genome reference, any read with more than half a read length mapped to the human reference genome is removed along with its paired mate in this subtraction step. Thus, a set of nonhuman sequences was generated after human sequence subtraction. In the second step, our algorithm determined whether the nonhuman sequences match any known viral sequences by searching a comprehensive database that includes all known viral sequences (Genome Information Broker for Viruses [GIB-V]; http://gib-v.genes.nig.ac.jp/) and quantified virus representation by a measure of the virus genome coverage (or overall count of mapped reads) to determine the existence of viruses in human samples. The algorithm subsequently excluded nontranscribed viral genome elements to eliminate/reduce the potential of nonsense reads or inclusion of nontranscribed viral genomic elements. The expression level of each viral transcript was measured by the normalized depth of coverage within each viral transcript. The cutoff of viral gene expression detection was empirically determined by profiling the distribution of viral gene expression levels across multiple cancer-associated viruses (e.g., HPV16, HPV33, and EBV) and multiple patient samples. Any viral expression level below the cutoff was treated as no expression.
Identification of virus integration sites.The genomes of viruses with expression levels detected in the previous steps were concatenated into a single genome, named chrVirus, with related annotation of each virus in refFlat format. A new hybrid reference genome, named hg19Virus, was built by combining hg19 and chrVirus. All PE reads without computational subtraction were again mapped to this reference (hg19Virus). If the PE reads were uniquely mapped with one end to hg19 and the other to chrVirus, the read pair was reported as a discordant read pair. All discordant reads were annotated by using the genes and viruses defined in the curated refFlat file. Our algorithm subsequently employed BLAT (27) to align the discordant reads against the hybrid reference genome (hg19Virus) and discards any discordant read pair if both discordant reads are aligned to the same gene or virus. It then clustered the remaining discordant read pairs that support the same integration (fusion) event (e.g., HBV-MLL4) and selected them as fusion candidates. A dynamic clustering procedure was implemented to accurately determine the exact fusion junction between a human gene and a virus. Specifically, the boundary for each discordant read cluster of a candidate fusion was estimated on the basis of discordant read mapping locations and orientations with fragment length distribution as a constraint of cluster size, which was measured by using the reads' genomic locations, excluding intronic sizes if mapped reads were located across adjacent exons in a candidate fusion. For the forward-aligned discordant reads in a fusion candidate, the clustering process started with the right-most read, and the genomic coordinate for the right-most read was used to define the in silico fusion junction, excluding outliers within the discordant read cluster. In order to remove outliers within a cluster, our algorithm implemented the robust “extreme studentized deviate” multiple-outlier procedure. If the outliers came from the right end of the cluster, the outliers were removed and the clustering process restarts with a new in silico fusion junction. If the outliers come from the left-most end of the cluster, the cluster size was reset with the in silico fusion junction intact by excluding the outlier reads. For the reverse-aligned discordant reads, the clustering process started with the left-most read, and the genomic coordinate for the left-most read was used to define the in silico fusion junction with the same outlier detection/removal processing step. For either side of the candidate fusion partner (gene or virus), this clustering process was performed independently. This dynamic clustering procedure accurately determined the exact fusion junction between a human gene and a virus. Meanwhile, an in silico sequence was generated using the consensus of reads within discordant read clusters for each fusion candidate to help the PCR primer design, which facilitates quick PCR validation.
RESULTS
RNA sequencing data from 3,775 malignant neoplasms in the TCGA database were interrogated for the presence and, as applicable, integration of site of viral transcripts. The overall median read number per sample was 132,183,938 (Table 1). None of the virus-positive cancers in this analysis harbored more than one type/strain of virus.
Neoplasms in The Cancer Genome Atlas Database surveyed for tumor-associated DNA viruses
Human papillomavirus detection in malignant cancers.We analyzed 239 squamous cell carcinomas of the head-and-neck region (HNSCC) available in the TCGA database. We detected HPV transcripts in 36 tumors as the following: 30 tumors with HPV16, 5 tumors with HPV33, and 1 tumor with HPV35. We also detected EBV in 1 tumor (discussed below). Among all cases with HPV transcripts, E7 was expressed in 22, E6 in 20, E1 in 17, and E4 in 8 tumors. In 24 tumors, HPV transcripts encoding key viral proteins/oncoproteins were integrated in the tumor genome, with the majority in association with known genes (Fig. 1 and 2). Tumors with HPV integration harbored the following types: 19 HPV16+, 4 HPV33+, and 1 HPV35+ (Table 2). Of the tumors with HPV integration, 18 have both E6 and E7 integration sites, 4 have only E7 integration sites, and 2 have only E6 integration sites.
Visualization of HPV16 integration breakpoints in the HPV16 genome. The frequency of integration breakpoints at different loci in the HPV16 genome is shown as a blue histogram. The scale bar indicates the number of tumors. The locations of the genes encoding HPV16 E6 (red), E7 (dark green), E1 (orange), E2 (purple), E4 (green), E5 (gray), L2 (light orange), and L1 (yellow) proteins are shown. Genomic positions are numbered.
Integration sites of HPV16 in head-and-neck squamous cell carcinoma tumors in the human genome (hg19). Chromosome numbers are shown (*, detected in two cases).
Head-and-neck squamous cell carcinoma tumors with integrated HPV transcripts
The clinicopathologic characteristics of HNSCC patients with available data (n = 35) are summarized in Tables 3 and 4. HPV status detected by our method correlated with perfect sensitivity and specificity with known clinicopathologic variables and with established methods for HPV detection (colorimetric in situ hybridization and/or p16ink4a expression) (28). For this HNSCC data set, the sensitivity for HPV16 detection was 100%, with 95% confidence intervals (CI) of 67.6 to 100%, and specificity for HPV16 detection was 100%, with 95% CI of 90.4 to 100%, as reported previously (39).
Clinicopathologic characteristics of patients with head-and-neck squamous cell carcinomaa
Continuous variables of clinicopathologic characteristics of patients with head-and-neck squamous cell carcinoma
We also detected HPV16 in two tumors of histologic types rarely associated with this virus. From a group of 219 lung squamous cell carcinoma tumors, one case harbored HPV16, where E1, E6, and E7 transcripts were highly expressed and integrated in NROB1 (also known as DAX1), a gene involved in steroidogenesis and cell cycle regulation. In this case, no E2, E4, E5, L1, or L2 transcripts were detected. Additionally, 1 of 253 endometrial carcinoma tumors harbored HPV16, where E1, E2, E4, E5, E6, and E7 transcripts were highly expressed, and no L1 or L2 transcripts were detected. The integration site of the latter could not be determined, because the TCGA endometrial carcinoma RNA-Seq data were not paired end.
Epstein-Barr virus detection in malignant cancers.We analyzed 71 cases of gastric carcinoma in the TCGA database and detected EBV transcripts in 4 (5.6%) tumors. Of the four tumors with unequivocal EBV association, all harbored transcripts encoding A73, RPMS1, BARF0, BALF3, BALF4, BALF5, LF1, LF2, and BILF1. One tumor additionally had significant levels of transcripts encoding BNLF2a/b, LMP1, and LMP2A (Table 5). Of note, unlike others who demonstrated EBNA1 expression in vitro in gastric carcinoma cell lines (29), we did not identify EBNA1 in any of the gastric carcinoma tumor samples with EBV transcripts. This might be due to biological variation between primary tumor samples and cell lines or low EBNA1 mRNA copy numbers. In the single head-and-neck squamous cell carcinoma tumor associated with EBV, the most abundant transcripts included those that encode BARF1/2, BdRF1, BMRF2, BLF1, BNLF2b, BBLF1, BMRF1, BLRF2, and BNLF2a. None of the tumors analyzed in this study had evidence of EBV integration into the host genome.
EBV transcripts detected in four gastric carcinoma tumors
Hepatitis B virus detection in malignant cancers.We analyzed 69 HCC tumors available in the TCGA database and detected HBV transcripts in 16 (23%) tumors. Eight of these patients had serologic evidence of HBV infection, and one patient was HBV (and hepatitis C virus) seronegative. Serologic data on the remaining patients were either negative or not available. Virus integration was identified in 15 (94%) of these tumors. Our data demonstrate frequent HBV integration within previously identified genes, namely, TERT (5 tumors) and MLL4 (3 tumors), suggesting that these sites are particularly susceptible to HBV insertion. None of the tumors had both MLL4 and TERT involvement. Other insertion sites were restricted to single cases (Table 6). Of the tumors with HBV integration, 11 have X protein integration sites, 8 have S protein integration sites, 6 have core/E antigen integration sites, 4 have pre-S protein integration sites, 4 have polymerase 1 integration sites, and 3 have polymerase 2 integration sites. Integration of two or more HBV genes was detected in 8 tumors, whereas in the remaining 7 tumors integration of only one HBV gene was detected. Interestingly, the latter group included the 3 tumors with MLL4 involvement and 2 with TERT involvement.
Hepatocellular carcinoma tumors with integrated HBV transcripts
We additionally identified HBV transcripts in one case among 460 clear-cell renal cell carcinoma tumors analyzed. In this tumor, as well as in one HCC tumor, we detected HBV S protein transcripts integrated into GLI2, generally considered a marker of activation of the sonic hedgehog signaling pathway, which has been shown to play a role in aggressive types of renal cell carcinoma (30–32) and in HCC (33, 34).
Malignant cancers without detectable DNA viruses.After confirming the accuracy of our tool (39) in malignancies that are known to be associated with DNA viruses, we turned our attention to screening all remaining TCGA tumors on which RNA-Seq data were available. Our analysis showed no evidence of transcribed viral elements in any of the following neoplasms: acute myeloid leukemia, breast carcinoma (ductal and lobular), colonic and rectal adenocarcinoma, lung adenocarcinoma, prostate adenocarcinoma, cutaneous melanoma, ovarian serous cystadenocarcinoma, papillary renal cell carcinoma, thyroid carcinoma, lower-grade glioma, and glioblastoma multiforme.
DISCUSSION
The ability to identify viral integration points within the host genome using RNA-Seq provides a useful tool to study DNA viruses in cancer. For instance, our findings are consistent with other reports that have demonstrated preferential HBV insertion into MLL4 and TERT in HCC (18, 19, 35). Although such a finding might result from an inherent susceptibility for insertional mutagenesis at these loci, the fact that these two genes are more recurrently involved in HCC might reflect a selection bias, whereby mutation of MLL4 or TERT confers a selective advantage that leads to tumor development. Support for the latter possibility is suggested by our data indicating that among 7 HCC tumors with HBV insertion at a single gene locus, 5 tumors (71%) had involvement of MLL4 or TERT (3 with MLL4 involvement, 2 with TERT involvement). This suggests that the threshold for malignant transformation is lower when these genes, particularly MLL4, are disrupted through HBV insertional mutagenesis.
The interplay between viral infections and the miRNA machinery is the subject of intensive investigation. In addition to the modulating effect of host miRNAs on viral gene expression, mounting evidence suggests the existence of viral mechanisms leading to reprogramming of host miRNA machinery (36). In a recent study by Wang and colleagues, the HBV X protein was demonstrated to trigger selective miRNA downregulation (37). Among the tumors with HBV insertion in this study, we identified two with involvement of miRNA genes (MIR4457 and MIR4424). The function of these miRNAs in HCC is unknown. Interestingly, both genes have X protein insertion sites, suggesting that insertional mutagenesis comprises another mechanism by HBV which modulates the host miRNA repertoire.
The pathogenetic role of DNA viruses has been well established in some cancers but remains unsettled in others. One of the main findings in our study is the absence of transcribed DNA viral transcripts in some of the most prevalent human cancers, including acute myeloid leukemia, cutaneous melanoma, low- and high-grade gliomas of the brain, and adenocarcinomas of the breast, colon and rectum, lung, prostate, ovary, kidney, and thyroid. In squamous cell carcinoma of the head and neck, our analysis achieved 100% sensitivity and specificity compared to established methods (28) for HPV detection. Furthermore, the vast majority of positive results across all cancers assessed were in line with expected outcomes (e.g., HPV in squamous cell carcinoma of the head and neck, HBV in hepatocellular carcinoma, etc.). On the basis of these results, the probability of false-negative results in this analysis is believed to be low. Indeed, Bexfield and Kellam (38) performed an elegant simulation to estimate the probability of detecting viral sequences based on the viral genome transcript sequence frequency and the number of sequence reads generated. On the basis of their model, 10 million sequence reads gives a >99.99% probability of detecting at least one viral sequence if every cell is infected and the viral genome transcript sequence frequency is 0.0001%. However, 10 million sequence reads gives a >60% probability of detecting at least one viral sequence if only 1 in 10 cells is infected. For our TCGA cohort, the median number of sequence reads is more than 100 million (Table 1). Thus, we have a >99.9% probability of detecting at least one viral sequence in the TCGA cohort even if only 1 in 10 cells is infected with 0.0001% viral genome transcript sequence frequency. However, the probabilities of detecting viral sequence decrease to approximately 84% with 20 million reads. Accordingly, and unless yet-unknown pathogenic DNA viruses are discovered in the future, the yield of future searches for DNA viruses in these types of cancers is likely to be limited.
One of the limitations of this study is the lack of available corresponding tissue samples for direct validation. For instance, further evaluation of the single case of clear cell renal cell carcinoma with HBV transcripts would have been warranted. In addition, it could be argued that RNA-Seq might not detect integration events that are biologically deleterious but do not result in expression of viral transcripts. The latter point, however, does not seem to be a common event (13–15, 18, 19, 26). However, while further validation of specific findings identified in this study is warranted, the overall outcome of the analysis highlights the utility of RNA-Seq in detecting tumor-associated DNA viruses and identifying viral integration sites that may unravel novel mechanisms of cancer pathogenesis.s
ACKNOWLEDGMENTS
This work was supported in part by TCGA grant U24CA143883, National Center for Research Resources grant UL1TR000371 (F.M.-B., J.N.W., X.S.), and The University of Texas MD Anderson Cancer Center support grant P30 CA016672. None of the funding sources had any role in study design, data collection, data analysis, data interpretation, or writing of the manuscript.
J.D.K. and X.S. conceived and designed the study and prepared the manuscript. X.S. designed and developed the bioinformatics algorithm. X.S., Y.C., H.Y., J.Z., and E.J.T. performed sequencing data analysis and statistical data analysis. N.M.T., M.D.W., F.M.-B., and L.J.M. reviewed findings for specific cancer types and also wrote and/or edited the manuscript. J.N.W. reviewed bioinformatics design and results and also wrote and/or edited the manuscript.
None of the authors has conflicts of interest pertaining to this study to declare.
FOOTNOTES
- Received 2 February 2013.
- Accepted 28 May 2013.
- Accepted manuscript posted online 5 June 2013.
- Address correspondence to Xiaoping Su, xsu1{at}mdanderson.org, or Joseph D. Khoury, jkhoury{at}mdanderson.org.
REFERENCES
- Copyright © 2013, American Society for Microbiology. All Rights Reserved.
