This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Samtools Depth And Gaps

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

samtools depth-of-coverage
• 6,676 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Juliofdiaz Juliofdiaz

1 answer

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);}'
• 2 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Good God man you are an awk wizard ~!

• 0 views • link
written by Zev.Kronenberg Zev.Kronenberg

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.

• 0 views • link
written by Lee Katz 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

• 0 views • link
written by Lee Katz Lee Katz

Log in to answer this question.

More posts like this