BWA and SAMtools

[seq_lover]

Hello Everyone,

I used BWA to align SOLiD mate pair reads (60,60) with parameters -n 8(total mismatch) -l 25 (seed) and -k 2 (mismatch in seed). I am getting a good mapping rate of around 65%. I will be using BFAST and Lifescope in future.

BWA outputs all the reads disregard of whether they were mapped, unmapped, mapped in pairs and other bitwise flags. To solve this problem I converted my SAM file to BAM file. As I am not interested in inversions or some unusual variants right now, I filtered out the SAM file so that it can be used for high confidence SNP and Indel calling. So I used:

samtools view -b -f 67 -f 31 -f 179 -f 115 old.bam > new.bam

67 and 31 (paired, mapped and properly paired) 179 and 115 (paired, mapped, properly mapped and both mapped reverse complimentary same strand)

Once I got the new.bam BAM, I sorted it and removed the duplicates usign samtools and then used mpileup to call for the SNPs and indels.

Below are my Yes or No questions:

1) This is my first time doing a NGS analysis. Am I doing things correctly? Is the order of steps I am performing correct?

2) As I only want to use high confidant reads I have filtered out all the unmapped, not properly paired reads. Do you think the flagwise bits I have used are correct.

Though I tried to remove the duplicates using samtools for my mate pair bam data but I can still see lot of mate-pair reads mapped to the same position as other mate-pair reads. Some people have suggested using Picard.

I think the problem could be because I used "trim 3'end" option in BWA. The reads that were duplicates before may not remain duplicates afterwards because the length of some reads got changed after trimming. Can anyone tell me how to resolve this issue.

PS: In future I am planning to try GATK but for now i think samtools variant calling can give me some idea about quality of mapping and my solid data.

Thanks.

Thanks

[swbarnes2]

It's not clear to me that bwa will treat two reads going in the same direction as properly paired. But I've never tried it, I could be wrong.

If all you want are the properly paired reads, do -f 2.

I'm not sure why you have to get rid of the unmapped reads anyway. Mpileup isn't going to use them, but someday, you may want to be able to look at them, so why not leave them in the file?

[seq_lover]

Thanks for the reply. This is my first analysis so i am going extra careful.

[sdriscoll]

It's not clear to me that bwa will treat two reads going in the same direction as properly paired. But I've never tried it, I could be wrong.

If all you want are the properly paired reads, do -f 2.

I'm not sure why you have to get rid of the unmapped reads anyway. Mpileup isn't going to use them, but someday, you may want to be able to look at them, so why not leave them in the file?

I like to remove unmapped reads from my BAM files as well. This is just to save disk space. When my lab does a sequencing run we get 500 GB of FASTQ data across 48 samples which balloons up to a lot more HD space once everything is aligned. If we ever care about unaligned reads I can re-run the alignments and parse them out.

[seq_lover]

I like to remove unmapped reads from my BAM files as well. This is just to save disk space. When my lab does a sequencing run we get 500 GB of FASTQ data across 48 samples which balloons up to a lot more HD space once everything is aligned. If we ever care about unaligned reads I can re-run the alignments and parse them out.

Hi,

If you use BWA can you comment on my questions. I mean "Am i using the correct flags" because when I used these flags my BAM file reduced from 18 gb to 4 gb. I dont care about the unmapped or single pair mapped reads.

[sdriscoll]

Can you explain why you don't want different reads to align in the same place in the genome? Normally people might exclude reads that align to multiple genome locations (repetitive sequences) but it's a normal and desirable part of RNA-Seq for multiple reads to align to the same locations within the genome. That's how we are able to measure gene expression and part of how we can consider SNPs to be significant.