Question: How To Get The Read Depth?
6
gravatar for madkitty
7.5 years ago by
madkitty590
Canada
madkitty590 wrote:

My boss said he wanted to get the read depth graph for each CNV we found. The only data I have is a file with the number of reads aligned to each position (let's say from position 1 to 500 on the genome) Let's say we have a CNV deletion identified from position 100 to 150 .. How am I suppose to output the read depth for that CNV ?

ADD COMMENTlink written 7.5 years ago by madkitty590

Thanks a lot for your answers guys! All of them are really useful :)

ADD REPLYlink written 7.5 years ago by madkitty590
3
gravatar for Stefano Berri
7.5 years ago by
Stefano Berri4.1k
Cambridge, UK
Stefano Berri4.1k wrote:

If you have a sorted and indexed bam file you could do

samtools view yourFile.bam chr1:1000-5000 | wc

this will give you the number of reads mapping over chromosome 1 from position 1000 to 5000. For comparison you should also do the same with the normal and, at least, normalise for total number of reads.

However, how did you figure out that you have a CNV in that region? What program did you use? Can't the program you used tell you?

also consider this is very crude, as the number of reads have noise.

ADD COMMENTlink written 7.5 years ago by Stefano Berri4.1k

Indeed, as your CNV calling algorithm probably used read-depth as part of its protocol, you should really be sure you have a good handle on what it is doing as well as what you expect to plot. That read depth is usually on a per nucleotide basis but you can also use the number of reads that overlap with your called interval as a proxy.

ADD REPLYlink written 7.5 years ago by DG7.1k

Thanks for your answer. I used Breakdancer to find CNVs, Breakdancer already gives me the number of reads aligned to the CNV found, but my boss doesn't believe breakdancer and wants to see a figure similar to the one posted by Chris Miller below...

ADD REPLYlink written 7.5 years ago by madkitty590
3
gravatar for Chris Miller
7.5 years ago by
Chris Miller21k
Washington University in St. Louis, MO
Chris Miller21k wrote:

You'll probably want to plot the read depth in small bins, so that you can visualize the changes more clearly. I'd recommend grabbing the CNV region, plus some flanks with samtools, then writing a little script that parses the output and returns the average depth of coverage in 100bp bins. Plot these values for the tumor and the normal across the region, and if all goes well, you'll see something like this:

enter image description here

Normal is on top, tumor is on bottom, and you can clearly see the deletion in the tumor.

ADD COMMENTlink written 7.5 years ago by Chris Miller21k

that's really nice! Exactly what I need, but for the script can I just use a package that already exists in R ? Are you aware of any ?

ADD REPLYlink written 7.5 years ago by madkitty590
1

For a simple implementation in python, try ngCGH. In R, see the GenomicRanges and Shortread packages and the coverage() method.

ADD REPLYlink written 7.5 years ago by Sean Davis25k
2
gravatar for Nikolay Vyahhi
7.5 years ago by
Nikolay Vyahhi1.3k
St. Petersburg, Russia
Nikolay Vyahhi1.3k wrote:

AFAIK read depth for CNV is coverage for this region in alignment.

E.g. if you identified that there is no gene/region (deletion), then it's depth = 0.

E.g. if you identified 2 copies of the same region, then depth will be about twice as usual depth (given uniformity of coverage).

ADD COMMENTlink written 7.5 years ago by Nikolay Vyahhi1.3k
1
gravatar for Sean Davis
7.5 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:

You might try:

samtools depth -r chr1:100-150 bamfile.bam

I agree with Stefano, though, that you probably need to understand what your CNV caller did before proceeding.

ADD COMMENTlink written 7.5 years ago by Sean Davis25k
1
gravatar for Ian
7.5 years ago by
Ian5.6k
University of Manchester, UK
Ian5.6k wrote:

If you can get a BED file (chr start end) of CNV regions you can use bedtools to get coverage of reads for each CNV region.

bedtools coverage -abam reads.bam -b cnv_regions.bed
ADD COMMENTlink written 7.5 years ago by Ian5.6k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1771 users visited in the last hour