How does BWA deal with reads that are shorter than expected?

[Jlap]

Hello all,

I am very new to bioinformatics so this might just be a stupid question....

I am using BWA and Samtools to map my sequence reads on a reference mitogenome, but the number of reads that can be mapped is very low - only about 100 reads. I suspect that BWA is only mapping the "complete" reads, so the reads that are 100bp (we did PE reads of 100bp on an Illumina HiSeq).

However, since we are working with highly degraded samples, we expect most of the fragments to be not longer than ~70bp. So if shorter reads are automatically kicked out, I am loosing >90% of my reads.

Does anyone know if it is the BWA script that kicks "shorter-than-expected" reads? Or how I can get to the shorter reads? Any help, suggestions, thoughts, comments, ect, are highly appreciated!

Cheers,
Johanna

[Richard Finney]

Can you show us your bwa commands ?
Did you index the "mitogenome" using the same BWA version that you aligned with? Mixing version <6 with >= 6 is wrong.

How many input reads?

Do the input reads look good ? Try hand BLATing a few hundred random reads using command line blat against your custom mitogenome? Not getting any really good hits ?

Try these commands to do a little QA ...

wc input.fastq # divide by 4 for number of reads

grep "NNNNNNN" input.fq | wc
# check for "NNNNN" (bad) reads; you might get a few near start/end

[arvid]

...
(we did PE reads of 100bp on an Illumina HiSeq).

However, since we are working with highly degraded samples, we expect most of the fragments to be not longer than ~70bp.
...

Did you clip adaptor sequences from your reads? If you mean that you have a lot of fragments around 70 bp, and you did paired end 100 bp sequencing, you would have read through your fragments into the adaptors; and these sequences might prevent your reads from aligning to the reference. If you didn't already, have a look at trying a paired end clipper/trimmer such as Trimmomatic.

[Jlap]

Thanks for the posts -

In answer to the first question, I don't think I'm using different versions of BWA - I was only introduced to it a couple of days ago.. Basically, I'm following this manual for the BWA commands: http://sourceforge.net/apps/mediawik...stall_software. (BWA/Samtools for dummies). So that's indexing the reference, sampe the reads, converting to .sam files and then to .bam to view with Tablet.

Thanks for the QA commands; the number of reads comes up around 370,000 per file (8 files per read). It doesn't find any NNNN reads - I can only assume those were already taken out before I got my hands on the data.

According to Arvid's suggestion, I've been trying to get Trimmomatic to work, but it keeps dropping 100% of the reads. I have the feeling this is because I am doing something wrong in creating the adapter fasta.. I've tried half a dozen different set-ups of the fasta (forward, reverse, merged, not merged, with /1 - /2, without, ect). I'm not sure how many other creative ideas I can apply here.. Does anybody have some pointers here?

[Jlap]

Just in case anyone is interested: I came across an adapter-trimming tool which performs really well to trim read-through adapters, and it's very easy to use. It's called Cutadapt and can be found here:

http://journal.embnet.org/index.php/embnetjournal/article/view/200

http://code.google.com/p/cutadapt/