Proportion of mitochondrial reads atac seq
1
0
Entering edit mode
12 days 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 • 274 views
ADD COMMENT
1
Entering edit mode
12 days ago
samtools idxstats yourBAMfile.bam | cut -f 1,3
ADD COMMENT
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?

ADD REPLY
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")'%'
ADD REPLY
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?

ADD REPLY
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.

ADD REPLY

Login before adding your answer.

Traffic: 1138 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