Long vector-plot/Coverage plot in R

Question

I really need your R skills here. Been working with this plot for several days now. I'm a R newbie, so that might explain it.

I have sequence coverage data for chromosomes (basically a value for each position along the length of every chromosome, making the length of the vectors many millions). I want to make a nice coverage plot of my reads. This is what I got so far:

Looks alright, but I'm missing y-labels so I can tell which chromosome it is, and also I've been having trouble modifying the x-axis, so it ends where the coverage ends. Additionally, my own data is much much bigger, making this plot in particular take extremely long time. Which is why I tried this HilbertVis plotLongVector. It works but I can't figure out how to modify it, the x-axis, the labels, how to make the y-axis logged, and the vectors all get the same length on the plot even though they are not equally long.

source("http://bioconductor.org/biocLite.R")
biocLite("HilbertVis")
library(HilbertVis)
chr1 <- abs(makeRandomTestData(len=1.3e+07)) 
chr2 <- abs(makeRandomTestData(len=1e+07)) 

par(mfcol=c(8, 1), mar=c(1, 1, 1, 1), ylog=T)

# 1st way of trying with some code I found on stackoverflow
# Chr1
plotCoverage <- function(chr1, start, end) { # Defines coverage plotting function.
  plot.new()
  plot.window(c(start, length(chr1)), c(0, 10))
  axis(1, labels=F) 
  axis(4)
  lines(start:end, log(chr1[start:end]), type="l")
}
plotCoverage(chr1, start=1, end=length(chr1)) # Plots coverage result.

# Chr2
plotCoverage <- function(chr2, start, end) { # Defines coverage plotting function.
  plot.new()
  plot.window(c(start, length(chr1)), c(0, 10))
  axis(1, labels=F) 
  axis(4)
  lines(start:end, log(chr2[start:end]), type="l")
}
plotCoverage(chr2, start=1, end=length(chr2)) # Plots coverage result.


# 2nd way of trying with plotLongVector
plotLongVector(chr1, bty="n", ylab="Chr1") # ylab doesn't work
plotLongVector(chr2, bty="n")

Then I have another vector called genes that are of special interest. They are about the same length as the chromosome-vectors but in my data they contain more zeroes than values.

genes_chr1 <- abs(makeRandomTestData(len=1.3e+07)) 
genes_chr2 <- abs(makeRandomTestData(len=1e+07)) 

These gene vectors I would like plotted as a red dot under the chromosomes! Basically, if the vector has a value there (>0), it is presented as a dot (or line) under the long vector plot. This I have not idea how to add! But it seems fairly straightforward.

Please help me! Thank you so much.

Answer 1

DISCLAIMER: Please do not simply copy and paste this code to run off the entire positions of your chromosome. Please sample positions (for example, as @Gx1sptDTDa shows) and plot those. Otherwise you'd probably get a huge black filled rectangle after many many hours, if your computer survives the drain.

Using ggplot2, this is really easily achieved using geom_area. Here, I've generated some random data for three chromosomes with 300 positions, just to show an example. You can build up on this, I hope.

# construct a test data with 3 chromosomes and 100 positions
# and random coverage between 0 and 500
set.seed(45)
chr <- rep(paste0("chr", 1:3), each=100)
pos <- rep(1:100, 3)
cov <- sample(0:500, 300)
df  <- data.frame(chr, pos, cov)

require(ggplot2)
p <- ggplot(data = df, aes(x=pos, y=cov)) + geom_area(aes(fill=chr))
p + facet_wrap(~ chr, ncol=1)

Answer 2

You could use the ggplot2 package.

I'm not sure what exactly you want, but here's what I did: This has 7000 random data points (about double the amount of genes on Chromosome 1 in reality). I used alpha to show dense areas (not many here, as it's random data).

library(ggplot2)
Chr1_cov <- sample(1.3e+07,7000)
Chr1 <- data.frame(Cov=Chr1_cov,fil=1)
pl <- qplot(Cov,fil,data=Chr1,geom="pointrange",ymin=0,ymax=1.1,xlab="Chromosome 1",ylab="-",alpha=I(1/50))
print(pl)

And that's it. This ran in less than a second. ggplot2 has a humongous amount of settings, so just try some out. Use facets to create multiple graphs.

---

The code beneath is for a sort of moving average, and then plotting the output of that. It is not a real moving average, as a real moving average would have (almost) the same amount of data points as the original - it will only make the data smoother. This code, however, takes an average for every n points. It will of course run quite a bit faster, but you will loose a lot of detailed information.

VeryLongVector <- sample(500,1e+07,replace=TRUE)

movAv <- function(vector,n){
    chops <- as.integer(length(vector)/n)
    count <- 0
    pos <- 0
    Cov <-0
    pos[1:chops] <- 0
    Cov[1:chops] <- 0
    for(c in 1:chops){
        tmpcount <- count + n
        tmppos <- median(count:tmpcount)
        tmpCov <- mean(vector[count:tmpcount])
        pos[c] <- tmppos
        Cov[c] <- tmpCov
        count <- count + n
    }

    result <- data.frame(pos=pos,cov=Cov)
    return(result)
}

Chr1 <- movAv(VeryLongVector,10000)
qplot(pos,cov,data=Chr1,geom="line")