Computing GC content and smoothed coverage in R

Edit: pre-processing the data to compute GC content and compute running mean with flexible window size. The genome sequence is required to compute the GC content. I assume it is present as DNAString object, use read.DNAStringSet to read in the genome from a fasta file or use the BSGenome package.

You will have to experiment a bit how to best load the data, e.g. maybe process only one chromosome at a time will help save memory.

library(Biostrings)

window.size <- 1000 # configure the window size

# create some test yeast data, the coordinates must fit with the coverage ofc, 
# so this might not work directly:
data(yeastSEQCHR1)
chrom <- DNAString(yeastSEQCHR1)

# compute gc content in a sliding window (as a fraction, if you want % : *100):
gc <- rowSums(letterFrequencyInSlidingView(chrom, window.size, c("G","C")))/window.size


# the coverage vector can be shorter than the sequence. 
# The length is equal to the largest end position of an alignment in the data:
pad.right <- function (rle, value=0, len) {
if (length(rle) > len) stop("length mismatch: len must be >= length(rle)!")
    if (len > length(rle)) 
      c(rle, Rle(value, len-length(rle))) 
    else rle
}
cvg <- pad.right(cvg, len=length(gc))

# we have to process the coverage vector to have the same window size:
cvgw <- runmean(x=cvg, k=window.size, endrule = c("constant")) 
# endrule is important to get a vector of same size, see ?runmean

Finally, do a scatter-plot and maybe add a lowess scatterplot smoothing curve to visualize the correlation:

plot(x=gc, y=cvgw)
lines(lowess(x=gc, y=cvgw), col=2)

should do the trick.

Plotting the GC content as a genome track

Once you are using R, the easiest way to plot coverage is possibly the GenomeGraphs package. Use the function makeBaseTrack to generate a coverage track. Given your coverage is in the cvg Rle object and gc is an Rle with the gc content, that looks like this:

library("GenomeGraphs")
baseCoverage <- as.numeric(cvg) # makeBaseTrack doesn't take Rle
gcContent <- as.numeric(gc)
 # contruct your gene models here: 
gene <- ...
pl <- list( 
      genes=gene, makeGenomeAxis(T,T,T),
      coverage=makeBaseTrack(base=1:length(baseCoverage), value=baseCoverage),
      gc=makeBaseTrack(base=1:length(gcContent), value=gcContent)
)


gdPlot(pl, min=1, max=1000) # restrict viewport

There are also multiple configurable functions in that package to plot the gene models, e.g.: makeAnnotationTrack or using biomaRt (from the package documentation):

mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# Next we can retrieve the gene structure of the gene of interest.
gene <- makeGene(id = "ENSG00000095203",
type = "ensembl_gene_id", biomart = mart)

Sorry if this is a silly suggestion, but have you looked at the data using the UCSC browser? It will easily present the two tracks for any given region. The browser can present BAM files (here) and there is the GC content track