背景知識
講課講到批量的logrank test可以告訴我們兩個組之間的生存率差別是否顯著廷粒，收到提問說：
老師，剛才的 生存分析中红且，展示了所有基因的p值坝茎，有沒有其他列可以判斷預(yù)后是好還是差？
這是一個好問題呀暇番，如果只看p值嗤放，只是知道差別是否顯著，那么到底是基因表達(dá)量高的組預(yù)后好壁酬，還是低的組預(yù)后好次酌？
搜索發(fā)現(xiàn)我們可以通過比較兩組的中位生存時間來量化恨课。但是這個不適用于曲線有交叉的情況。
沒關(guān)系岳服，我們就只看p＜0.05的基因就可以啦剂公。
怎么得到兩組的中位生存時間呢
一段最少的代碼，用內(nèi)置數(shù)據(jù)lung來做
library(survminer)
library(survival)
fit <- survfit(Surv(time, status) ~ sex, data = lung)
summary(fit)$table[,'median']
## sex=1 sex=2
## 270 426
plot(fit)
用表達(dá)矩陣?yán)锏幕騺碜?/h4>
使用我的tinyarray包里的示例數(shù)據(jù)
library(tinyarray)
str(exprSet_hub1)
## num [1:8, 1:177] 8.03 19.13 17.63 12.06 17.64 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:8] "CXCL8" "FN1" "COL3A1" "ISG15" ...
## ..$ : chr [1:177] "TCGA-3A-A9IO-01A" "TCGA-US-A774-01A" "TCGA-HZ-A49H-01A" "TCGA-FB-A4P5-01A" ...
str(meta1)
## 'data.frame': 177 obs. of 4 variables:
## $ sample : chr "TCGA-3A-A9IO-01A" "TCGA-US-A774-01A" "TCGA-HZ-A49H-01A" "TCGA-FB-A4P5-01A" ...
## $ event : int 0 1 0 1 0 1 0 0 1 0 ...
## $ X_PATIENT: chr "TCGA-3A-A9IO" "TCGA-US-A774" "TCGA-HZ-A49H" "TCGA-FB-A4P5" ...
## $ time : int 1942 695 491 179 228 183 289 392 476 1854 ...
寫了個函數(shù)吊宋，輸入數(shù)據(jù)是整理好的表達(dá)矩陣和臨床信息纲辽，表達(dá)矩陣的列名和臨床信息表格的行名是一一對應(yīng)的。
which_better = function(exprSet_hub,meta,pvalue_cutoff = 1){
x = surv_KM(exprSet_hub,meta,pvalue_cutoff = pvalue_cutoff)
diffs = apply(exprSet_hub,1,function(g){
#g = 'ICAM1'
gene = ifelse(g>median(g),"high","low")
fit <- survfit(Surv(time, event) ~ gene, data = meta)
me = summary(fit)$table[,'median']
diff =ifelse(as.numeric( me[1]-me[2])>0,"high","low")
return(diff)
})
re = data.frame(p = x,
better_group = diffs[names(x)])
return(re)
}
which_better(exprSet_hub1,meta1)
## p better_group
## ICAM1 0.0218004 high
## FN1 0.4983822 low
## COL3A1 0.5428363 high
## MMP9 0.6977895 low
## CXCL10 0.7190313 low
## COL1A2 0.7417349 high
## ISG15 0.9598463 low
## CXCL8 0.9599565 high
如果只想看p<0.05的基因那么就加上參數(shù)：
which_better(exprSet_hub1,meta1,pvalue_cutoff = 0.05)
## p better_group
## ICAM1 0.0218004 high
等我有空璃搜，把這個函數(shù)再精裝修一下放進(jìn)tinyarray里面拖吼。
畫圖檢驗成果
p = exp_surv(exprSet_hub1,meta1)
library(patchwork)
wrap_plots(p,nrow = 2)
感謝提問的小朋友，祝大家科研順利O(∩_∩)O这吻。
