Show Overall Coverage

[Firebird]

Hi,

is there a possibility to get an overview picture of the coverage? I can do this with IGV only for parts of a chromosome when I am zoomed in very deep.

Further more, is there a possiblity to filter out regions or exons which don't reach a coverage threshold?
Thanks

[Bruins]

Hi,

Are you familiar with R? I once got some R code I think from this forum (I neglected to write that down ).

I don't know what files you have available but I had a samtools pileup so I ran:

cut -f 2,8 pileup.out > pileup.cov

to extract the coverage for each position and then in R:

data <- read.table (file="pileup.cov", sep="\t", header=F)
colnames(data) <- c("position", "coverage")
depth <- mean(data[,"coverage"])
window <- 101
rangefrom <- 0
rangeto <- length(data[,"position"])
data.smoothed<-runmed(data[,"coverage"],k=window)
png(file="coverage.png",width=1900,height=1000)
plot(x=data[rangefrom:rangeto,"position"],y=data.smoothed[rangefrom:rangeto],pch=".", cex=1,xlab="bp position",ylab="depth",type="l")
dev.off()

Then check the png file for a quick coverage overview.

I did this for a small insect genome. Perhaps you should split per chromosome if human.

Hope that helps,
cheers

[Firebird]

No, unfortunatly I'm not familiar with R.

I can try it with the code u provided. Is it ok when I just change the path to my file?

Thank

[simonandrews]

SeqMonk can display an overview of a quantitation over the whole genome, for any of the quantitation types it supports. Is this anything close to what you're after?

[drio]

bruins approach does not account for regions without coverage (pileup does not show them). You need to incorporate your reference genome to account for that.

[dan]

Try Tablet:

http://bioinf.scri.ac.uk/tablet/

[Bruins]

@drio: true, but the plot can be useful to quickly see if there are peaks in the coverage and the overall depth of the coverage. An example. (sorry for the size...) It shows that there are some extreme peaks and that mostly the coverage is very low.


@Firebird: I forgot to mention the cut command is a shell command. Er...

Perhaps (if you're still interested in creating the quick plot and you have a samtools pileup file) you can split the pileup per chromosome and read in all columns. Then copy paste into R:

setwd("/path/to/your/workingdirectory")
data <- read.table (file="pileup_chr1.cov", sep="\t", header=F)

As for the path, I'm a Linux-person. If you use Windows you know better than me what the path should look like.
Now I'm not so sure about the column naming and if R accepts less column names than there are columns. If it complains about it, just add the correct number of column names after column 8 coverage.

colnames(data) <- c("chr", "position", "3", "4", "5", "6", "7", "coverage")

depth <- mean(data[,"coverage"])
window <- 101
rangefrom <- 0
rangeto <- length(data[,"position"])
data.smoothed<-runmed(data[,"coverage"],k=window)
png(file="coverage_chr1.png",width=1900,height=1000)
plot(x=data[rangefrom:rangeto,"position"],y=data.smoothed[rangefrom:rangeto],pch=".", cex=1,xlab="bp position",ylab="depth",type="l")
dev.off()

You can close R after this, it has written a png file called coverage_chr1.png.

[honey]

Hi Simon,

http://seqanswers.com/forums/showthread.php?t=8460&highlight=deep+coverage

Can it be used to show RNA-seq coverage if so what would be the steps I will be following in Seqmonk, I have tophat out put bam file.
Thanks.

[simonandrews]

Hi Simon,

http://seqanswers.com/forums/showthread.php?t=8460&highlight=deep+coverage

Can it be used to show RNA-seq coverage if so what would be the steps I will be following in Seqmonk, I have tophat out put bam file.
Thanks.

A lot of this is covered in the tutorial video we made for RNA-Seq analysis in SeqMonk, but the basic steps are:

1. Import the BAM file as single end (even if your data is paired end) and choose to split spliced reads
2. Use the feature probe generator to make probes over mRNA features, and split these into subfeatures (exons)
3. Use the base pair quantitation to quantitate your data, and ensure that you correct for probe length (so as to get comparable measures for probes of different lengths)

Your data will be now be quantitated by the read density over each exon, and you can zoom right out to see the pattern of expression across the currently selected chromosome. To see the pattern over the whole genome simply select the dataset you want to see in the data view (the folders in the top left), and the genome view in the top right will show you a coverage plot for the whole genome (similar to the one attached to an earlier note in this thread).

Hope this helps