samtools idxstats not removing ChrM
1
1
Entering edit mode
13 months ago
Jen ▴ 70

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 • 1.1k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

looks good :)

ADD REPLY
1
Entering edit mode
13 months ago

You mitochondrial genome is named MT not chrM...

ADD COMMENT

Login before adding your answer.

Traffic: 2793 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6