Customize IGV coverage track in R

Question

[lkmklsmn]

Hi Biostars,

I want to recreate the IGV coverage track in R. I want to be able to customize and edit a coverage plot. All I need for this is per base coverage of a given gene. I have bam input file.

Thanks

Answer 1

[poisonAlien]

Okay. Is this what you want? Check out Gviz Bioconductor package specifically section 4.9: Alignment track.

For ex:

library(Gviz) #load library
alTrack=AlignmentsTrack("foo_sorted.bam",isPaired=T) #Read bam file
plotTracks(alTrack,from=87072069,to=87073068,chromosome="chr4") #plot the required region; plots coverage, alignment section
plotTracks(alTrack,from=87072069,to=87073068,chromosome="chr4",type="coverage") #plot only coverage, no depth or alignment section

Answer 2

[Caddymob]

GATK DepthOfCoverage or Bedtools genomecov are probably your best bet

[lkmklsmn]

I used this:

coverageBed -abam test.bam -b test.bed -d -split -hist > coverageTest.txt

I can generate a coverage track from this in R but the coverage quantity does not agree with what I see on IGV. IGVs numbers are lower.

Answer 3

[Varun Gupta]

Try this

genomeCoverageBed -ibam test.bam -g my.genome -d > perbase_cov.txt

where my.genome is the genomic file.

Genome files must be tab-delimited and are structured as follows (this is an example for C. elegans):

chrI 15072421
chrII 15279323
...
chrX 17718854
chrM 13794

2nd column is the length of the chromosome.

Refer bedtools

Varun

Answer 4

[dariober]

Since you have "RNA-seq" among your tags I guess this is what you want to plot. In this case I would use bedtools genomeCovergeBed as others have suggested, but with the -split option to correctly represent reads spanning different exons. From the docs:

bedtools genomecov will, by default, screen for overlaps against the entire span of a spliced/split BAM alignment or blocked BED12 feature. When dealing with RNA-seq reads, for example, one typically wants to only screen for overlaps for the portions of the reads that come from exons (and ignore the interstitial intron sequence). The -split command allows for such overlaps to be performed

Answer 5

[lkmklsmn]

Thanks for the responses. I actually ended up using:

samtools depth file.bam > coverage.bam

This command gave me the exact numbers I was seeing in IGV on a per-nucleotide basis

[lkmklsmn]

Maybe a little more explanation since I didnt really end up using R.

Here is my pseudocode form within R:

system("samtools depth file.bam > coverage.bam")
tab<-read.table("coverage.bam")
barplot(tab[,3])

[dariober]

Take care that samtools depth reports coverage capped to max 8000 I think. Also if you have RNA-Seq make sure reads spanning introns are not counted within the intron.