Canadian Bioinformatics Workshops

Canadian Bioinformatics Workshops www.bioinformatics.ca

Module #: Title of Module 2

Module 5Metatranscriptomics John Parkinson Analysis of Metagenomic Data June 24-26, 2015

Learning Objectives of Module • At the end of this module the student will have an appreciation of the opportunities and challenges of metatranscriptomics by: • Understanding the capabilities of metatranscriptomics • Gaining an appreciation of sample collection and experimental design • Learning important steps in data processing • Being able to process a simple metatranscriptomic dataset

Overview • What is metatranscriptomics and RNA-Seq? • Experimental design, sample collection and preparation • Processing of reads • Filters • Assembly • Dereplication • Functional annotation • Taxonomic annotation • Statistical analysis • Visualization

Metagenomics focuses on community makeup Metagenomics is typically focussed on community makeup, identifying the proportion of sequences (DNA or 16S rRNA) that map to different taxa • Although metagenomics can inform on the species present it has several drawbacks • Does not provide much in the way of functional context, i.e. which genes are important • Different sets of species have been shown to result in similar functional potential • Relying on 16S sequencing provides species level information, but due to HGT the gene complements of different strains of the same species can vary considerably Qin et al Nature 2010

Metatranscriptomics focuses on community activity Metatranscriptomics exploits RNA-Seq to determine which genes and pathways are being actively expressed within a community Genes involved in pathways associated with cell wall biogenesis Metatranscriptomics can reveal active functions (knowing the taxa responsible is unimportant) It can also reveal which taxa are responsible for the active functions Relative abundance and contributing taxa Relative abundance

Metatranscriptomics through RNA Seq Species D Species B Species C Species A Fragment and sequence Extract RNA RNA-Seq is the unbiased sequencing of an RNA sample to yield a digital readout of the relative expression of transcripts within a sample Typically applied to organisms with a reference (sequenced) genome, microbiome applications face a number of challenges 6 3 3 1 2 2 6 4 3 Gene A1 Gene A2 Gene A3 Gene A4 Gene B1 Gene B2 Gene B3 Gene C1 Align reads to known transcripts to obtain relative expression Gene D4 III IV

Metatranscriptomics: Challenges In a typical RNA-Seq experiment applied to a single eukaryotic organism, mRNA is isolated. After fragmentation and sequencing, reads are mapped to a reference genome using standard software such as MAQ and BWA to provide: 1) support that the transcript is expressed; 2) the relative abundance of the transcript; and 3) the presence and abundance of isoforms • For microbiome samples we have the following problems: • Lack of a polyAsignal makes it difficult to isolate bacterial mRNA and resulting in (massive) rRNA contamination • Environmental microbiome samples lack reference genomes making it difficult to map reads back to their source transcripts III IV

A typical metatranscriptomic analytical pipeline 1. Obtain RNA 3. Generate Reads 2. Prep for sequencing 7. Assemble 5. Remove rRNA reads 6. Remove host reads 4. Remove low quality 7. Identify bacterial transcript 8. Map to pathways Gene A Gene C 9. Sample comparisons Gene B Gene D

1. Sample collection and RNA extraction RNA quality deteriorates rapidly Best: Process immediately to extract RNA then store at -80 Next best: Snap freeze in liquid nitrogen and store at -80 Avoid use of RNALater – it lyses some cells and can interfere with RNA extraction kits (e.gAmbionRiboPure Bacteria Kit / LifeTechTrizol Plus) Number of biological replicates “At least 2!” Power analyses challenging Cost per sample (20 million reads) ~$300-$400

2. Preparing sample for sequencing Bacterial mRNA’s lack a polyA tail so how to remove abundant rRNA species? Once RNA has been extracted, several kits are available to remove rRNA – need 500ng-2.5ug RNA/sample Ribo-Zero (Illumina) provides reasonable success mRNA Host rRNA Adaptor / Low quality Host mRNAs can also prove challenging – can also be informative!

2. Preparing sample for sequencing Comparisons of storage and RNA treatments Temperature (4oC > -20oC) Storage Days (Fresh > 1 week) Kits (RiboZero > Mxpress) RNA later (+ > -)

3. Generating reads How many reads are “enough”? While 454 and MiSeq provide long reads useful for annotation HiSeq (or NextSeq) provide sequencing depth and offer possibility of multiplexing ~5 million mRNA reads provide 90-95% of ECs in a microbiome With kits yielding mRNA read rates of ~25%, this suggests 20 million/sample mRNA

4.-5. Read processing - filtering Trimmomatic uses a sliding window approach from the 5` end to identify low quality regions which are then trimmed from the 3` end. Reads < 36 bp are discarded To identify reads derived from mRNA bioinformatics pipelines need to be in place that remove contaminating reads: Low quality - Trimmomatic Adaptors – Trimmomatic & Cross_Match Host - BWA rRNA – BLAT / Infernal

6. Read processing - Assembly contig2 contig1 attagcggcgattttcggcgatcttatcttgatctgggcgcgtatcggtagcgtagcgattcgtagc attagcggcgattttcggcgatcttatcttgatctgggcgcgtatcggtagcgtagcgattcgtagc attagcggcgattttcggcgatcttatcttgatctgggcgcgtatcggtagcgtagcgattcgtagc Assembly improves annotation accuracy Trinity appears to provide best performance in terms of reads that can be annotated Chimera’s, misassembled contigs, can become a problem due to reads derived from orthologs from different species

7. Read processing – functional annotation Functional annotations rely on the identification of sequence similarity searches using tools such as: BWA, BLAT and BLAST While fast, BWA and BLAT rely on near perfect matches to reference genomes However sequence diversity is huge! Every time a new strain of S. galactae was sequenced, new genes are discovered that were not present in existing strains Diversity at the nucleotide level is >> protein level

7. Read processing – functional annotation One solution is to work in peptide space and use BLASTX to search protein databases - this is very time consuming and requires cloud/cluster computing Other solutions include MBLAST and USEARCH (issues over quality and cost) Even with BLAST many reads remain unannotated Can be improved with longer read length

7. Quality of BLASTX matches Typical match to a 71bp read E-Value = 39 (NOT e-39) Species looks about right though % of Read Length Summary for all matches % ID of match

7. Read processing – converting mappings to expression To normalize expression levels to account for differences in gene length, read counts are converted to Reads per kilobase of transcript mapped (RPKM) Expression is biased for gene length (longer transcripts should have more reads) to normalize, reads are converted to Reads per Kilobase of transcript per million reads mapped RPKMgeneA = 109CgeneA / NL CgeneA = number of reads mapped to geneA N = total number of reads L = length of transcript in units of Kb RPKM # 1 8 6 24 12 8 Several software tools available to do mapping and calculate normalized expression measurements across different samples including Bowtie and Cufflinks

8. Read processing – taxonomic annotation Previous studies have shown that microbiomes can vary significantly in taxonomic contributions while yielding similar functionality Knowing which species are providing which functions may not be that important On the other hand, assigning RNA reads to taxa may reveal critical functions contributed by keystone taxa, can also help in binning for assembly How do we assign taxonomic information? The human microbiome consortium 2012

8. Read processing – taxonomic annotation Alignment based methods such as BLAST and BWA can fail where we lack suitable reference genomes – particularly for short read datasets where assignments may be ambiguous Compositional methods (e.g. nt frequency, codon bias) offer alternative strategies NBC Nearest neighbours methods then try to assign a sequence to the genome with the closest distribution Here a sequences is classified into frequencies of 3-mers MetaCV RITA

AA NB AA CD AA 1NN AA GMM classifier construction 8. Read processing – Gist BWA NT NB NT CD NT 1NN NT GMM class definitions generate class distributions Gist is a computational pipeline for accurate assignment of reads to individual species Integrates several methods, but uniquely assigns different weights to methods for each genome Can also take in expected sequence distributions (e.g. based on 16S rRNA surveys) taxonomy tags taxonomy database reshape into strain-specific records (DELIN) environment description 16S counts other abundances log normalization taxonomic assignment genomes (*.ffn) class priors FragGeneScan mRNA amino acid sequences infer protein sequence Classification 1 select candidate hits for detailed processing for each read initial confidence measurements with T-test for hit specificity CalcTaxonInclusion Classification 2 final assignment scores select best hits CalcTaxonInclusion final confidence measurements (T-test for hit specificity)

8. Read processing – taxonomic annotation A key advance in Gist is finding the set of weights for each feature/method that best discriminates between different genomes – here C. diff 630 is reliant on Naïve Bayes analysis of amino acid frequencies, whereas S. agalactiae A909 relies more on a Gaussian mixture model of amino acid frequencies

8. Read processing – taxonomic annotation Gist outperforms MetaCV which features a large number of ‘unclassifiable’ NBC works very well if you have reference genomes, but due to reliance on large kmers, can be prone to genetic diversity

mRNA expression is not equivalent to rRNA abundance • Comparisons between rRNA and mRNA read distributions reveal subtle differences in read abundance • Artefact due to biases in rRNA sequencing? • Abundance != expression?

9. Visualizing results Once reads have been assigned to transcripts, we can explore their function through e.g. functional categories such as Gene Ontology or COG However these types of analyses tend not to be very informative as the categories are very broad Transcription Energy production and metabolism Lipid transport and metabolism

9. Visualizing results Beyond focusing on broad functional categories, we can also start to undertake systems based analyses Signalling networks Genes and proteins do not operate in isolation but form parts of interconnected functional modules Metabolic pathways By placing a bacterial transcripts within these functional contexts, we can understand how the microbiome functions at a molecular level Protein complexes

9. Visualizing results By mapping abundance of reads that map to E. coli homologs, protein interaction networks can be used as a template to provide more mechanistic insights – can we identify modules of genes that are coordinately expressed

9. Visualizing results How do we identify those systems which are differentially regulated across samples (e.g. systems up-regulated in disease vrs. health)? Development of a statistically robust platform (1) Gene set enrichment analysis (GSEA) requires discrete groups of functionally related genes (2) Network approach; identify groups of genes based on functional interactions AND their differential expression GSEA may define trpA-E as a single group of related genes, but group as a whole may not be significantly upregulated A network based approach may define a group based on their interactions AND expression

9. Visualizing results Metabolic pathways are among the most highly conserved and best characterized systems MG-RAST and MEGAN are automated metagenomic annotation tools that rely on KEGG A major problem with KEGG pathway definitions is that the boundaries of pathways are arbitrarily defined and links between pathways (i.e. functional relationships) can be lost

Nucleotide metabolism 9. Visualizing results Carbohydrate metabolism Metabolic network of a mouse gut microbiome Combining taxonomic and functional information allows us to identify contributions of individual taxa within a microbiome Lipid metabolism Amino acid metabolism The interconversion of pyruvate to acetyl CoA is largely mediated by proteobacteria Piecharts indicate relative contribution of each taxon to an enzymatic activity

9. Visualizing results In addition to representing protein interaction or metabolic relationships through network visualization, Cytoscape offers the ability to represent nodes (proteins) as pie charts or donuts, allowing the visualization of taxonomic contributions to discrete functions Genes involved in pathways associated with cell wall biogenesis Relative abundance and contributing taxa Relative abundance

10. Statistical considerations • There is no dedicated software or statistical tool for statistical comparisons of metatranscriptomic datasets • Number of biological replicates? (at least two, preferably at least three but these are expensive experiments) • Power analyses? • Differential expression of individual genes • Gene set enrichment analyses • Ultimately metatranscriptomics could be viewed as hypothesis generating requiring subsequent targeted validation

10. Statistical considerations • While there are no dedicated tools for metatranscriptomicsanalyses, tools used for RNA Seq offer potential • DESeq, EdgeR, ALDEx • Alternatively simply rely on fold change • Challenges include which genes to include (minimum RPKM?) • Differentially expressed genes can be subsequently analysed through Gene Set Enrichment Approaches based on (for example) hypergeometric distributions

We are on a Coffee Break & Networking Session

Canadian Bioinformatics Workshops