Hello,
I have a sorted bam file (bwa mem alignment). Now I need to mark and remove duplicates using samtools. Before it asks me to do collate and fixmate then sort and then markdup. I was wondering if bam output (without sorting) is already name sorted (bwa version 0.7.17)? So I could skip unnecessary sorting, collating steps and after alignment pass directly to fixmate?
1 answer
The alignment output (if without further manipulation) is grouped by name since fastq is grouped by name. Therefore you can pipe this directly into fixmate, e.g.:
bwa mem (...options) | samtools fixmate -m - - | samtools sort -o sorted.bam
The -m option of fixmate adds the mate tags which are required for the markdup command.
The file must also be coordinate-sorted whereas fixmate requires files to be grouped or sorted by name.
samtools markdup sorted.bam markdup.bam
So if I dont need any of the files just the last without duplicates, would this work in snakemake?
"""/Tools/bwa-0.7.17/bwa mem -t {threads} {ref} {input.fastq1} {input.fastq2} | /Tools/samtools-1.10/samtools fixmate -m - - | /Tools/samtools-1.10/samtools sort - | /Tools/samtools-1.10/samtools markdup -r {output.bam}"""
PS. thanks ATpoint, you are always replying to my posts.
Alignment is sorted by name since fastq is sorted by name
Actually, this assumption has gotten me into trouble on the odd occasion. For example BamUtils will complain that the raw output from an aligner is NOT sorted by name, if read-pairs are together, e.g. if a read-pair has 2 alignments, bamUtils wants to see them in the order R1,R1,R2,R2 (not R1,R2,R1,R2 like straight from the aligner). Unfortunately you still need to sort the raw alignment using samtools sort -n it to get it working. Others have speculated here that in theory, read-order can get slightly jumbled if the aligner is multithreading due to race conditions, depending on how it implements multithreading (not sure if this applies to BWA tho).
Minor correction: actually BamUtils clipOverlap --readName does work with R1,R2,R1,R2 sorting, as long as alignments are strictly alphanumerically sorted by the read-pair name (which the raw aligner output is generally not).
I am actually running fixmate and then sort by coordinates right straight after alignment, and it is working fine :) here in the manual (http://www.htslib.org/doc/samtools-markdup.html#AUTHOR) they say: "Typically the fixmate step would be applied immediately after sequence alignment and the markdup step after sorting by chromosome and position. Thus no additional sort steps are normally needed."
I see your point. I edited my answer now saying files are grouped by name. You are right that infact sorted by name would also mean that there is indeed a strict sort order, I was rather referring to the fact that reads will always be R1/R2 after paired-end alignment and this is all fixmate needs.
Log in to answer this question.
-
Why did the QV value decrease significantly after using YAHS for HiC scaffolding?
written by xinguok794 •
-
BWA alt-aware pipeline troubleshooting
written by Duy •
-
Samtools pipe in job script to modify BAM files
written by Cory •
-
Input and output of picard markduplicate for deepvariant
written by hello_bioinf •
-
paired reads have different names BWA MEM after samtools bam > fastq
written by robert.murphy •
-
SAMtools markdup and insert size in Qualimap
written by Egelbets •
-
samtools duplicate reads removing
written by Arslan Tariq •
-
Samtools Parse error at line 1
written by storm1907 •
-
remove duplicates with markdup
written by dimitrischat
-
Calling SNPs using bcftools mpileup
written by jaafari.omid
BWA does not provide sorted bam file as output. If you want to check whether the aligned file is sorted or not use the following command.