Algorithms for Multisample Read Binning Student Gabriel Ilie

Algorithms for Multisample Read Binning Student: Gabriel Ilie Major advisor: Ion Măndoiu Associate Advisor: Sanguthevar Rajasekaran Associate Advisor: Yufeng Wu University of Connecticut November 2013

Outline • • • Motivation Previous approaches Algorithms Results Ongoing work

Human microbiome • • The microorganisms inhabiting our bodies comprise the microbiome. Microbes outnumber our own cells by 10 to 1 [Wooley et al, 2010]. Diabetes, obesity, cancer and even attractiveness to mosquitoes seem to correlate with changes in our microbiome [Turnbaugh et al, 2010]. To better understand our own condition we also need to understand the composition of the microbial communities inhabiting our bodies and how they interact not just between themselves but also with their habitat (e. g. the host organism).

Single organism studies • • • they rely on clonal cultures most microorganism cannot be cultivated suffer from amplification bias microbes do not live in single species communities members of these communities interact not just with each other but also with their habitats, which includes the host organism Data obtained from clonal cultures is highly biased and does not capture a true picture of microbial life.

Transcriptomics • • Genomes only provide information about the potential function of organisms. Having a gene does not mean that this gene is also expressed inside the host or during a particular condition. The transcriptome is the set of all RNA molecules produced in one or a population of cells. In order to understand the physiology of the microorganisms we need to know their transcriptome.

Metatranscriptomics • • • The transcriptome of a community is the union over the transcriptomes of each of its members. In metatranscriptomic studies: o bulk RNA is extracted directly from environmental samples o the RNA is reverse-transcribed o the resulting RNA-Seq libraries are sequenced Metatranscriptomic studies are essential in order to understand the physiology of the microbiome.

Challenges of working with metatranscriptomic data 1. The volume of sequencing data is several orders of magnitude larger than single organisms. 2. Reads can come from hundreds of different species, each with a different abundance level. 3. In addition to having a range in the abundance levels of the microorganisms, genes are expressed at drastically different levels. 4. Genes usually have multiple isoforms. Metatranscriptomics has to deal with all of the challenges of metagenomics (1 and 2) plus some extra challenges (3 and 4), therefore algorithms devised for metagenomic data can also be applied to metatranscriptomic data.

(Meta)Transcriptome assembly types Genome independent reconstruction (de novo): • de Bruijn k-mer graph ❖ Velvet (2008) ❖ Trinity (2011) Genome guided reconstruction: • • • spliced read mapping exon identification splice graph ❖ Cufflinks (2010) ❖ TRIP (2012)

Outline • • • Motivation Previous approaches Algorithms Results Ongoing work

Clustering reads into bins • • • Analysis of environmental samples is difficult. To simplify the assembly process, many metagenomic tools have been developed to cluster the reads into bins (i. e. species). Algorithms developed for binning metagenomic reads can also be applied to (meta)transcriptomic reads (bins represent transcripts instead of species).

Types of reads binning algorithms Genome dependent ❖ Compost. Bin(2008) ❖ Metacluster(2012) DNA composition patterns. G+C content, dinucleotide frequencies vary amongst species. Drawbacks: o achieve reasonable performance only for long reads (800~1000 bp, [Wu et al, 2011]) o NGS technologies produce short reads • • • Genome independent ❖ Abundance. Bin(2011) ❖ Multi. Bin(2011) K-mer frequencies are usually linearly proportional to a genome’s abundance. Sufficiently long k-mers are usually unique. Works with short sequencing reads. Drawbacks o group together reads from different species if they have close abundance levels o do not perform well on species with low abundances • •

Metacluster (2012) ❖ Two round unsupervised binning algorithm: ➢ the first round clusters high-abundance reads ■ filter reads with 16 -mer that appear less than T times ■ the reads are grouped based on shared 36 -mers ■ the groups are clustered using 5 -mer distributions ➢ the second round clusters the remaining reads (lowabundance) ■ those with unique 16 -mers are discarded ■ the reads are grouped based on shared 22 -mers ■ the groups are clustered using 4 -mer distributions ❖ Advantages: ➢ better binning of reads from low abundance species ➢ uses both techniques: k-mer abundances and DNA composition patterns

Multi. Bin (2011) • • • Processes multiple samples (N > 1) of the same microbial community. Clusters the reads into b bins (b is the number of species). Binning algorithm: o all reads are pooled together o a graph G=(V, E) is generated § V is the set of reads § edges connect reads with substantial overlap (50 bp) o greedily partition the vertices into a set, the tags, s. t. each read is either a tag or affiliated with one which substantially overlaps it o cluster the set of tags using Vt=(ct 1, ct 2, …, ct. N), cti=number of reads from sample i which substantially overlap tag t o each non-tag read is assigned to the same bin as its affiliated tag ❖ b needs to be known in advance or estimated somehow. ❖ Uses abundance differences from any of the samples to tell low abundance species apart.

Outline • • • Motivation Previous approaches Algorithms Results Ongoing work

Our approach We propose a novel method for unsupervised abundance-based multiple samples reads binning algorithm

Pseudo-exons are defined as substrings whose k-mers and (k+1)mers appear in the same transcripts with the same multiplicity. T 1 A B C T 2 A B C T 3 A B C T 4 F E Exon signatures: A: T 1 x 1 T 2 x 2 T 3 x 1 B: T 1 x 1 T 2 x 1 T 3 x 1 C: T 1 x 1 T 2 x 1 T 3 x 1 D: T 1 x 1 T 3 x 1 E: T 3 x 1 T 4 x 1 F: T 3 x 1 T 4 x 1 D A E F D Pseudo-exons: A BC D E F Notice that exons E and F form different pseudo-exons because in T 4 they are in reverse order compared to T 3.

K-mer counting • • we count k-mers and (k+1)-mers using Jellyfish (2011) Jellyfish was designed for shared memory parallel computers with more than one core, it uses several lock-free data structures and multi-threading to count k-mers much faster than other tools formally, for a given value of k, we count the number of occurrences of all k-mers in each of the N samples our algorithms assume we have strand nonspecific data, therefore the counts of complementary k-mers are summed together as they are indistinguishable we store k-mers in canonical form (the smaller value lexicographically between a k-mer and its reverse complement) we combine the counts over all the samples into a list of Ndimensional vectors the maximum value for k supported by Jellyfish is 31

De Bruijn graph • We construct the de Bruijn graph; vertices are k-mers and edges are (k+1)-mers. • • • Because we store vertices in canonical form, each vertex represents two k-mers: itself and its reverse-complement. We add an edge between any two k-mers if and only if there is a (k+1)-mer in the reads such that its prefix of length k in canonical form matches one of the vertices, while its suffix of length k in canonical form matches the other one. We will sometimes use vertices to refer to k-mers and edges to refer to (k+1)-mers.

De Bruijn graph ACG CGT CGC GCG CGA TCG ATC GAT Vertices ACG CGC CGA ATC Edges ACGC CGCG ACGA ATCG ● The lists of k-mers and (k+1)-mers define an implicit representation of the de Bruijn graph, therefore we don’t need to construct it explicitly. ● Relative to the canonical form of a vertex, for each edge we can say whether it is incoming or outgoing: ○ if the k-mer matches the 5’ end of either forms of an edge then that is an outgoing edge ○ else it is an incoming edge

Error removal • • A common approach is to remove k-mers which have counts lower than t >= 1. We found that even for t = 1, we lose too much information, because removing unique k-mers compromises the results for ultralow abundance transcripts. We found that “tip removal” and “bubble removal” give much better results [Zerbino et al, Velvet, 2008]. These methods use the structure of the de Bruijn graph instead of coverage information to remove k-mers affected by sequencing errors.

Tip and bubble error removal • • • When a read contains a sequencing error o the first few k-mers may be correct, until they start to overlap the position where the error occurred o this creates a branch going out of the “correct” path o this new branch will either end in a leaf (creating a “tip”), or if the read is long enough, the k-mers will stop overlapping the error and the branch will merge back into the path (creating a “bubble”) A “tip” is a chain of nodes that is disconnected on one end o we expect the majority of tips to have a maximum length of 2 k o removing tips is straightforward; we remove all tips which have a length up to some threshold o removing a tip does not disrupt the connectivity of the graph Implementing “bubble” removal is still an ongoing work.

Partitioning the de Bruijn graph • • • From the de Bruijn graph we want to extract putative pseudo-exons. These putative pseudo-exons, if we ignore self-edges, correspond to paths or chordless cycles in the de Bruijn graph. We use the structure of the graph and the vectors with the counts to do the partitioning.

Partitioning the de Bruijn graph

Partitioning the de Bruijn graph We have the following cases: If a vertex has indegree and outdegree equal to at most 1 (it is on a path), we do nothing. If a vertex has outdegree (indegree) greater than 1, o then we remove all outgoing (incoming) edges from that vertex o however, we keep the most abundant edge if the ratio between its abundance and the sum of the abundances of the other edges is higher than a threshold 0 < e < 1 The value of e should be close to 1 (e. g. 0. 97) • • •

Our approach We propose a novel method for unsupervised abundance-based multiple samples reads binning algorithm

Outline • • • Motivation Previous approaches Algorithms Results Ongoing work

Test data - error free • • GNF Atlas [Su et al, 2004] is a dataset which contains information about the expression levels of a set of genes in several human tissues From this dataset we extracted the expression levels of 19, 371 genes in 10 human tissues We used only one isoform per gene We simulated 30 million error free RNA-Seq pairedreads of length 50 from this dataset using a tool called Grinder (2012)

Test data - with sequencing errors • • Grinder was very useful for simulating the error free data, however when we wanted to introduce errors its long running time became an issue. Instead, we simulated sequencing errors by using the error free data. We simulated only one type of errors, substitutions, because these are the most common type found in Illumina datasets. We introduced substitutions into the error free reads with a probability of 0. 1%, 0. 5% and 1% per base.

K-mer counts in the simulated data transcripts #30 -mers #31 -mers reads error 0% 0. 1% 49, 572, 543 49, 611, 691 0. 5% 1% k #correct k-mers #incorrect k-mers #missing kmers percentage of incorect k-mers 30 49, 512, 839 0 59, 704 0% 31 49, 546, 279 0 65, 412 0% 30 49, 508, 364 109, 380, 690 64, 179 68. 84% 31 49, 540, 750 108, 528, 925 70, 941 68, 66% 30 49, 483, 820 404, 415, 622 88, 723 89. 1% 31 49, 510, 299 402, 714, 823 101, 392 89. 05% 30 49, 433, 000 708, 460, 569 139, 543 93. 48% 31 49, 446, 451 706, 531, 643 165, 240 93. 46%

Efficiency of different error removal techniques Error removal method error #correct 31 mers #incorrect 31 mers #missing 31 mers none 0% 49, 546, 279 0 65, 412 non-unique in at least 1 sample 0% 46, 520, 486 0 3, 091, 205 non-unique in at least 1 sample 0. 1% 46, 257, 210 6, 906, 069 3, 354, 481 remove tips <= 21 over the union of the samples 0. 1% 49, 502, 470 21, 772, 808 109, 221 remove tips <= 60 over the union of the samples 0. 1% 49, 236, 255 12, 120, 069 375, 436 remove tips <= 21 sample-by-sample and over the union of the samples 0. 1% 49, 127, 456 10, 998, 956 484, 235 remove tips <= 60 sample-by-sample and over the union of the samples 0. 1% 48, 707, 567 5, 304, 983 904, 124

Results for the graph partitioning error removal technique edge removal threshold #pseudo-exons >= 50 bp and do not have wrong 31 -mers #30 -mers %transcriptome covered by col 5 0% none 1 60, 285 49, 299, 665 99. 5% 0% none 1 65, 051 49, 239, 335 99. 3% 0. 1% remove tips <= 60 over the union of the samples 0. 97 416, 853 60, 335 37, 815, 256 76. 3% 0. 1% remove tips <= 60 sampleby-sample and over the union of the samples 0. 97 228, 841 77, 010 46, 699, 963 94. 2% The first row (green) uses all k-mers, the counts are computed from the transcripts.

Outline • • • Motivation Previous approaches Algorithms Results Ongoing work

Ongoing work • • • We believe that bubble removal will help us get rid of most of the erroneous k-mers which are still present in the data after tip removal. We want to incorporate an error correction algorithm, SEECER (2013), to correct the erroneous reads before doing starting the k-mer counting. Currently our algorithms do not take into account strand-specific RNA-Seq data. Optimizing the algorithms to take advantage of this information, when available, represents another opportunity to improve the results of this approach.

Q&A