SciLifeLab Prize: you still have time to participate!

Usable human gene expression data and annotation-agnostic methods

In the last decade RNA sequencing (RNA-seq) has become the predominant assay for measuring gene expression. RNA-seq allows us to measure all expressed genes, improve gene and transcript annotation, and measure expressed sequences that otherwise are excluded in microarray studies. Typical RNA-seq analysis starts with a matrix containing the number of RNA-seq reads overlapping each gene for each sample (Love et al. 2016; Law et al. 2016). To compute such a matrix, you first map the raw data to the genome with aligners such as TopHat2 and HISAT (Kim et al. 2013; Kim, Langmead, and Salzberg 2015) and then use tools such as Rsubread and HTSeq (Y, GK, and W 2013; Anders, Pyl, and Huber 2014) to construct the read count matrix. After measuring enough samples, you can determine which genes are differentially expressed between two or more groups.

When I started my graduate studies with Jeff Leek, the most commonly used methods for differential expression analysis were DESeq (Anders and Huber 2010) and edgeR (Robinson, McCarthy, and Smyth 2010). Leek and colleagues put together a set of gene count matrices in a project called ReCount (Frazee, Langmead, and Leek 2011). ReCount allowed researchers to access several datasets without having to run the whole processing pipeline. ReCount then helped the development of new methods such as DESeq2 (Love, Huber, and Anders 2014). It was also Leek and colleagues that decided to look at RNA-seq data in a manner less dependent on potentially incomplete gene annotation. With the increase in size of RNA-seq projects, they thought it would be feasible to assess the differential expression signal at base pair resolution in an approach called differentially expressed regions (DER) finder (Frazee et al. 2014).

It was around then that I started to work with Andrew Jaffe, who had RNA-seq data from the human brain. Jaffe and colleagues hypothesized that the gene annotation for the human brain might be incomplete and were interested in applying the DER finder approach to their data. With them, I modified the DER finder approach to make it computationally feasible to analyze these human brain samples at base resolution (Jaffe et al. 2015). We indeed found that the human brain transcriptome was not fully annotated and identified intronic DERs that were present in other brain datasets, including from mouse, and showed these regions were enriched for genetic risk for brain disorders. Further computational improvements to the DER finder approach yielded the derfinder (Collado-Torres, Nellore, Frazee, et al. 2017) Bioconductor (Huber et al. 2015) package, reducing computing time from days to hours.

While I was working on derfinder, Abhinav Nellore was developing Rail-RNA (Nellore, Collado-Torres, et al. 2016) with the goal of aligning large RNA-seq datasets. We determined that we only needed coverage data instead of alignments to use derfinder and to perform gene level analyses. Coverage data files are much smaller than alignment files, which meant that we could greatly reduce our storage costs when using Rail-RNA and allowed us to think as big as possible. Initially we looked at 21,504 human RNA-seq samples and explored the landscape of exon-exon junctions (Nellore, Jaffe, et al. 2016). The reads spanning exon-exon junctions provided by Rail-RNA can be used for a second type of annotation-agnostic analysis as these reads provide information about exon skipping, alternative donor/acceptor sites and novel events. With a fairly conservative filter we determined that 18.6% of exon-exon junctions were missing in the annotation (Nellore, Jaffe, et al. 2016). This observation along with the fact that exon-exon junctions can be used to determine differential transcript usage gave strength to this second annotation-agnostic RNA-seq analysis method.

Researchers around the world have shared their raw data via the Sequence Read Archive or via large consortiums such as the Genotype-Tissue Expression study (GTEx) and the Cancer Genome Atlas (TCGA). However, it still is complicated to align and compute count matrices for these datasets. We thought that if the public data was uniformly processed that we could make this large body of human expression data more reusable. We were well positioned to carry out this project and just needed resources. Luckily, Jeff Leek, Ben Langmead and others got together and funded the implementation of Rail-RNA to analyze all the human public RNA-seq data available at the time, over 70,000 samples.

Figure 1: Overview of the data available in recount2. Reads aligned to the reference genome can be used to compute a base pair coverage curve, identify exon-exon junctions and compute gene and exon count matrices. The reads spanning exon-exon junctions (jx) are used to compute a third count matrix that includes un-annotated junctions (jx 3 and 4). Annotation-agnostic expressed regions can be determined from the coverage data.

We computed gene count matrices for these 70,000 samples and made them available as the recount2 resource that can be accessed at jhubiostatistics.shinyapps.io/recount/ as well as through the recount (Collado-Torres, Nellore, Kammers, et al. 2017) Bioconductor package, thus representing a major improvement over ReCount. To take advantage of the data as much as possible, we also released exon and exon-exon junction count matrices as well as the coverage data files required for the DER finder approach, thus enabling both annotation-agnostic methods we developed that complement traditional methods. Figure 1 illustrates the type of information that is available via recount2.

My colleagues and I believe that recount2 will fuel the development of new analytical methods and greatly increase the usability of the public human gene expression data researchers have collected over the years. In the meantime, we are working hard on improving the data in recount2 (Ellis, Collado-Torres, and Leek 2017; Wilks et al. 2017) and exploring how to leverage results from different RNA-seq analytical approaches (Jaffe et al. 2017).

Ph.D. thesis abstract

Leonardo Collado-Torres’ thesis work is centered around the development of R software packages for analyzing RNA sequencing (RNA-seq) and ChIP sequencing (ChIP-seq) high throughput genomic data. The first chapter describes the derfinder Bioconductor package which implements the DER Finder approach for identifying differentially expressed regions with RNA-seq data in an annotation-agnostic manner. The second chapter shows how derfinder can be applied to ChIP-seq data to identify differentially bounded regions. The third chapter describes the regionReport Bioconductor package for producing HTML or PDF reports from region-based genomic analyses, such as the derfinder analyses described in the previous chapters. The last thesis project Leonardo Collado-Torres carried out was the development of the recount2 resource and accompanying recount Bioconductor package using the Rail-RNA results from processing all the public human RNA-seq data at the time, which spans over 70,000 samples.

Affiliation

Lieber Institute for Brain Development, Johns Hopkins Medical Campus, Baltimore, Maryland, 21205, USA.

Education

• 2011-2016. Ph.D. in Biostatistics at Johns Hopkins Bloomberg School of Public Health, Baltimore, US. Advised by Jeff Leek and Andrew Jaffe.
• 2005-2009. Bachelor in Genomic Sciences (LCG) at the National Autonomous University of Mexico (UNAM), Cuernavaca, MX.
• 2002-2005. High school at ITESM Campus Cuernavaca, Cuernavaca, MX.

Academic and professional awards

• 2011: Awarded CONACyT Mexico scholarship for PhD studies outside Mexico.
• 2009: Summa Cum Laude for bachelor in Genomic Sciences studies at LCG-UNAM.
• 2005: Best high school average (\(\sim\) 200 students): awarded ITESM system 90% scholarship for college studies, declined to join LCG-UNAM.

Professional experience

• 2016-current. Data Scientist with Andrew Jaffe lab at LIBD, Baltimore. At LIBD, Dr. Collado-Torres is part of the Data Science team which goals include better understanding and characterizing genomics signatures in the human brain, including DNA methylation and gene expression. Leonardo helps mentor other team members, provides support for LIBD projects and is advancing his academic career as part of Andrew Jaffe’s lab.
• 2009-2011. Bioinformatician with Enrique Morett lab at IBT-UNAM, Cuernavaca, MX. Identified transcriptions start sites and transcription units in Escherichia coli and Geobacter sulfurreducens with RNA-seq data. Developed the BacterialTranscription R package.
• 2009-2011. Scientific executive at Winter Genomics in Cuernavaca, MX. Responsible for recruiting and hiring new personnel, overseeing and supervising bioinformaticians, training new employees, writing research reports and presenting them to colleagues, and organizing all scientific projects.
• 2007-2009. Undergraduate research assistant at Guillermo Davila’s lab at CCG-UNAM, Cuernavaca, MX. Determined bacteriophage ecological groups by developing a method based on codon distribution of all phage sequenced genomes. Joint work with Sur Herrera Paredes.
• 2006-2007. Undergraduate research assistant at Roberto Kolter’s lab at Harvard, Boston, US. Supervisor: Elizabeth Shank. Carried out screenings to identify bacteria that activate the production of exopolysaccharide through the activation of the gene tasA in Bacillus subtilis.