Abstract
The recent advantage obtained by next generation sequencing allows a depth investigation of a new “old” kind of noncoding transcript, the circular RNAs. Circular RNAs are nontranslated RNAs, typically nonpolyadenylated, with a resistance to exonucleases that gives them the ability to be more stable than the common linear RNA isoforms. We used a bioinformatic detection tool (CIRCexplorer) to research predictive circRNAs from the next generation sequenced data of five samples of ductal in situ carcinoma (DCIS) and matched adjacent invasive ductal carcinoma (IDC). Furthermore, we also investigated the circular RNAs expressed in MCF7, an invasive breast ductal carcinoma cell line. We described the genomic context of the predicted circular RNAs and we address the hypothetical possible functional roles. This study showed a perspective of a panel of predictive circRNAs identified and the function that circRNAs could exert.
1. Introduction
Annually more than a million new diagnoses worldwide for breast cancer (BC) emerge, becoming the most common female malignancy [1]. The term “breast cancer” refers to a complex disease, displaying distinctive histologic and genetic features. Dramatic changes occur in gene expression during the transformation of normal breast tissues into ductal in situ carcinoma (DCIS), much less important, and somewhat surprisingly, in the transition from DCIS to invasive ductal carcinoma (IDC). In fact, those from the same histological subtype share very similar genetic, epigenetic alterations and gene expression patterns. DCIS and matched adjacent IDC (synchronous DCIS and IDC) have remarkably similar genetic and copy number profile [2, 3]. Therefore, none of the current biomarkers, or genomic classifications, can predict if, and when, DCIS will progress to IDC. In principle, a better understanding of DCIS progression could highlight ways to arrest invasion at the preinvasive stage [4]. From coding to noncoding, high-throughput next generation sequencing (NGS) has been analyzing in depth the entire genome at base-pair resolution. The noncoding panorama is a real opportunity to reveal new biomarkers or potential targets: in this landscape long noncoding RNAs (lncRNAs) are an emerging class of key regulatory RNAs that do not code for protein and are not translated [5]. lncRNAs are crucial players in various key biological processes that include dosage compensation, genomic imprinting, chromatin organization, gene regulation, and alternative splicing [6]. lncRNA mechanism of action need still to be widely revealed [7]. More recently, the research interest is focused on circular RNAs (circRNAs), a large class of long noncoding RNAs circularized and discovered from a handful of transcribed genes more than twenty years ago [8]. Nowadays, there are different types of circRNAs [9] and are known to be abundant, conserved, and stable in cytoplasm and even in blood. At the time of the discovery these RNA molecules had generally been considered to be of low abundance and likely representing errors in splicing. CircRNAs are nontranslated RNAs, typically nonpolyadenylated, with a resistance to exonucleases that gives them the ability to be more stable than the common linear RNA isoforms, also increasing proper folding by imposing structural constraints on the RNA [10]. Their formation is due to a non-random back-splice event that involves in a covalent junction between the 3′ and 5′ ends, providing covalently closed continuous loops configuration, in opposition to the theory that for decades considering the circRNAs as the result of gene rearrangements or splicing artifacts producing rare isoforms [10, 11]. Several studies revealed that circRNAs are usually composed of 1 to 5 exons, where each exon can be up to three times longer than the average expressed exon. This perhaps suggests that the probability to circularizing for an exon is directly proportional to the length of that exon [12]. Although the formation mechanism and the cellular function have been completely understood, the function and the expressions are not clear. The majority of circRNAs often exhibit tissue/developmental-stage-specific expression, described recently for the brain [13]. CircRNAs are also deregulated in cancer cell lines [11], but not many studies have been performed yet on cancer patients. Recently, some studies have been providing evidences about the interaction between circRNA and miRNA, supporting a “sponge effect” of the circRNAs as posttranscriptional regulators of other ncRNAs, such as miRNAs [10, 14]. Based on these findings, we will define the panel of expressed circRNAs, in DCIS and IDC samples. We predicted strong circRNAs candidates, enforcing their existence with an independent false discovery polyadenylated subtraction, and then we studied the possible roles of selected circRNAs.
2. Material and Methods
2.1. Data
Publicly available RNA-seq data of total RNA from ten tandem DCIS/IDC samples from five patients afflicted with DCIS that are synchronous with IDC within the same breast (GSE66301, downloaded from the Sequence Read Archive (SRA)) (Illumina GAII, paired-end sequencing (2 × 75 nt length)). GSE66301 series reports data from 6 pairs of DCIS-IDC samples, we excluded one pair that has one sample that not respected the high media level (>88%) of the mapping quality reads percentage, showing the 49%. Reference genome human hg19 (February 2009, GRCh37), was downloaded from the UCSC genome browser (http://genome.ucsc.edu/). For CIRCexplorer we used STAR as aligner software [15]. We also used the fastq.gz samples of the MCF7 cell line. These are public and available at the ENCODE repository site (http://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/): poly(A) and nonpoly(A) selected RNA-seq data, generated with paired-end sequencing (2 × 75 nt length), after the trimming quality control steps.
2.2. Data Quantification and Statistical Analysis
For circRNAs, each junction required supported from at least two independent reads within the sample. The number of reads for each circRNAs identified was used as measure of expression. The raw read counts were normalized to sequencing depth by dividing the raw number with the number of total reads mapped to and multiply for 109 (SRPBM: reads per billion mapped reads). RPM (reads per million) was used for the mRNAs. The circBase track on genome browser was consulted to verify the annotation and if the circRNAs was showed by other previous experiments [16]. For the enrichment analysis of the potential disease involved with expression of circRNAs we used the Circ2Trait tool [17]. miRTarbase (http://miRTarBase.mbc.nctu.edu.tw/) was used to find the genes involved in EMT and targeted by the human miR-200b, miR-200c, and miR.429. Custom pipelines to identify: different and common circRNAs, and to create graphics, were performed using R software (http://www.R-project.org/).
2.3. Bioinformatic Pipeline
The raw total RNA sequencing data, after quality control procedures (trimming), was used for running the bioinformatic detection pipeline: CIRCexplorer. We used this algorithm because it was one of the mostly used and cited in many sounded scientific research; and in comparison to other four pipelines, it has shown high accuracy and a good sensitivity [18]. Briefly, CIRCexplorer was focused on the identification of junction reads from back-spliced exons and intron lariats (the terminal parts of all unmapped reads that are aligned independently to the genome) [19]. In parallel, we analyzed circular RNAs expression of the total RNA-seq of MCF7, an invasive breast ductal carcinoma cell line. This analysis revealed us the candidates circular RNAs expressed in polyadenylated (polyA+) and not polyA+ selected samples. Previous studies have reported that circRNAs were nonpolyadenylated [11] or RNase R-resistant [12], then we expected that the circRNAs discovered using poly(A)+ selection could be the results of miscellaneous amplification products, misalignment on references genomes or false positives, probably not well filtered by the quality control software. All the circRNAs, predicted in the polyA+ samples, were considered false positive and then used as criteria of exclusion. Basically, we applied a false discovery polyadenylated subtraction from the output lists of all the circRNAs candidates of each samples. Moreover, the circRNAs were included by at least two supporting reads. Following these selection criteria this procedure would gave us a list of high-confidence circRNAs. Here we listed an example of the command parameters used for “CIRCexplorer pipeline”: //Alignment STAR - -genomeDir “hg19_index” - -readFilesIn “fastq1” “fastq2” - -readFilesCommand zcat - -runThreadN 8 - -chimSegmentMin 20 - -chimJunctionOverhangMin 20 - -outFileNamePrefix $STARout_dir 1≫“log_file” 2≫“err_file” //convert chimeric to file.txt./star_parse.py Chimeric.out.junction fusion_junction.txt //search circRNAs./CIRCexplorer.py -j fusion_junction.txt -g “GENOME.fa” -r “ANNOTATION.txt”.
3. Results and Discussion
3.1. The Genomic Context of the Predicted Circular RNAs in Breast Cancer Samples
DCIS invasive progression has been characterized by investigating the circRNAs expression in tandem DCIS/IDC model using five samples patients afflicted with DCIS that were synchronous with IDC within the same breast (GSE66301) [20]. We obtained for each sample the list of circRNAs candidates from algorithms’ output and refined by the false positive subtraction (see Section 2 for details). As expected, there were differences in terms of the total amount of circRNAs identified in each sample: the CIRCexplorer found 1938 raw circRNAs with one read. The total candidate distinct circular RNAs that passed our criteria were 756 in the total of the breast cancer: 552 expressed in DCIS and 208 in IDC (Supplementary Table 1 in Supplementary Material available online at http://dx.doi.org/10.1155/2016/4503840). Characterizing the genomic context of the predicted circRNAs, we firstly described the genomic localization of this class of noncoding RNAs. In DCIS the 80% (443/552) of circular RNAs belonged to the “exonic circRNAs” class derived from back-spliced and the 20% belonged to ciRNAs (circular intronic RNA). ciRNAs were a form of intron lariats that have been circularized at the branchpoint 2′–5′ linkage; normally degraded from the 3′ end up to the branchpoint, ciRNAs were stabilized in some way not yet clearly described. IDC showed a different profile: 55% (112/204) were ciRNAs, the remaining were “exonic circRNAs”. The majority of the circular RNAs identified had a length under the 12500 nucleotides in both the classes studied (Figures 1(a) and 1(b) and Supplementary Figure 1). The circRNAs production can involve more than one exon, we found up to 16 exons involved both in DCIS and IDC. As expected, there was a significant positive Spearman’s rank correlation between the length and the number of exons involved in the circularization in DCIS (rho = 0.844, 1 tail value < 2.2e − 16) and IDC too (rho = 0.630, 1 tail value < 2.2e − 16). The chromosomes (chr) most represented in terms of circRNAs were showed in Figures 1(c) and 1(d): chr1, chr2, and chr3 in DCIS, while in IDC there were chr2, chr3, chr10, and chr19.
(a)
(b)
(c)
(d)
3.2. Common and Different circRNAs Were Expressed in Synchronous DCIS and IDC Samples
Two circular RNAs, belonging to two different genes, FBXW4 and PSMA7, were found in common between the IDC samples and MCF7 cell line. These data did not reflect the gene expression profiled studies and they may suggest that these discrepancies could be due to the different technology used, or that circRNAs would be mostly related to the environment where the cancer cells were growing. A panel of 18 circRNAs were both expressed in at least one DCIS and IDC sample (Table 1 and Supplementary Figure 2). Only one coding gene, SEC62 (also named TLOC1), showed two circular RNAs that differed in terms of number of the exons involved. These two circRNAs were located in the chromosome region 3q26.2, that it is interestingly frequent amplified in several cancers [21]. The longest in terms of nucleotides (11414 nts) was the hsa-circ-0001358 (alias hsa-circ-001803), while the second was hsa-circ-0122662, 8920 nts long. There were no significant differences in terms of expression between these two circRNAs comparing the levels of expression of DCIS and IDC classes (Wilcoxon-Mann-Whitney, value = 0.8). The SEC62 mRNA was expressed more than the circular RNAs isoforms (Figure 2(b)) and the biological trend of expression suggested a positive correlation between the mRNA and the two circular RNAs (Figures 2(c) and 2(d)), but this hypothesis was not significantly confirmed by the Spearman rank test. To figure out the possible involvements of these circRNAs as “sponges” for the miRNAs, we investigated if there were possible target binding sites in their sequences. Table 2 listed five possible miRNAs that can interact with the hsa-circ-0001358, no interactions were found for the hsa-circ-0122662. The miRNAs interacting with the circRNA were obtained using the Starbase human Pan cancer tool [22]. The alignment scores of each miRNAs were showed in Supplementary Figure 3. Of note, a decreased expression of MicroRNA-200 family (in Table 2 were listed three members: miR-200b-3p, miR-200c-3p, and miR-429-3p) was associated to migration and invasion of breast cancer cells in vitro in several reports in literature [23, 24]. Therefore, we investigated the most important genes involved in Epithelial Mesenchymal Transitions (EMT) and that were targeted by the microRNA-200 family to figure out a putative circRNA/miRNA/mRNA axis involved in DCIS to IDC progression (Supplementary Table 2). No significant correlations were found investigating the hsa-circ-0001358 and the ZEB1 (Zinc Finger E-Box Binding Homeobox 1), ZEB2 (Zinc Finger E-Box Binding Homeobox 2), and VIM (Vimentin) mRNAs levels. Intriguingly, two targets of the miRNAs showed an opposite behavior: BMI-1 (B lymphoma Mo-MLV insertion region 1 homolog) positively correlated (rho = 0.886 and values = 0.02) with the circ-0001358, instead FN1 (Fibronectin 1) showed a negative correlation (rho = 0.886 and values = 0.02). Lastly, using circ2traits, we searched the measure of likelihood of the hsa-circ-0001358, to figure out possible association with the 174 disease reported in the circ2traits database. Breast cancer was the third on twenty-four significantly associated disease ( value = 0.002, Supplementary Table 3). Interestingly, we found many circRNAs expressed only in DCIS or IDC, that could represent a difference in terms of selected dynamic expression. This data will be in depth investigated for a better understanding of the DCIS progression to invasive.
(a)
(b)
(c)
(d)
4. Conclusion
We have described a panel of predicted circRNAs investigating a pilot small cohort of samples having DCIS and IDC synchronous. We have to take into account two limitations of this study: (i) since the majority of circRNAs were expressed at very low levels, and we analyzed total RNA sequencing data without the use of RNase R, an enzyme that digests linear RNAs but preserves circRNAs [25], the number of circRNAs identified was lower respect the media of other reports in literature. Therefore only an RNase R treatment will permit to enrich circRNAs and then to further verify the existence of the predicted circular RNAs. (ii) The second weakness point was the small size of this cohort. To further address the potential effect of the predicted circRNAs, it will be fundamental enlarge the cohort and perform molecular biological experiments that might permit to validate these new findings. Altogether these data try also to shed light on the possible involvement of a circRNA/miRNA/mRNA axis, but these data are still preliminary and in need of a further validation. In the next future, this class of noncoding transcripts, still largely uncharted, will be revealed important in biomedical studies, especially in those seeking to find useful marker for diagnosis and prognosis and maybe to develop new therapeutics targets.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
Stefano Volinia was supported by Associazione Italiana Ricerca sul Cancro (IG 13585 and IG 17063), Ministero della Istruzione Università e Ricerca (PRIN 2010), and EPIGEN.
Supplementary Materials
The Supplementary File was composed by three Supplementary Tables and the legends of the Supplementary Figures.
Supplementary Table 1 showed the numbers of circRNAs predicted.
Supplementary Table 2 listed Spearman correlation rho values for each of the EMT genes studied in the cohort.
Supplementary Table 3 listed the diseases enriched for the candidate circRNA using the circ2traits database.
Supplementary Figure 1 showed the plot of the density presence of predicted circRNAs.
Supplementary Figure 2 showed the heatmap of log values of 18 expressed circRNAs.
Supplementary Figure 3 showed the alignment of each miRNAs against the candidate circRNA.