This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Samtools depth to get the read depth

I am using samtools depth to get read depth at each position from the bam file which I have after doing alignment. I used the command: samtools depth -a -H input.bam -o output_duplicates

Since my bam file comprise more than one sample, how samtools depth will calculate the read depth? The output file I have consist of three columns with the first column chromosome, second position and third read depth.

Thank you in advance

sequencing rna-seq
• 4,305 views • link
written by bandanaschapagain • bandanaschapagain

What do you see for the following?

$ samtools --version
• 0 views • link
written by GenoMax GenoMax

Hi @GenoMax samtools that I have is of 1.3.1 version. I see: samtools 1.3.1 Using htslib 1.3.1

• 0 views • link
written by bandanaschapagain • bandanaschapagain

Hi,

Knowing that -H option is not so important as it is for header, I just ran the code without it. for RG in SAMPLE; do samtools view -u -r ${RG} input.bam | samtools depth -a -o {RG}.depth.txt; done

Now the error I have is: depth: invalid option -- 'o'

Is it because of the old version of samtools I am using?

• 0 views • link
written by bandanaschapagain • bandanaschapagain

Is it because of the old version of samtools I am using?

yes

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

you can always redirect the output instead of using -o

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Install latest samtools using conda if you are not able to upgrade the main install.

• 0 views • link
written by GenoMax GenoMax

-I did install latest version of samtools: 

samtools --version=  samtools 1.11
Using htslib 1.11

Then I tried to split the original .bam file by using one sample

samtools view -h -r Sampl1 input.bam > Sample1.bam

Then I used 

samtools depth -a -H  Sampl1.bam -o  Sample1_bam_depth

But I see that in Sample1_bam_depth it's only printing header as 

#CHROM POS Sample1.bam

I don't understand what I am missing at this point.

• 0 views • link
updated by GenoMax GenoMax • written by bandanaschapagain • bandanaschapagain

Is there data in Sample1.bam aside from the header? Is your read-group called Sampl1 in your file.

• 0 views • link
written by GenoMax GenoMax

I found the error. It was the RG I was missing up, and it just printed header to my bam file. Now, I am able to run samtools depth. Thank you.

• 0 views • link
written by bandanaschapagain • bandanaschapagain

1 answer

use option -r or -R of samtools view

  -r STR   only include reads in read group STR [null]
  -R FILE  only include reads with read group listed in FILE [null]

and pipe into samtools depth

for RG in RGSAMPLE1 RGSAMPLE2 RGSAMPLE1 ; do samtools view -u -r ${RG} in.bam | samtools depth -a -H -o {RG}.depth.txt  ; done
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Hi @Pierre Lindenbaum I did as per your suggetion, and I got the error as depth: invalid option -- 'H' So I changed my code to parse just one sample. But I got the same error.

My code for one sample is: for RG in Sample1; do samtools view -u -r ${RG} input.bam | samtools depth -a -H -o {RG}.depth.txt; done

With the code samtools view -u -r Sample1 input.bam > Sample1.bam the ouput file is not empty. I am really stuck at this point.

Thank you for your help.

• 0 views • link
written by bandanaschapagain • bandanaschapagain

got the error as depth: invalid option -- 'H'

your version samtools is too old. check your version number.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this