Hi all,
I want to remove all mitochondrial reads from my BAM files, but after doing so the number of mapped reads increased. Here's what I did after indexing:
samtools idxstats .bam | cut -f 1 | grep -v chrM | xargs samtools view -b .bam > clean.bam
However comparing with wc -l says:
samtools view -F 0x904 -c clean.bam
19452291
samtools view -F 0x904 -c .bam
18328474
The second number matches the summary output by tophat2 (fr-frstrand)
any suggestions?
2 answers
Your approach looks pretty strange to me, and I would instead do the following:
samtools view -h yourbam.bam | grep -v 'chrM' | samtools view -b > yourbamwithoutchrM.bam
Hm, you are probably right. I'm not sure why OP would want to remove reads mapping to chrM in the first place. But I think the event you describe is rare and discordant reads like that aren't really desirable.
That said, I'm not sure if there is a good reason to start tampering with bam files.
I found the issue! this is with samtools 1.2 with reads mapped to gencode v25
the increased number is due to chr1 doubling in number after sorting by position (default). I still dont know why this happens though.
no chrM:
chr1 248956422 5140060 0
all:
chr1 248956422 2570030 0
maybe it has something to do with the reads named, GL00000x.x...
Log in to answer this question.
-
Seeking advice on methods to calculate ENCODE ChIP-seq data standard M_DISTINCT
written by kalavattam •
-
samtools idxstats not removing ChrM
written by Jen •
-
paired-end reads, idxstats, and samtools and grep
written by Bogdan
-
Differing number of reads for read 1 and 2 in fastq's from a subset bam file
written by michael •
-
How to plot alignment statistics; how many reads are mapped to the genome?
written by anamaria
-
Aligned read loss after SAM to BAM conversion?
written by saddamhusain77 •
-
Samtools idxstats fails to load index
written by lvande •
-
Is there a way to run samtools sort, index, and view on multiple files at once?
written by ely.j.fish •
-
Different number of total reads in different alignments using the same original readset
written by malteherold •
-
filter mapped and unmapped bam file
written by mgadrianam •
I am using the exact same command, except that I put chrM in quotation marks. Also, don't forget the -h flag in the last
view, as many tools need the BAM header to work properly (or work at all).He's outputting to BAM (-b) so it will always have a header. It's only SAM where it's optional for some reason.