2

我真的需要你的 R 技能。几天来一直在处理这个情节。我是 R 新手，所以这可以解释它。

我有染色体的序列覆盖数据（基本上是每个染色体长度上每个位置的值，使向量的长度达到数百万）。我想为我的阅读制作一个很好的覆盖图。这是我到目前为止得到的： 在此处输入图像描述

看起来不错，但我缺少 y 标签，所以我可以知道它是哪条染色体，而且我在修改 x 轴时遇到了麻烦，所以它在覆盖范围结束的地方结束。此外，我自己的数据要大得多，这使得这个情节特别需要很长时间。这就是我尝试这个 HilbertVis plotLongVector 的原因。它可以工作，但我不知道如何修改它、x 轴、标签、如何记录 y 轴以及向量在图上的长度都相同，即使它们的长度不一样。

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")

然后我有另一个特别感兴趣的称为基因的载体。它们与染色体向量的长度大致相同，但在我的数据中，它们包含的零比值多。

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

我想将这些基因载体绘制为染色体下方的一个红点！基本上，如果向量在那里有一个值 (>0)，它会在长矢量图下显示为一个点（或线）。这个我不知道如何添加！但这似乎相当简单。

请帮我！太感谢了。

4

2 回答 2

4

免责声明：请不要简单地复制和粘贴此代码来运行您染色体的整个位置。请采样位置（例如，如@Gx1sptDTDa 所示）并绘制它们。否则，如果您的计算机在排水管中幸存下来，您可能会在许多小时后得到一个巨大的黑色填充矩形。

使用ggplot2，使用 确实很容易实现geom_area。在这里，我为 300 个位置的三个染色体生成了一些随机数据，只是为了展示一个例子。我希望你可以在此基础上再接再厉。

# 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)

ggplot2_geom_area_coverage_plot

于 2013-01-31T17:38:15.513 回答
1

您可以使用 ggplot2 包。

我不确定你到底想要什么，但这就是我所做的： 在此处输入图像描述 这有 7000 个随机数据点（实际上是染色体 1 上基因数量的两倍）。我使用 alpha 来显示密集区域（这里不多，因为它是随机数据）。

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)

就是这样。这运行不到一秒钟。ggplot2 有大量的设置，所以尝试一下。使用构面创建多个图表。

下面的代码是一种移动平均线，然后绘制它的输出。它不是真正的移动平均线，因为真正的移动平均线将具有（几乎）与原始移动平均线相同数量的数据点 - 它只会使数据更平滑。但是，此代码对每 n 个点取平均值。它当然会运行得更快一些，但你会丢失很多详细信息。

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")

在此处输入图像描述

于 2013-01-31T16:38:51.340 回答