A workflow for filtering

[Ackia]

Hi!

First of all thanks to the community who often solve my problems before I even have time to ask the question.

I have a problem though. I’m working with datasets from 454 sequencing and I know one organism to be present and want to filter out those reads to find possible other organisms by analyzing the remaining reads.

From what I have managed to grasp running around in these forums and other forums it should be an easy task if I could extract the data from an alignment. MAQ as well as Bowtie have these functions built into them (to dump the unaligned sequences) but since I’m using 454 they won’t work.

So my question is, is this doable with BWA and samtools? I think it should be doable but I’m left with mostly incomplete clues in the forums and I’m a bit lost right now.
I know that samtools view with the flag –f should be able to dump the reads but I’m unsure of the flag and even more unsure on how to handle the output (ideally I will have an output that I can incorporate into MIRA for assembly (or velvet / Ray)).

[nilshomer]

samtools view -f 4 <in.bam> > unmapped.sam

You can then use Picard's SamToFastq.jar utility.

See http://picard.sourceforge.net/explain-flags.html, http://picard.sourceforge.net
and http://samtools.sourceforge.net

[Ackia]

samtools view -f 4 <in.bam> > unmapped.sam

You can then use Picard's SamToFastq.jar utility.

See http://picard.sourceforge.net/explain-flags.html, http://picard.sourceforge.net
and http://samtools.sourceforge.net

Great! Then I was not totally out on thin ice at least.

Is there another way to convert between sam files and FastQ? Also, since Mira for example will whine if I do not provide it with 3 files for assembly (.fa .qual and .xml) I think it might be simpler to run Mira first (the datasets are easily handled with the current cluster anyway so computing power is not a problem) then map it using BWA, convert and BLAST the unmapped results. Should filter out the problematic reads/contigs from the dataset

What I'm trying to emulate is close to the PathSeq workflow but using different tools (since they use the illumina platform and MAQ but my datasets are 454).

Best Regards

[pcs_murali]

Hi Ackia,

You can samtools JDK toolkit and extract the information.

Chandra