This is a test version of Biostars. For the public version, visit https://www.biostars.org.
samtools idxstats not removing ChrM

I am trying to remove ChrM from my ChIP-seq data. Below is my pipeline for one sample up to where I am having the issue (samtools idxstats). The output file from samtools idxstats is the same size as the input so it doesn't look like it is removing anything. I've used this before with success so I have no idea what I am doing wrong this time around. Any help is appreciated.

#Align to genome using bowtie2
bowtie2 -p 8 -q --local \
-x chipseq/reference_data/mm10/mm10 \
-U chipseq/raw_data/PPARg/PPARg_Input_1.fastq.gz \
-S chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_unsorted.sam

#Changing file format from SAM to BAM 
samtools view -h -S -b \
-o chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_unsorted.bam \
chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_unsorted.sam

#Sort BAM 
sambamba sort -t 8 \
-o chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam \
chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_unsorted.bam 

#Index BAM (need .bai for remove ChrM reads)
samtools index chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam

#Remove reads mapping to mitochondria
samtools idxstats chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam | cut -f 1 | grep -v chrM | xargs samtools view -b chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam > chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted_woChrM.bam
samtools chip-seq

can we see the output from this? samtools idxstats chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam

conversely just as an experiment try isolating the mitochondrial alignments: samtools idxstats chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam | cut -f 1 | grep chrM | xargs samtools view -b chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam > chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted_onlyChrM.bam

i'm not sure this is doing what you expect?

you should try filtering your bam directly and then running idx stats

samtools view -h your.bam | awk '{if($3 != "chrM"){print $0}}' | samtools view -Shb - > filter.bam

enter image description here

I think I got it to work using "MT" rather than "ChrM" (code below)

samtools idxstats chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam | cut -f 1 | grep -v MT | xargs samtools view -b chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted.bam > chipseq/results/bowtie2/PPARg/PPARg_Input_Rep1_aln_sorted_woMT.bam

enter image description here

1 answer

You mitochondrial genome is named MT not chrM...

Log in to answer this question.