Coverage from Cram input
Created by: Schmurf87
Hi, I'm trying to retrieve average coverage over chr20 from a cram file.
I used the following perl script:
use strict;
use warnings;
use Bio::DB::HTS;
my $hts= Bio::DB::HTS->new(-bam => "aln.cram",
-fasta =>"ref.fa",
);
my ($coverage) = $sam->get_features_by_location(-type=>'coverage',-seq_id=>"chr20");
my @data = $coverage->coverage;
my $total;
for (@data) { $total += $_ }
my $average_coverage = $total/@data;
With it and my cram file, I got an average coverage of 0.15, which seems particularly low. I wanted to check if I got the same average coverage with a bam file. To do so, I converted my cram file to bam using the following command:
samtools view -b -h aln.cram -o aln.bam -T ref.fa
I ran the same perl code for this bam and get an average coverage of 5.16 ...
Why is there a difference and is there any way to get the same average coverage result?
Thank you in advance for your answer.