Mapping Of RNA-Seq Reads Research Articles

Identification and Analysis of RNA Editing Sites in the Chloroplast Transcripts of Aegilops tauschii L.

RNA editing is an important way to convert cytidine (C) to uridine (U) at specific sites within RNA molecules at a post-transcriptional level in the chloroplasts of higher plants. Although it has been systematically studied in many plants, little is known about RNA editing in the wheat D genome donor Aegilops tauschii L. Here, we investigated the chloroplast RNA editing of Ae. tauschii and compared it with other wheat relatives to trace the evolution of wheat. Through bioinformatics prediction, a total of 34 C-to-U editing sites were identified, 17 of which were validated using RT-PCR product sequencing. Furthermore, 60 sites were found by the RNA-Seq read mapping approach, 24 of which agreed with the prediction and six were validated experimentally. The editing sites were biased toward tCn or nCa trinucleotides and 5′-pyrimidines, which were consistent with the flanking bases of editing sites of other seed plants. Furthermore, the editing events could result in the alteration of the secondary structures and topologies of the corresponding proteins, suggesting that RNA editing might impact the function of target genes. Finally, comparative analysis found some evolutionarily conserved editing sites in wheat and two species-specific sites were also obtained. This study is the first to report on RNA editing in Aegilops tauschii L, which not only sheds light on the evolution of wheat from the point of view of RNA editing, but also lays a foundation for further studies to identify the mechanisms of C-to-U alterations.

RNAEditor: easy detection of RNA editing events and the introduction of editing islands.

RNA editing of adenosine residues to inosine ('A-to-I editing') is the most common RNA modification event detectible with RNA sequencing (RNA-seq). While not directly detectable, inosine is read by next-generation sequencers as guanine. Therefore, mapping RNA-seq reads to their corresponding reference genome can detect potential editing events by identifying 'A-to-G' conversions. However, one must exercise caution when searching for editing sites, as A-to-G conversions also arise from sequencing errors as well as mutations. To address these complexities, several algorithms and software products have been developed to accurately identify editing events. Here, we survey currently available methods to analyze RNA editing events and introduce a new easy-to-use bioinformatics tool 'RNAEditor' for the detection of RNA editing events. During the development of RNAEditor, we noticed editing often happened in clusters, which we named 'editing islands'. We developed a clustering algorithm to find editing islands and included it in RNAEditor. RNAEditor is freely available at http://rnaeditor.uni-frankfurt.de. We anticipate that RNAEditor will provide biologists with an easy-to-use tool for studying RNA editing events and the newly defined editing islands.

OR45 RNA sequencing of two lymphoblastoid cell lines reveals novel microrna transcripts of the MHC

Aim MicroRNA (miRNA) are short (∼22 bp) single stranded RNA transcripts that bind to and regulate the posttranscriptional expression of targeted mRNA transcripts. There are 12 known miRNA transcribed from within the MHC, which have been previously identified by RNA sequencing of various primary cell types. This approach to identify novel miRNA is however inherently flawed since it relies upon significant sequence concordance between the reference genome and the generated reads, which is particularly problematic for transcripts originating from within the highly polymorphic MHC. To address this gap, we sequenced short RNA transcripts (less than 50 bp) from two homozygous cell lines with fully characterized MHC haplotypes, PGF and COX, in order to identify and characterize novel microRNA (miRNA) transcripts from within the MHC. Methods Total RNA was extracted from two lymphoblastoid cell lines, PGF and COX. Size selected transcripts (less thank 50 bp) were isolated for library construction and sequenced on an Illumina HiSeq 2500. Single-end reads from PGF and COX were mapped to their respective MHC haplotype reference sequence using NovoAlign. Novel miRNA were identified using a custom computational analysis pipeline and validated by qPCR. Results In total, our analysis has identified 67 and 47 high confidence, significantly expressed, novel miRNAs from RNA-Seq data generated from PGF and COX lymphoblastoid cell lines, respectively. Of these, 25 identified novel miRNA were found to be expressed in both PGF and COX cells, while 42 and 22 novel miRNA were found to exhibit haplotype specific expression patterns in PGF and COX, respectively. Conclusions The ability to de novo assemble the MHC is critical to ensure accurate RNA-Seq read mapping within the highly variable MHC. Our results demonstrate the utility of an accurate MHC reference assembly in combination with a custom developed computational analysis pipeline to identify and characterize novel miRNA transcripts from within the MHC of two homozygous cell lines with fully characterized MHC reference sequences. Our analysis revealed haplotype specific as well as commonly expressed novel miRNA transcripts, identifying novel regulatory elements encoded within the MHC.

RapMap: a rapid, sensitive and accurate tool for mapping RNA-seq reads to transcriptomes

Motivation: The alignment of sequencing reads to a transcriptome is a common and important step in many RNA-seq analysis tasks. When aligning RNA-seq reads directly to a transcriptome (as is common in the de novo setting or when a trusted reference annotation is available), care must be taken to report the potentially large number of multi-mapping locations per read. This can pose a substantial computational burden for existing aligners, and can considerably slow downstream analysis.Results: We introduce a novel concept, quasi-mapping, and an efficient algorithm implementing this approach for mapping sequencing reads to a transcriptome. By attempting only to report the potential loci of origin of a sequencing read, and not the base-to-base alignment by which it derives from the reference, RapMap—our tool implementing quasi-mapping—is capable of mapping sequencing reads to a target transcriptome substantially faster than existing alignment tools. The algorithm we use to implement quasi-mapping uses several efficient data structures and takes advantage of the special structure of shared sequence prevalent in transcriptomes to rapidly provide highly-accurate mapping information. We demonstrate how quasi-mapping can be successfully applied to the problems of transcript-level quantification from RNA-seq reads and the clustering of contigs from de novo assembled transcriptomes into biologically meaningful groups.Availability and implementation: RapMap is implemented in C ++11 and is available as open-source software, under GPL v3, at https://github.com/COMBINE-lab/RapMap.Contact: rob.patro@cs.stonybrook.eduSupplementary information: Supplementary data are available at Bioinformatics online.

QuickRNASeq lifts large-scale RNA-seq data analyses to the next level of automation and interactive visualization

BackgroundRNA sequencing (RNA-seq), a next-generation sequencing technique for transcriptome profiling, is being increasingly used, in part driven by the decreasing cost of sequencing. Nevertheless, the analysis of the massive amounts of data generated by large-scale RNA-seq remains a challenge. Multiple algorithms pertinent to basic analyses have been developed, and there is an increasing need to automate the use of these tools so as to obtain results in an efficient and user friendly manner. Increased automation and improved visualization of the results will help make the results and findings of the analyses readily available to experimental scientists.ResultsBy combing the best open source tools developed for RNA-seq data analyses and the most advanced web 2.0 technologies, we have implemented QuickRNASeq, a pipeline for large-scale RNA-seq data analyses and visualization. The QuickRNASeq workflow consists of three main steps. In Step #1, each individual sample is processed, including mapping RNA-seq reads to a reference genome, counting the numbers of mapped reads, quality control of the aligned reads, and SNP (single nucleotide polymorphism) calling. Step #1 is computationally intensive, and can be processed in parallel. In Step #2, the results from individual samples are merged, and an integrated and interactive project report is generated. All analyses results in the report are accessible via a single HTML entry webpage. Step #3 is the data interpretation and presentation step. The rich visualization features implemented here allow end users to interactively explore the results of RNA-seq data analyses, and to gain more insights into RNA-seq datasets. In addition, we used a real world dataset to demonstrate the simplicity and efficiency of QuickRNASeq in RNA-seq data analyses and interactive visualizations. The seamless integration of automated capabilites with interactive visualizations in QuickRNASeq is not available in other published RNA-seq pipelines.ConclusionThe high degree of automation and interactivity in QuickRNASeq leads to a substantial reduction in the time and effort required prior to further downstream analyses and interpretation of the analyses findings. QuickRNASeq advances primary RNA-seq data analyses to the next level of automation, and is mature for public release and adoption.

Complete Mitochondrial Genome ofHelicoverpa zea(Lepidoptera: Noctuidae) and Expression Profiles of Mitochondrial-Encoded Genes in Early and Late Embryos

The mitochondrial genome (mitogenome) of the bollworm, Helicoverpa zea (Boddie), was assembled using paired-end nucleotide sequence reads generated with a next-generation sequencing platform. Assembly resulted in a mitogenome of 15,348 bp with greater than 17,000-fold average coverage. Organization of the H. zea mitogenome (gene order and orientation) was identical to other known lepidopteran mitogenome sequences. Compared with Helicoverpa armigera (Hübner) mitogenome, there were a few differences in the lengths of gaps between genes, but the lengths of nucleotide overlaps were essentially conserved between the two species. Nucleotide composition of the H. zea mitochondrial genome was very similar to those of the related species H. armigera and Helicoverpa punctigera Wallengren. Mapping of RNA-Seq reads obtained from 2-h eggs and 48-h embryos to protein coding genes (PCG) revealed that all H. zea PCGs were processed as single mature gene transcripts except for the bicistronic atp8 + atp6 transcript. A tRNA-like sequence predicted to form a hammer-head-like secondary structure that may play a role in transcription start and mitogenome replication was identified within the control region of the H. zea mitogenome. Similar structures were also found within the control regions of several other lepidopteran species. Expression analysis revealed significant differences in levels of expression of PCGs within each developmental stage, but the pattern of variation was similar in both developmental stages analyzed in this study. Mapping of RNA-Seq reads to PCG transcripts also identified transcription termination and polyadenylation sites that differed from the sites described in other lepidopteran species.

Abstract PR16: New tools for mapping genetic modifiers of cancer risk in the tumor microenvironment

Abstract The majority of heritable breast cancer risk is unknown. One potential source of “missing heritability” is genetic modifiers in the tumor microenvironment (TME). Although genetic modifiers in the TME have long been suspected, they have rarely been studied and are largely unknown. Here, we used two new techniques: the Consomic Xenograft Model (CXM) and species-specific RNA-seq (SSRS) to map genetic modifiers in the TME. In CXM, human breast cancer xenografts are implanted in immunodeficient consomic rat strains and tracked for tumor progression. Because the rat strains vary by one chromosome (i.e., consomic), whereas the malignant tumor cells do not differ, any observed changes in tumor phenotypes are due to genetic modifiers in the TME and can be localized to one chromosome. The SSRS method uses probabilistic mapping of RNAseq reads to a joint human and rat transcriptome to assess differential expression (DE) in malignant (human) tumor cells and the nonmalignant (rat) TME. Validation of SSRS revealed >99.4% specificity in calling human or rat reads, which was significantly better than conventional RNA-seq. Using CXM, we found that BN-derived genetic variant(s) on rat chromosome 3 significantly reduced growth of MDA-MB-231-Luc (231Luc+) tumors by 49% (P<0.05) in the SS.BN3IL2Rγ CXM strain compared with parental SSIL2Rγ. This coincided with a 3.1-fold (P<0.001) decrease in blood vascular invasion by 231Luc+ tumor cells and 7.3-fold (P<0.05) lower metastatic burden in the lungs in SS.BN3IL2Rγ compared with SSIL2Rγ, despite a paradoxical 27% (P<0.05) increase in blood vessel density (BVD) in SS.BN3IL2Rγ rats. The tumor-associated blood vessels in SS.BN3IL2Rγ rats appeared collapsed and dysfunctional, possibly explaining the decreased tumor growth and metastasis, despite increased BVD. Lymphatic vasculature and lymphogenous metastasis were completely unaffected by the SS.BN3IL2Rγ background, suggesting that the causative variant(s) on BN rat chromosome 3 are vascular cell-type specific. We used SSRS to begin identifying the TME-specific mediators on rat chromosome 3 (RNO3) that inhibit growth and hematogenous metastasis of human 231Luc+ breast cancer xenografts implanted in the SS.BN3IL2Rγ. Compared with SSIL2Rγ tumors, we identified a network of 539 DE transcripts in the TME of SS.BN3IL2Rγ rats, of which 28% (150 genes) reside on RNO3, which was significantly higher (>4-fold; P<0.001) than any other rat chromosome. Moreover, a two-sample Kolmogorov-Smirnov test revealed that the difference in distributions of adjusted p-values for RN03 versus the rest of the genome was highly significantly higher for DE genes (P=3.152e-08) or DE transcripts (P=3.441e-16). Compared with other rat chromosomes, RNO3 also had by far the highest incidence of alternative isoform usage (91% of all instances). Pathway analysis of DE genes using DAVID revealed that the two most significant GO clusters were extracellular matrix (49 genes; P<10-20) and blood vessel development (43 genes; P<10-17), which recapitulated the vascular defects observed in the SS.BN3IL2Rγ tumors. Collectively, our data demonstrate that CXM and SSRS can be used to detect genetic modifiers in the TME. Citation Format: Michael J. Flister, Alexander Stoddard, Shirng-Wern Tsaih, Angela Lemke, Jozef Lazar, Howard Jacob. New tools for mapping genetic modifiers of cancer risk in the tumor microenvironment. [abstract]. In: Proceedings of the AACR Special Conference on Translation of the Cancer Genome; Feb 7-9, 2015; San Francisco, CA. Philadelphia (PA): AACR; Cancer Res 2015;75(22 Suppl 1):Abstract nr PR16.

Abstract PR08: New tools for mapping genetic modifiers of cancer risk in the tumor microenvironment

Abstract The majority of heritable breast cancer risk is unknown. One potential source of “missing heritability” is genetic modifiers in the tumor microenvironment (TME). Although genetic modifiers in the TME have long been suspected, they have rarely been studied and are largely unknown. Here, we used two new techniques: the Consomic Xenograft Model (CXM) and species-specific RNA-seq (SSRS) to map genetic modifiers in the TME. In CXM, human breast cancer xenografts are implanted in immunodeficient consomic rat strains and tracked for tumor progression. Because the rat strains vary by one chromosome (i.e., consomic), whereas the malignant tumor cells do not differ, any observed changes in tumor phenotypes are due to genetic modifiers in the TME and can be localized to one chromosome. The SSRS method uses probabilistic mapping of RNAseq reads to a joint human and rat transcriptome to assess differential expression (DE) in malignant (human) tumor cells and the nonmalignant (rat) TME. Validation of SSRS revealed >99.4% specificity in calling human or rat reads, which was significantly better than conventional RNA-seq. Using CXM, we found that BN-derived genetic variant(s) on rat chromosome 3 significantly reduced growth of MDA-MB-231-Luc (231Luc+) tumors by 49% (P<0.05) in the SS.BN3IL2Rγ CXM strain compared with parental SSIL2Rγ. This coincided with a 3.1-fold (P<0.001) decrease in blood vascular invasion by 231Luc+ tumor cells and 7.3-fold (P<0.05) lower metastatic burden in the lungs in SS.BN3IL2Rγ compared with SSIL2Rγ, despite a paradoxical 27% (P<0.05) increase in blood vessel density (BVD) in SS.BN3IL2Rγ rats. The tumor-associated blood vessels in SS.BN3IL2Rγ rats appeared collapsed and dysfunctional, possibly explaining the decreased tumor growth and metastasis, despite increased BVD. Lymphatic vasculature and lymphogenous metastasis were completely unaffected by the SS.BN3IL2Rγ background, suggesting that the causative variant(s) on BN rat chromosome 3 are vascular cell-type specific. We used SSRS to begin identifying the TME-specific mediators on rat chromosome 3 (RNO3) that inhibit growth and hematogenous metastasis of human 231Luc+ breast cancer xenografts implanted in the SS.BN3IL2Rγ. Compared with SSIL2Rγ tumors, we identified a network of 539 DE transcripts in the TME of SS.BN3IL2Rγ rats, of which 28% (150 genes) reside on RNO3, which was significantly higher (>4-fold; P<0.001) than any other rat chromosome. Moreover, a two-sample Kolmogorov-Smirnov test revealed that the difference in distributions of adjusted p-values for RN03 versus the rest of the genome was highly significantly higher for DE genes (P=3.152e-08) or DE transcripts (P=3.441e-16). Compared with other rat chromosomes, RNO3 also had by far the highest incidence of alternative isoform usage (91% of all instances). Pathway analysis of DE genes using DAVID revealed that the two most significant GO clusters were extracellular matrix (49 genes; P<10-20) and blood vessel development (43 genes; P<10-17), which recapitulated the vascular defects observed in the SS.BN3IL2Rγ tumors. Collectively, our data demonstrate that CXM and SSRS can be used to detect genetic modifiers in the TME. Citation Format: Michael J. Flister, Alexander Stoddard, Shirng-Wern Tsaih, Angela Lemke, Jozef Lazar, Howard Jacob. New tools for mapping genetic modifiers of cancer risk in the tumor microenvironment. [abstract]. In: Proceedings of the AACR Special Conference on Computational and Systems Biology of Cancer; Feb 8-11 2015; San Francisco, CA. Philadelphia (PA): AACR; Cancer Res 2015;75(22 Suppl 2):Abstract nr PR08.

RNA-Seq Analysis of Differential Splice Junction Usage and Intron Retentions by DEXSeq.

Alternative splicing is an important biological process in the generation of multiple functional transcripts from the same genomic sequences. Differential analysis of splice junctions (SJs) and intron retentions (IRs) is helpful in the detection of alternative splicing events. In this study, we conducted differential analysis of SJs and IRs by use of DEXSeq, a Bioconductor package originally designed for differential exon usage analysis in RNA-seq data analysis. We set up an analysis pipeline including mapping of RNA-seq reads, the preparation of count tables of SJs and IRs as the input files, and the differential analysis in DEXSeq. We analyzed the public RNA-seq datasets generated from RNAi experiments on Drosophila melanogaster S2-DRSC cells to deplete RNA-binding proteins (GSE18508). The analysis confirmed previous findings on the alternative splicing of the trol and Ant2 (sesB) genes in the CG8144 (ps)-depletion experiment and identified some new alternative splicing events in other RNAi experiments. We also identified IRs that were confirmed in our SJ analysis. The proposed method used in our study can output the genomic coordinates of differentially used SJs and thus enable sequence motif search. Sequence motif search and gene function annotation analysis helped us infer the underlying mechanism in alternative splicing events. To further evaluate this method, we also applied the method to public RNA-seq data from human breast cancer (GSE45419) and the plant Arabidopsis (SRP008262). In conclusion, our study showed that DEXSeq can be adapted to differential analysis of SJs and IRs, which will facilitate the identification of alternative splicing events and provide insights into the molecular mechanisms of transcription processes and disease development.

Abstract 3217: NextGen strategies for mapping genetic modifiers in the tumor microenvironment

Abstract The majority of heritable breast cancer risk is unknown. One potential source of “missing heritability” is genetic modifiers in the tumor microenvironment (TME). Although genetic modifiers in the TME have long been suspected, they have rarely been studied and are largely unknown. Here, we used two new techniques: the Consomic Xenograft Model (CXM) and species-specific RNA-seq (SSRS) to map genetic modifiers in the TME. In CXM, human breast cancer xenografts are implanted in immunodeficient consomic rat strains and tracked for tumor progression. Because the rat strains vary by one chromosome (i.e., consomic), whereas the malignant tumor cells do not differ, any observed changes in tumor phenotypes are due to genetic modifiers in the TME and can be localized to one chromosome. The SSRS method uses probabilistic mapping of RNAseq reads to a joint human and rat transcriptome to assess differential expression (DE) in malignant (human) tumor cells and the nonmalignant (rat) TME. Validation of SSRS revealed >99.4% specificity in calling human or rat reads, which was significantly better than conventional RNA-seq. Using CXM, we found that BN-derived genetic variant(s) on rat chromosome 3 significantly reduced growth of MDA-MB-231-Luc (231Luc+) tumors by 49% (P<0.05) in the SS.BN3IL2Rγ CXM strain compared with parental SSIL2Rγ. This coincided with a 3.1-fold (P<0.001) decrease in blood vascular invasion by 231Luc+ tumor cells and 7.3-fold (P<0.05) lower metastatic burden in the lungs in SS.BN3IL2Rγ compared with SSIL2Rγ, despite a paradoxical 27% (P<0.05) increase in blood vessel density (BVD) in SS.BN3IL2Rγ rats. The tumor-associated blood vessels in SS.BN3IL2Rγ rats appeared collapsed and dysfunctional, possibly explaining the decreased tumor growth and metastasis, despite increased BVD. We used SSRS to begin identifying the TME-specific mediators on rat chromosome 3 (RNO3) that inhibit growth and hematogenous metastasis of human 231Luc+ breast cancer xenografts implanted in the SS.BN3IL2Rγ. Compared with SSIL2Rγ tumors, we identified a network of 539 DE transcripts in the TME of SS.BN3IL2Rγ rats, of which 28% (150 genes) reside on RNO3, which was significantly higher (>4-fold; P<0.001) than any other rat chromosome. Moreover, a two-sample Kolmogorov-Smirnov test revealed that the difference in distributions of adjusted p-values for RN03 versus the rest of the genome was highly significantly higher for DE genes (P = 3.152e-08) or DE transcripts (P = 3.441e-16). Compared with other rat chromosomes, RNO3 also had by far the highest incidence of alternative isoform usage (91% of all instances). Pathway analysis of DE genes using DAVID revealed that the two most significant GO clusters were extracellular matrix (49 genes; P<10−20) and blood vessel development (43 genes; P<10−17), which recapitulated the vascular defects observed in the SS.BN3IL2Rγ tumors. Collectively, our data demonstrate that CXM and SSRS can be used to detect genetic modifiers in the TME. Citation Format: Michael Flister, Angela Lemke, Michael Dwinell, Carmen Bergom, James Shull, Howard Jacob. NextGen strategies for mapping genetic modifiers in the tumor microenvironment. [abstract]. In: Proceedings of the 106th Annual Meeting of the American Association for Cancer Research; 2015 Apr 18-22; Philadelphia, PA. Philadelphia (PA): AACR; Cancer Res 2015;75(15 Suppl):Abstract nr 3217. doi:10.1158/1538-7445.AM2015-3217

An evaluation of public genomic references for mapping RNA-Seq data from Chinese hamster ovary cells.

While RNA-Seq is increasingly used as the method of choice for transcriptome analysis of mammalian cell culture processes, no universal genomic reference for mapping RNA-Seq reads from CHO cells has been reported. In previous publications, de novo transcriptomes assembled using these RNA-Seq reads were subsequently used for mapping. Potential caveats with this approach include the incomplete coverage and the non-universal nature of the de novo assemblies, leading to challenges in comparing results across studies. In order to facilitate future RNA-Seq studies in CHO cells, we performed a comprehensive evaluation of four public genomic references for CHO cells hosted by the NCBI Reference Sequence Database (RefSeq), including two annotated genomes released in 2012 and 2014 and their accompanying transcriptomes. Each genome showed significantly higher mapped rates compared to its accompanying transcriptome. Furthermore, higher mapped rates in deep intra-genic regions, especially within exons, were observed for the more recent genome release (2014) compared to the older one (2012), indicating that the 2014 genome was the preeminent reference among the four. Sequential addition of human and mouse genomes increased the total mapped rate to 87.3 and 89.7%, respectively, from 73.5% using the 2014 Chinese hamster genome alone. Thus, the sequential combination of the 2014 RefSeq Chinese hamster genome, the Ensembl human genome (h38), and the Ensembl mouse genome (m38) was suggested as the most effective strategy for mapping RNA-Seq data from CHO cells.

From Gigabyte to Kilobyte: A Bioinformatics Protocol for Mining Large RNA-Seq Transcriptomics Data.

RNA-Seq techniques generate hundreds of millions of short RNA reads using next-generation sequencing (NGS). These RNA reads can be mapped to reference genomes to investigate changes of gene expression but improved procedures for mining large RNA-Seq datasets to extract valuable biological knowledge are needed. RNAMiner—a multi-level bioinformatics protocol and pipeline—has been developed for such datasets. It includes five steps: Mapping RNA-Seq reads to a reference genome, calculating gene expression values, identifying differentially expressed genes, predicting gene functions, and constructing gene regulatory networks. To demonstrate its utility, we applied RNAMiner to datasets generated from Human, Mouse, Arabidopsis thaliana, and Drosophila melanogaster cells, and successfully identified differentially expressed genes, clustered them into cohesive functional groups, and constructed novel gene regulatory networks. The RNAMiner web service is available at http://calla.rnet.missouri.edu/rnaminer/index.html.

Transcriptome responses to Ralstonia solanacearum infection in the roots of the wild potato Solanum commersonii.

BackgroundSolanum commersonii is a wild potato species that exhibits high tolerance to both biotic and abiotic stresses and has been used as a source of genes for introgression into cultivated potato. Among the interesting features of S. commersonii is resistance to the bacterial wilt caused by Ralstonia solanacearum, one of the most devastating bacterial diseases of crops.ResultsIn this study, we used deep sequencing of S. commersonii RNA (RNA-seq) to analyze the below-ground plant transcriptional responses to R. solanacearum. While a majority of S. commersonii RNA-seq reads could be aligned to the Solanum tuberosum Group Phureja DM reference genome sequence, we identified 2,978 S. commersonii novel transcripts through assembly of unaligned S. commersonii RNA-seq reads. We also used RNA-seq to study gene expression in pathogen-challenged roots of S. commersonii accessions resistant (F118) and susceptible (F97) to the pathogen. Expression profiles obtained from read mapping to the S. tuberosum reference genome and the S. commersonii novel transcripts revealed a differential response to the pathogen in the two accessions, with 221 (F118) and 644 (F97) differentially expressed genes including S. commersonii novel transcripts in the resistant and susceptible genotypes. Interestingly, 22.6% of the F118 and 12.8% of the F97 differentially expressed genes had been previously identified as responsive to biotic stresses and half of those up-regulated in both accessions had been involved in plant pathogen responses. Finally, we compared two different methods to eliminate ribosomal RNA from the plant RNA samples in order to allow dual mapping of RNAseq reads to the host and pathogen genomes and provide insights on the advantages and limitations of each technique.ConclusionsOur work catalogues the S. commersonii transcriptome and strengthens the notion that this species encodes specific genes that are differentially expressed to respond to bacterial wilt. In addition, a high proportion of S. commersonii-specific transcripts were altered by R. solanacearum only in F118 accession, while phythormone-related genes were highly induced in F97, suggesting a markedly different response to the pathogen in the two plant accessions studied.Electronic supplementary materialThe online version of this article (doi:10.1186/s12864-015-1460-1) contains supplementary material, which is available to authorized users.

A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification.

BackgroundRNA-Seq has become increasingly popular in transcriptome profiling. One aspect of transcriptome research is to quantify the expression levels of genomic elements, such as genes, their transcripts and exons. Acquiring a transcriptome expression profile requires genomic elements to be defined in the context of the genome. Multiple human genome annotation databases exist, including RefGene (RefSeq Gene), Ensembl, and the UCSC annotation database. The impact of the choice of an annotation on estimating gene expression remains insufficiently investigated.ResultsIn this paper, we systematically characterized the impact of genome annotation choice on read mapping and transcriptome quantification by analyzing a RNA-Seq dataset generated by the Human Body Map 2.0 Project. The impact of a gene model on mapping of non-junction reads is different from junction reads. For the RNA-Seq dataset with a read length of 75 bp, on average, 95% of non-junction reads were mapped to exactly the same genomic location regardless of which gene models was used. By contrast, this percentage dropped to 53% for junction reads. In addition, about 30% of junction reads failed to align without the assistance of a gene model, while 10–15% mapped alternatively. There are 21,958 common genes among RefGene, Ensembl, and UCSC annotations. When we compared the gene quantification results in RefGene and Ensembl annotations, 20% of genes are not expressed, and thus have a zero count in both annotations. Surprisingly, identical gene quantification results were obtained for only 16.3% (about one sixth) of genes. Approximately 28.1% of genes’ expression levels differed by 5% or higher, and of those, the relative expression levels for 9.3% of genes (equivalent to 2038) differed by 50% or greater. The case studies revealed that the gene definition differences in gene models frequently result in inconsistency in gene quantification.ConclusionsWe demonstrated that the choice of a gene model has a dramatic effect on both gene quantification and differential analysis. Our research will help RNA-Seq data analysts to make an informed choice of gene model in practical RNA-Seq data analysis.Electronic supplementary materialThe online version of this article (doi:10.1186/s12864-015-1308-8) contains supplementary material, which is available to authorized users.

Allelic mapping bias in RNA-sequencing is not a major confounder in eQTL studies.

BackgroundRNA sequencing (RNA-seq) is the current gold-standard method to quantify gene expression for expression quantitative trait locus (eQTL) studies. However, a potential caveat in these studies is that RNA-seq reads carrying the non-reference allele of variant loci can have lower probability to map correctly to the reference genome, which could bias gene quantifications and cause false positive eQTL associations. In this study, we analyze the effect of this allelic mapping bias in eQTL discovery.ResultsWe simulate RNA-seq read mapping over 9.5 M common SNPs and indels, with 15.6% of variants showing biased mapping rate for reference versus non-reference reads. However, removing potentially biased RNA-seq reads from an eQTL dataset of 185 individuals has a very small effect on gene and exon quantifications and eQTL discovery. We detect only a handful of likely false positive eQTLs, and overall eQTL SNPs show no significant enrichment for high mapping bias.ConclusionOur results suggest that RNA-seq quantifications are generally robust against allelic mapping bias, and that this does not have a severe effect on eQTL discovery. Nevertheless, we provide our catalog of putatively biased loci to allow better controlling for mapping bias to obtain more accurate results in future RNA-seq studies.Electronic supplementary materialThe online version of this article (doi:10.1186/s13059-014-0467-2) contains supplementary material, which is available to authorized users.

Assessment of the impact of using a reference transcriptome in mapping short RNA-Seq reads.

RNA-Seq has become increasingly popular in transcriptome profiling. The major challenge in RNA-Seq data analysis is the accurate mapping of junction reads to their genomic origins. To detect splicing sites in short reads, many RNA-Seq aligners use reference transcriptome to inform placement of junction reads. However, no systematic evaluation has been performed to assess or quantify the benefits of incorporating reference transcriptome in mapping RNA-Seq reads. In this paper, we have studied the impact of reference transcriptome on mapping RNA-Seq reads, especially on junction ones. The same dataset were analysed with and without RefGene transcriptome, respectively. Then a Perl script was developed to analyse and compare the mapping results. It was found that about 50–55% junction reads can be mapped to the same genomic regions regardless of the usage of RefGene model. More than one-third of reads fail to be mapped without the help of a reference transcriptome. For “Alternatively” mapped reads, i.e., those reads mapped differently with and without RefGene model, the mappings without RefGene model are usually worse than their corresponding alignments with RefGene model. For junction reads that span more than two exons, it is less likely to align them correctly without the assistance of reference transcriptome. As the sequencing technology evolves, the read length is becoming longer and longer. When reads become longer, they are more likely to span multiple exons, and thus the mapping of long junction reads is actually becoming more and more challenging without the assistance of reference transcriptome. Therefore, the advantages of using reference transcriptome in the mapping demonstrated in this study are becoming more evident for longer reads. In addition, the effect of the completeness of reference transcriptome on mapping of RNA-Seq reads is discussed.

Integration of mapped RNA-Seq reads into automatic training of eukaryotic gene finding algorithm.

We present a new approach to automatic training of a eukaryotic ab initio gene finding algorithm. With the advent of Next-Generation Sequencing, automatic training has become paramount, allowing genome annotation pipelines to keep pace with the speed of genome sequencing. Earlier we developed GeneMark-ES, currently the only gene finding algorithm for eukaryotic genomes that performs automatic training in unsupervised ab initio mode. The new algorithm, GeneMark-ET augments GeneMark-ES with a novel method that integrates RNA-Seq read alignments into the self-training procedure. Use of ‘assembled’ RNA-Seq transcripts is far from trivial; significant error rate of assembly was revealed in recent assessments. We demonstrated in computational experiments that the proposed method of incorporation of ‘unassembled’ RNA-Seq reads improves the accuracy of gene prediction; particularly, for the 1.3 GB genome of Aedes aegypti the mean value of prediction Sensitivity and Specificity at the gene level increased over GeneMark-ES by 24.5%. In the current surge of genomic data when the need for accurate sequence annotation is higher than ever, GeneMark-ET will be a valuable addition to the narrow arsenal of automatic gene prediction tools.

Inducible Defenses Stay Up Late: Temporal Patterns of Immune Gene Expression inTenebrio molitor

The course of microbial infection in insects is shaped by a two-stage process of immune defense. Constitutive defenses, such as engulfment and melanization, act immediately and are followed by inducible defenses, archetypically the production of antimicrobial peptides, which eliminate or suppress the remaining microbes. By applying RNAseq across a 7-day time course, we sought to characterize the long-lasting immune response to bacterial challenge in the mealworm beetle Tenebrio molitor, a model for the biochemistry of insect immunity and persistent bacterial infection. By annotating a hybrid de novo assembly of RNAseq data, we were able to identify putative orthologs for the majority of components of the conserved insect immune system. Compared with Tribolium castaneum, the most closely related species with a reference genome sequence and a manually curated immune system annotation, the T. molitor immune gene count was lower, with lineage-specific expansions of genes encoding serine proteases and their countervailing inhibitors accounting for the majority of the deficit. Quantitative mapping of RNAseq reads to the reference assembly showed that expression of genes with predicted functions in cellular immunity, wound healing, melanization, and the production of reactive oxygen species was transiently induced immediately after immune challenge. In contrast, expression of genes encoding antimicrobial peptides or components of the Toll signaling pathway and iron sequestration response remained elevated for at least 7 days. Numerous genes involved in metabolism and nutrient storage were repressed, indicating a possible cost of immune induction. Strikingly, the expression of almost all antibacterial peptides followed the same pattern of long-lasting induction, regardless of their spectra of activity, signaling possible interactive roles in vivo.

OMIGA: Optimized Maker-Based Insect Genome Annotation

Insects are one of the largest classes of animals on Earth and constitute more than half of all living species. The i5k initiative has begun sequencing of more than 5,000 insect genomes, which should greatly help in exploring insect resource and pest control. Insect genome annotation remains challenging because many insects have high levels of heterozygosity. To improve the quality of insect genome annotation, we developed a pipeline, named Optimized Maker-Based Insect Genome Annotation (OMIGA), to predict protein-coding genes from insect genomes. We first mapped RNA-Seq reads to genomic scaffolds to determine transcribed regions using Bowtie, and the putative transcripts were assembled using Cufflink. We then selected highly reliable transcripts with intact coding sequences to train de novo gene prediction software, including Augustus. The re-trained software was used to predict genes from insect genomes. Exonerate was used to refine gene structure and to determine near exact exon/intron boundary in the genome. Finally, we used the software Maker to integrate data from RNA-Seq, de novo gene prediction, and protein alignment to produce an official gene set. The OMIGA pipeline was used to annotate the draft genome of an important insect pest, Chilo suppressalis, yielding 12,548 genes. Different strategies were compared, which demonstrated that OMIGA had the best performance. In summary, we present a comprehensive pipeline for identifying genes in insect genomes that can be widely used to improve the annotation quality in insects. OMIGA is provided at http://ento.njau.edu.cn/omiga.html .

Allelic mapping bias in RNA-sequencing is not a major confounder in eQTL studies

RNA sequencing (RNA-seq) is the current gold-standard method to quantify gene expression for expression quantitative trait locus (eQTL) studies. However, a potential caveat in these studies is that RNA-seq reads carrying the non-reference allele of variant loci can have lower probability to map correctly to the reference genome, which could bias gene quantifications and cause false positive eQTL associations. In this study, we analyze the effect of this allelic mapping bias in eQTL discovery. We simulate RNA-seq read mapping over 9.5 M common SNPs and indels, with 15.6% of variants showing biased mapping rate for reference versus non-reference reads. However, removing potentially biased RNA-seq reads from an eQTL dataset of 185 individuals has a very small effect on gene and exon quantifications and eQTL discovery. We detect only a handful of likely false positive eQTLs, and overall eQTL SNPs show no significant enrichment for high mapping bias. Our results suggest that RNA-seq quantifications are generally robust against allelic mapping bias, and that this does not have a severe effect on eQTL discovery. Nevertheless, we provide our catalog of putatively biased loci to allow better controlling for mapping bias to obtain more accurate results in future RNA-seq studies.