9
$\begingroup$

I am trying to use samtools depth (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage at every position:

cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1
3088286401

I would like to know how to run samtools depth so that it produces 3,088,286,401 entries when run against a GRCh38 bam file:

samtools depth -b $bedfile -a $inputfile

I tried it for a few bam files that were aligned the same way, and I get differing number of entries:

3087003274
3087005666
3087007158
3087009435
3087009439
3087009621
3087009818
3087010065
3087010408
3087010477
3087010481
3087012115
3087013147
3087013186
3087013500
3087149616

Is there a special flag in samtools depth so that it reports all entries from the bed file?

If samtools depth is not the best tool for this, what would be the equivalent with sambamba depth base?

sambamba depth base --min-coverage=0 --regions $bedfile $inputfile

Any other options?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ Out of interest, does -aa at least give the same number of bases each time? -a -a (or -aa) output absolutely all positions, including unused ref. sequences There is also an old but relevant discussion about -a and BED behaviour on github/samtools/#374 $\endgroup$
    – Sam Nicholls
    Jun 7, 2017 at 9:59
  • $\begingroup$ I tried with -a and with -aa, and given a bed file, it seems to produce the same number of entries. $\endgroup$
    – 719016
    Jun 7, 2017 at 12:35
  • 2
    $\begingroup$ Two samtools commits from May 2016; Added 30ish new tests for depth/mpileup -a and -aa and Depth/Mpileup -a fixes. Fixes #374 are somewhat enlightening. My gut feeling is that although possible, using -a/-aa with a BED file is not entirely supported with samtools depth (i.e. has potentially undefined behaviour) and perhaps you would want to investigate another tool as you've suggested? $\endgroup$
    – Sam Nicholls
    Jun 7, 2017 at 13:43
  • 1
    $\begingroup$ I dont work with human genomes, so excuse my naive question...why do you need to provide a bed file if you want every position? $\endgroup$
    – Nathan S. Watson-Haigh
    Jun 7, 2017 at 14:07
  • $\begingroup$ If every reference defined in the BAM file is touched by at least one read, then -a and -aa will produce the same output. The second option, -aa is useful if you want to produce zero-depth information for sequences that are completely uncovered. $\endgroup$
    – gringer ♦
    Jun 7, 2017 at 22:24

2 Answers 2

Reset to default
6
$\begingroup$

You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.

bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov"

You can also provide a BED file if you just want to calculate in the target region.

Improve this answer
$\endgroup$
2
$\begingroup$

Just replace the -a option with -aa:

samtools depth -b $bedfile -aa $inputfile

I see that you're using the GRCh38 human reference genome build, which includes alternate scaffolds that represent a wider variety of genomic variation in the human genome. This was quite a substantial change from previous reference genome builds, and many bioinformatics tools still haven't been updated to accept the GRCh38 reference.

The discrepancy in entry size is probably because there are alternate contigs that are uncovered by the reads that have been mapped. This is expected, but can cause a bit of confusion.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes, the reason why I use GRCh38 bed file is to have the depth values for only the karyotypic chromosomes, rather than the scaffolds and decoys. The -aa option gave the same output as the -a option. $\endgroup$
    – 719016
    Jun 8, 2017 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.