绘制斜率和截距的回归图(格子或ggplot2)

Question

我有一个微阵列数据集，我对其进行了limma lmFit()测试。如果你以前没听说过，这是一个强大的线性模型包，用来测试>20k基因的差异基因表达。你可以提取斜率，从模型中截取每一个基因。

我的问题是:给定一个斜率和截距值表，如何将一个地块(如果需要的话，可以选择ggplot2**'s** geom_abline**，** lattice**'s** panel.abline**，)与其相应的斜率和截距匹配？*

我的表(称为"slopeInt")的截距为第1列，斜率为第2列，并有与基因名称相对应的行名。他们的名字是这样的：

"202586_at"   "202769_at"   "203201_at"   "214970_s_at" "219155_at"

这些名称与另一个表(“数据”)中的基因名称相匹配，其中包含了一些关于我的样本的详细信息(我有24个具有不同in和时间/治疗组合的样本)和基因表达值。

它是长格式的，基因名称(如上面所示)每24行重复一次(同一基因的不同表达水平，我的每一个样本)：

ID    Time  Treatment  Gene_name    Gene_exp
...   ...   ...        ...          ...

我总共有8个我感兴趣的基因，我的Data$Gene_name中的名字与我的slopeInt表的行名相匹配。我也可以将两个表合并在一起，这不是问题。但是，我尝试了以下两种方法，通过适当的回归，为我的每一个基因提供有图表的图表，但没有结果：

使用ggplot2

ggplot(Data, aes(x = Time, y = Gene_exp, group = Time, color = Treatment)) +
facet_wrap(~ Gene_name, scales = "free_x") +
geom_point() +
geom_abline(intercept = Intercept, slope = Time), data = slopeInt) +
theme(panel.grid.major.y = element_blank())`

同时也使用Lattice

xyplot(Gene_exp ~ Time| Gene_name, Data, 
   jitter.data = T,
   panel = function(...){
     panel.xyplot(...)
     panel.abline(a = slopeInt[,1], b = slopeInt[,2])},
   layout = c(4, 2))

我在实际的geom_abline()和panel.abline()参数中尝试过多种其他方法，包括一些for循环，但我对R没有经验，我无法让它工作。我也可以有一个宽格式的数据文件(每个基因的单独列)。

任何帮助和进一步的指导将不胜感激！

下面是一个可重复的示例的一些代码：

Data <- data.frame(
 ID = rep(1:24, 8),
 Time = (rep(rep(c(1, 2, 4, 24), each = 3), 8)),
 Treatment = rep(rep(c("control", "smoking"), each = 12), 8),
 Gene_name = rep(c("202586_at", "202769_at", "203201_at", "214970_s_at",
   "219155_at", "220165_at", "224483_s_at", "227559_at"), each = 24),
 Gene_exp = rnorm(192))

slopeInt <- data.frame(
 Intercept = rnorm(8),
 Slope = rnorm(8))


row.names(slopeInt) <- c("202586_at", "202769_at", "203201_at",
"214970_s_at", "219155_at", "220165_at", "224483_s_at", "227559_at")

Answer 1

对于格子，这应该是可行的。

xyplot(Gene_exp ~ Time| Gene_name, Data, slopeInt=slopeInt,
   jitter.data = T,
   panel = function(..., slopeInt){
     panel.xyplot(...)
     grp <- trellis.last.object()$condlevels[[1]][which.packet()]
     panel.abline(a = slopeInt[grp,1], b = slopeInt[grp,2])
   },
   layout = c(4, 2)
)

在生成示例数据之前在下面的图中使用set.seed(15)

​

​

这里的“诀窍”是使用trellis.last.object()$condlevels来确定我们目前所处的条件块。然后，我们使用这些信息从我们现在通过参数传入的附加数据中提取正确的斜率信息。我认为有一种更优雅的方法来确定条件变量的当前值，但如果有，那么此时我就记不起来了。

Answer 2

如果您将Gene_name指定为slopeInt中的一个列，那么它将按照我所理解的那样工作。还请注意对ggplot调用的其他一些更改。

slopeInt$Gene_name <- rownames(slopeInt)
ggplot(Data, aes(x = Time, y = Gene_exp, color = Treatment)) +
   facet_wrap(~ Gene_name, scales = "free_x") +
   geom_point() +
   geom_abline(aes(intercept = Intercept, slope = Slope), data = slopeInt) +
   theme(panel.grid.major.y = element_blank())