Samtools Depth And Gaps

Question

[Juliofdiaz]

IS there a way that samtools depth include gaps? Right now every time I run samtools depth the data just skips the base positions without coverage. Thanks

Answer 1

[Pierre Lindenbaum]

you can try to fill the gaps with following awk command:

samtools depth file.bam |\
awk 'BEGIN { prev_chr="";prev_pos=0;} { if($1==prev_chr && prev_pos+1!=int($2)) {for(i=prev_pos+1;i<int($2);++i) {printf("%s\t%d\t0\n",$1,i);}} print; prev_chr=$1;prev_pos=int($2);}'

[Zev.Kronenberg]

Good God man you are an awk wizard ~!

[Lee Katz]

Any quick way to also include the beginning/ends of the reference sequence? For example in mine, bases 1-6 have zero coverage and the output starts when it is at base 7 of the chromosome. When it gets to a depth of zero in the middle of the chromosome, it reports it thanks to the awk command.

[Pierre Lindenbaum]

use bedtools to get that information.

[Lee Katz]

I got something bioperl-ish that seems to work

samtools depth $sorted |\
  perl -MBio::Perl -e 'while(<>){chomp; @F=split /\t/;$depth{$F[0]}{$F[1]}=$F[2];} $in=Bio::SeqIO->new(-file=>"'$ref'");while($seq=$in->next_seq){$max{$seq->id}=$seq->length;} for $chr(keys(%max)){for $i(1..$max{$chr}){$depth{chr}{$i}||=0; print join("\t",$chr,$i,$depth{$chr}{$i})."\n";}}'\
  > "$sorted.depth"

where $ref is a reference genome in fasta format defined in the shell environment, $sorted is a sorted bam