Proportion of mitochondrial reads atac seq

Question

[conarial]

I have atac seq data and would like to calculate (or view) the proportion of mitochondrial DNA in our sample. So how much mitochondrial DNA reads vs the whole number of reads. As far as I could see, samtools flagstats does not help me there. Could anyone help me out?

Answer 1

[benformatics]

samtools idxstats yourBAMfile.bam | cut -f 1,3

[conarial]

okay, so this gives me the chromosomes and the number of mitochondrial reads, right? So I would just have to add them add and divide by the total number of reads?

[ATpoint]

What @benformatics says, just automated, assuming name is chrM:

mtReads=$(samtools idxstats $1 | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats $1 | awk '{SUM += $3} END {print SUM}')
echo 'chrM reads:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'

[conarial]

Another beginner question, but how do I find out how the mitochondrial dna is connotated, because with "samtools idxstats <my_bam_file>" and look at column1, I only get the scaffold points from "sca1" to "sca787", so I cannot grep "chrM" ... can you help me there?

[benformatics]

You need to look at your reference FASTA (or the resource from where it was derived) and find out which one is the mitochondria. If you have GTF you could try to see which scaffold the mitochondrial genes are present in.