Samtools depth to get the read depth
1
0
Entering edit mode
3.2 years ago

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 • 2.4k views
ADD COMMENT
0
Entering edit mode

Hi, I highly appreciate any help at this point.

ADD REPLY
0
Entering edit mode

What do you see for the following?

$ samtools --version
ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

I have is of 1.3.1 version

yes it's too old. I added the option -H in https://github.com/samtools/samtools/pull/937 sept 2018.

while 1.3.1 is 2016 https://github.com/samtools/samtools/releases/tag/1.3.1

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
1
Entering edit mode

you can always redirect the output instead of using -o

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

-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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

yes

ADD REPLY
1
Entering edit mode
3.2 years ago

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
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

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

ADD REPLY

Login before adding your answer.

Traffic: 1960 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6