我有一个资料框(merged_COIN_plink),其中包含 181 个 RSID 编号的列,例如(rs2807367_G)-(仅显示其中两个的示例)。
ID | 表型 | rs2807367_G | rs2807376_A | 事件 | 生存 | 附言 | RS | 白细胞 |
---|---|---|---|---|---|---|---|---|
001 | -9 | 1 | 0 | 1 | 349 | 2 | 1 | 8.8 |
我正在尝试使用 Rstudio 上的生存包做一个 coxph。我已经设法为每个 RSID 单独执行此操作,但是对于 181,我想知道是否有一种方法可以并行化它或运行它,以便它自动通过每个 RSID。
这是个人 coxph 的代码:
coxph(Surv(merged_COIN_plink$SURVIVAL, merged_COIN_plink$EVENT) ~ rs2807367_G PS RS WBC, data= merged_COIN_plink)
我浏览了其他帖子,但它们似乎都令人困惑。我想知道是否可以使用诸如 RS* 之类的通配符,但我不确定这在 Rstudio 上是否容易做到。我还认为在 unix 上您可以列出 RS 编号并通过该串列运行代码以选择每个 RSID,但我不知道这在 Rstudio 上是否可行。
我试图从以 RS 开头的列名中列出一个串列,但这似乎无法正常作业:
rs_list <- merged_COIN_plink[grep("^rs",colnames(merged_COIN_plink)),]
View(rs_list)
我也不确定 for 回圈是否可以作业,但无法确定名称不同的列名。
structure(list(ID = c("100002", "100003", "100004", "100005",
"100006", "100007", "100008", "100010", "100011", "100012", "100013",
"100014", "10004", "1002", "1003", "1004", "1005", "1006", "1007",
"1008", "1010", "101001", "101002", "101003", "101004"), PHENOTYPE = c(-9L,
-9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L,
-9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L), rs2807367_G = c(1L,
0L, 2L, 2L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L,
0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L), rs34963268_C = c(1L, 1L, 1L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 1L, 0L, 2L,
0L, 0L, 1L, 0L, 0L, 0L), EVENT = c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L,
1L, 1L), SURVIVAL = c(349L, 384L, 283L, 671L, 674L, 285L, 224L,
687L, 571L, 495L, 510L, 302L, 159L, 44L, 85L, 347L, 604L, 447L,
1230L, 444L, 1260L, 758L, 392L, 379L, 188L), PS = c(2L, 0L, 0L,
0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 2L, 1L, 1L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L), RS = c(1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L,
0L), WBC = c(8.8, 8.1, 9.3, 8.9, 7.2, 6.7, 11.6, 10.7, 6.1, 12.9,
10.1, 9.1, 6.8, 13.3, 13.5, 10.9, 8.7, 11.4, 9.8, 8.9, 8, 11.3,
6, 5.6, 8.8)), row.names = c(NA, 25L), class = "data.frame")
uj5u.com热心网友回复:
for
回圈迭代当然是有效的,但是在您的(真实)资料中有这幺多列可能需要很长时间。一个可能更“优雅”的答案:)
# Get names of columns with RSID numbers (could also use `grep` or some other way
RSIDcols <- colnames(merged_COIN_plink[,3:4])
# Define formulas
formulas <- sapply(RSIDcols,
function(x) as.formula(paste('Surv(SURVIVAL, EVENT) ~ ', paste(x, " PS RS WBC"))))
# Run models
models <- lapply(formulas, function(x) {coxph(x, data = merged_COIN_plink)})
这会将每个结果输出到一个串列 ( models
) 中,然后可以通过 RSID 名称使用以下内容访问该串列:
names(models) <- RSIDcols # rename list elements to RSID name
# Access individual model results using `[[]]`
models[["rs34963268_C"]]
# > models[["rs34963268_C"]]
# Call:
# coxph(formula = x, data = merged_COIN_plink)
#
# coef exp(coef) se(coef) z p
# rs34963268_C -0.11518 0.89120 0.36030 -0.320 0.7492
# PS 0.72277 2.06012 0.43391 1.666 0.0958
# RS -1.02832 0.35761 0.54250 -1.896 0.0580
# WBC 0.05065 1.05195 0.11965 0.423 0.6721
#
# Likelihood ratio test=6.04 on 4 df, p=0.1962
# n= 25, number of events= 21
uj5u.com热心网友回复:
使用 for 回圈的解决方案应该足以满足您的需要,因为您没有太多列可以回圈并且速度似乎不太可能成为问题。我也认为很容易看到发生了什么。
这种方式将每个模型存盘在一个串列中。
# identify the columns you need
columns = colnames(dat)[grep("rs", colnames(dat))]
# create a list to store the results in
res = list()
# loop over the columns and perform the regression and save results in the list
for (i in seq_along(columns)) {
res[[i]] = coxph(Surv(SURVIVAL, EVENT) ~ get(columns[i]) PS RS WBC, data = dat)
}
然后,您可以制作一个包含您喜欢的任何值的表格,如下所示:
data.frame(nevent = sapply(res, function(x) x$nevent), columns)
nevent columns
1 21 rs2807367_G
2 21 rs34963268_C
0 评论