GitLab

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.