Separating reads from different organisms

[christinehicks]

My RNAseq data contains reads from two different organisms (a host and a pathogen). I wish to separate these reads so I can work on both species separately. I plan to do a de novo assembly on the host reads, which requires the reads to be in fastq format, not bam format.

In order to separate the reads, I have aligned to the pathogen genome using tophat2, which produced an unmapped.bam file containing the host reads. I used samtools sort to sort the bam file, then ran BEDtools bamtofastq to recreate the fastq files containing only the host reads. However, when I run bed tools I always get an empty output file. Looking at the error messages, this is because there are unpaired reads in the bam file. I ran tophat2 with the 'no-mix' parameter, which should prevent the reads from being separated but I still get unpaired reads in my unmapped.bam file.

I'm pretty stumped, does anyone have any experience in using bedtools bamtofastq? Or maybe a better suggestion for separating reads from two different organisms?

Thank you!

[GenoMax]

Give BBMap a try as an alternative: http://seqanswers.com/forums/showthread.php?t=41818

[Brian Bushnell]

BBMap allows you to output in fastq, and it has various output streams including "outu" for only unmapped reads and "outm" for only mapped reads. You can also use the "po" flag (meaning "paired only") to specify that unpaired reads should be considered unmapped - specifically, "po=f" (the default) will send mapped but unpaired reads to "outm" while "po=t" will consider unpaired reads as unmapped and thus send them to "outu".

I slightly messed up the commands in the thread mentioned by Genomax, and have corrected them. This is the corrected version:

(index)
bbmap.sh ref=reference.fasta -Xmx29g

(map)
bbmap.sh in=reads.fq outm=mapped.fq outu=unmapped.fq -Xmx29g

That will split it into two fastq files. If you change the output file extensions to ".sam" instead of ".fq", you will get sam output. If you have paired reads in 2 files, you would do this:

bbmap.sh in1=reads1.fq in2=reads2.fq outm1=mapped.fq outm2=mapped2.fq outu1=unmapped1.fq outu2=unmapped2.fq -Xmx29g

The "-Xmx29g" indicates how much memory java is allowed to use; set that at ~85% of your physical memory.

If you have assemblies for both the organism and pathogen, you can use BBSplit instead of BBMap; in that case, just reply here if you need help. BBSplit works fine with RNA-seq data but its defaults are set for DNA so they would need adjustment.

[khughitt]

Hi Christine,

By "no-mix", do you mean "--no-mixed"? It might be helpful to post the entire command you are using.

Also, Tophat ships with a conversion tool, bam2fastx, which you might want to try for converting unmapped.bam back to a FASTQ file.

Finally, not necessarily related to your question, but are you doing any QC/trimming before attempting to map? This is generally useful to save pathogen reads from going to unmapped.bam when they have low quality bases on the ends, etc.

HTH,
Keith

[Michael.Ante]

You may try to combine both genomes in one reference file and align all reads against it. This maybe bypasses the problem of homologous regions in both species.
Afterwards, you can separate the bam-file species-wise and create the fastq-file as Keith already mentioned.