disentangling target genome and endosymbiont at read level

[PHSchi]

Hi!

Recently I got data at my hands were a single lane of GAIIx was sequenced from genomic DNA. >6GB output, all high qual. However, GC plot of reads shows two distinct peaks (larger at 37% -> target genome, smaller at around 50%). Seeing this and knowing the source of the DNA the second peak seems to come from an endosymbiont (or bacterial contamination). When I assemble with velvet (already tested cc=50 and large kmers) or Ray I get a genome of around 2MB (far to small) with bad cegma and also none of the stuff that should be in there, although blast hits for the right organism. Questions is: how to separate the endosymbiont from the target, possibly at read level?

Any help highly appreciated.

[maubp]

I'd try splitting the reads by GC, and assembling the high GC and low GC pools separately.

[PHSchi]

my first though, but

Hi!

That was my first thought as well, but then I loose all the reads of high GC from the target genome, i.e. they will be included in the other one - it is a curve of GC content and thus the target genome has to have regions with ~ 50% GC as well, whereas the possible endosymbiont has have reads of lower GC as well. I could take all the reads from a certain level of GC upwards and hope to assemble a bacterium out of those, then extract all the reads that went into the bacterial genome from the complete set of reads and hope to end up with more or less pure target genome, but is this sensible and feasible?!?

[Rockx]

I am also having the same problem.

What programs should I be using to separate reads based on GC content? A search of these forums only revealed replies such as "there are many programs that do this" but with no examples.

Any help would be appreciated, cheers!

[maubp]

What is your favourite scripting language? e.g. BioPerl/Biopython/etc would all make it easy to write a quick script to filter FASTQ on GC content. Also, do you have paired end data - and if so presumably you might want to filter at the pair level? That makes things a little more complicated...

[swbarnes2]

The other way to split is to align all the reads to your target organism stringently, and the velvet only the unmapped reads. Then, either figure out what your mystery contaminant is from the velvet, or include the velvet contigs in your genome alongside your desired organism, so that the reads will align to that, instead of beign forced somewhere in your target organism genome.

[Rockx]

Thanks maubp and swbarnes2. Indeed, I did end up editing the DynamicTrim perl script to include and option for GC trimming, this deals with paired end data fine.

swbarnes2, thanks for this tip. However, I 'm unable to do this as I am assembling de novo. Makes things a bit tougher.

[koadman]

You could try running an assembly pipeline designed explicitly to deal with mixes of organisms. metAMOS seems to be one such option:

https://github.com/treangen/metAMOS/wiki

It uses metagenome taxonomy analysis to figure out which organism each scaffold group comes from and creates a separate assembly fasta file for each organism. Looks like it's under very active development at the moment.

[polyatail]

There are a number of tools out there that attempt to cluster or classify reads or contigs by sequence-intrinsic properties (i.e. k-mers, protein domains). Check out TETRA, WebCarma, TACOA or PhyloPythia.

[koadman]

The authors of PhyloPythia have an interesting comparison of the nucleotide composition-based methods to a sequence identity/homology-based method (MEGAN) in the 50 pages of supplemental material for this 1.5 page paper:

http://www.nature.com/nmeth/journal/v8/n3/full/nmeth0311-191.html

I didn't notice whether they ran a nucleotide or amino acid blast search for MEGAN, but either way, it seems that using homology information gives pretty darn good results compared to the composition methods (among which PhyloPythiaS seems to be superior).