To calculate the length of aligned genome

Question

[sambioinfo2018]

Dear All,

I have a bam file aligned to a reference genome. I would like to know if there is a way to calculate the over all length of aligned genome.

Thanks in advance. Sam.

Answer 1

[Bastien Hervé]

genomecov should do what you are asking for

Edit :

May take a while :

output=$(bedtools genomecov -ibam in_file.bam -d -split) && (echo "scale=2; $(echo "$output" | cut -f3 | grep -vc "0")/$(echo "$output" | wc -l)" | bc)

[sambioinfo2018]

Thanks Bastien. Is this the command to do so: bedtools genomecov -ibam bamfile.bam -bga -g genomelength.txt

[Bastien Hervé]

More something like this (get the depth at each position):

bedtools genomecov -ibam input_file.bam -d -split > output_depth.txt

From the output file (output_depth.txt), count the number of line with NOT 0 in the 3rd column (with awk for example) divided by the number of line in the file. You will get the percentage of reference positions covered by at least one read.