Proportion of mitochondrial reads atac seq
1
0
Entering edit mode
10 weeks ago
conarial ▴ 10

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?

samtools seq atac • 377 views
1
Entering edit mode
10 weeks ago
samtools idxstats yourBAMfile.bam | cut -f 1,3
1
Entering edit mode

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?

0
Entering edit mode

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")'%'
0
Entering edit mode

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?

0
Entering edit mode

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.