如何获得每个基因的 Cox p 值?
How to get Cox p-value for each gene?
如果你运行下面的代码,你将得到一个数据框real.dat
,其中包含 20531 个基因的 1063 个样本。有 2 个额外的列名为 time
和 event
,其中 time
是生存时间,event
在 1
和 [= 的情况下是 death
18=] 在 censored
.
的情况下
lung.dat <- read.table("genomicMatrix_lung")
lung.clin.dat <- read.delim("clinical_data_lung")
# For clinical data, get only rows which do not have NA in column "X_EVENT"
lung.no.na.dat <- lung.clin.dat[!is.na(lung.clin.dat$X_EVENT), ]
# Getting the transpose of main lung cancer data
ge <- t(lung.dat)
# Getting a vector of all the id's in the clinical data frame without any 'NA' values
keep <- lung.no.na.dat$sampleID
# getting only the samples(persons) for which we have a value rather than 'NA' values
real.dat <- ge[ge[, 1] %in% keep, ]
# adding the 2 columns from clinical data to gene expression data
keep_again <- real.dat[, 1]
temp_df <- lung.no.na.dat[lung.no.na.dat$sampleID %in% keep_again, ]
# naming the columns into our gene expression data
col_names <- ge[1, ]
colnames(real.dat) <- col_names
dd <- temp_df[, c('X_TIME_TO_EVENT', 'X_EVENT')]
real.dat <- cbind(real.dat, dd)
# renaming the 2 new added columns
colnames(real.dat)[colnames(real.dat) == 'X_TIME_TO_EVENT'] <- 'time'
colnames(real.dat)[colnames(real.dat) == 'X_EVENT'] <- 'event'
我想获取上述数据框中每个基因的单变量 Cox 回归 p 值。我怎样才能得到这个?
您可以从here下载数据。
Edit:
很抱歉没有说清楚。我已经尝试使用 survival
库中的 coxph
函数来获取它。但即使对于一个基因,它也会显示以下错误 -
> coxph(Surv(time, event) ~ HIF3A, real.dat)
Error in fitter(X, Y, strats, offset, init, control, weights = weights, :
NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights, :
Ran out of iterations and did not converge
这就是为什么我没有提供更小的可重现示例的原因。
你真的要对20531个基因中的每个基因做单变量回归吗??
疯狂猜测您的数据结构(因此根据帮助中的示例创建一个虚拟集),并猜测您要使用以下玩具示例做什么......
library("survival")
?coxph ## to see the examples
## create dummy data
test <- list(time=c(4,3,1,1,2,2,3),
event=c(1,1,1,0,1,1,0),
gene1=c(0,2,1,1,1,0,0),
gene2=c(0,0,0,0,1,1,1))
## Cox PH regression
coxph(Surv(time, event) ~ gene1, test)
coxph(Surv(time, event) ~ gene2, test)
您可能希望使用以下内容获取 CI 和更多信息。
summary(coxph(...))
希望代码的可重现性足以帮助您澄清问题
如果你运行下面的代码,你将得到一个数据框real.dat
,其中包含 20531 个基因的 1063 个样本。有 2 个额外的列名为 time
和 event
,其中 time
是生存时间,event
在 1
和 [= 的情况下是 death
18=] 在 censored
.
lung.dat <- read.table("genomicMatrix_lung")
lung.clin.dat <- read.delim("clinical_data_lung")
# For clinical data, get only rows which do not have NA in column "X_EVENT"
lung.no.na.dat <- lung.clin.dat[!is.na(lung.clin.dat$X_EVENT), ]
# Getting the transpose of main lung cancer data
ge <- t(lung.dat)
# Getting a vector of all the id's in the clinical data frame without any 'NA' values
keep <- lung.no.na.dat$sampleID
# getting only the samples(persons) for which we have a value rather than 'NA' values
real.dat <- ge[ge[, 1] %in% keep, ]
# adding the 2 columns from clinical data to gene expression data
keep_again <- real.dat[, 1]
temp_df <- lung.no.na.dat[lung.no.na.dat$sampleID %in% keep_again, ]
# naming the columns into our gene expression data
col_names <- ge[1, ]
colnames(real.dat) <- col_names
dd <- temp_df[, c('X_TIME_TO_EVENT', 'X_EVENT')]
real.dat <- cbind(real.dat, dd)
# renaming the 2 new added columns
colnames(real.dat)[colnames(real.dat) == 'X_TIME_TO_EVENT'] <- 'time'
colnames(real.dat)[colnames(real.dat) == 'X_EVENT'] <- 'event'
我想获取上述数据框中每个基因的单变量 Cox 回归 p 值。我怎样才能得到这个?
您可以从here下载数据。
Edit:
很抱歉没有说清楚。我已经尝试使用 survival
库中的 coxph
函数来获取它。但即使对于一个基因,它也会显示以下错误 -
> coxph(Surv(time, event) ~ HIF3A, real.dat)
Error in fitter(X, Y, strats, offset, init, control, weights = weights, :
NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights, :
Ran out of iterations and did not converge
这就是为什么我没有提供更小的可重现示例的原因。
你真的要对20531个基因中的每个基因做单变量回归吗??
疯狂猜测您的数据结构(因此根据帮助中的示例创建一个虚拟集),并猜测您要使用以下玩具示例做什么......
library("survival")
?coxph ## to see the examples
## create dummy data
test <- list(time=c(4,3,1,1,2,2,3),
event=c(1,1,1,0,1,1,0),
gene1=c(0,2,1,1,1,0,0),
gene2=c(0,0,0,0,1,1,1))
## Cox PH regression
coxph(Surv(time, event) ~ gene1, test)
coxph(Surv(time, event) ~ gene2, test)
您可能希望使用以下内容获取 CI 和更多信息。
summary(coxph(...))
希望代码的可重现性足以帮助您澄清问题