---

微信公眾號：http://blog.genesino.com
生存分析指根據(jù)試驗或調(diào)查得到的數(shù)據(jù)對生物或人的生存時間進(jìn)行分析和推斷，研究生存時間和結(jié)局與眾多影響因素間關(guān)系及其程度大小的方法，也稱生存率分析或存活率分析。常用于腫瘤等疾病的標(biāo)志物篩選、療效及預(yù)后的考核。

簡單地說，比較兩組或多組人群隨著時間的延續(xù)，存活個體的比例變化趨勢。活著的個體越少的組危險性越大，對應(yīng)的基因?qū)膊∮绊懺酱?#xff0c;對應(yīng)的藥物治療效果越差。

生存分析適合于處理時間-事件數(shù)據(jù)，如下

生存時間數(shù)據(jù)有兩種類型：

• 完全數(shù)據(jù) (complete data)指被觀測對象從觀察起點到出現(xiàn)終點事件所經(jīng)歷的時間; 一般用狀態(tài)值1或TRUE表示。
• 截尾數(shù)據(jù) (consored data)或刪失數(shù)據(jù)，指在出現(xiàn)終點事件前，被觀測對象的觀測過程終止了。由于被觀測對象所提供的信息是不完全的，只知道他們的生存事件超過了截尾時間。截尾主要由于失訪、退出和終止產(chǎn)生。一般用狀態(tài)值0或FALSE表示。
• TCGA中的臨床數(shù)據(jù)標(biāo)記也符合這個規(guī)律，在下面軟件運(yùn)行時也可修改狀態(tài)值的含義, 但一般遵循這個規(guī)律。

生存概率 (survival probability)指某段時間開始時存活的個體至該時間結(jié)束時仍然存活的可能性大小。

生存概率=某人群活過某段時間例數(shù)/該人群同時間段期初觀察例數(shù)。

生存率 (Survival rate)，用S(t)表示，指經(jīng)歷t個單位時間后仍存活的概率，若無刪失數(shù)據(jù)，則為活過了t時刻仍然存活的例數(shù)/觀察開始的總例數(shù)。如果有刪失數(shù)據(jù)，分母則需要按時段進(jìn)行校正。

生存分析一個常用的方法是壽命表法。

壽命表是描述一段時間內(nèi)生存狀況、終點事件和生存概率的表格，需計算累積生存概率即每一步生存概率的乘積 (也可能是原始生存概率)，可完成對病例隨訪資料在任意指定時點的生存狀況評價。

R做生存分析

R中做生存分析需要用到包survival和survminer。輸入數(shù)據(jù)至少兩列，存活時間和生存狀態(tài)，也就是測試數(shù)據(jù)中的Days.survial和vital_status列。如果需要比較不同組之間的差異，也需要提供個體的分組信息，如測試數(shù)據(jù)中的PAM50列。對應(yīng)TCGA的數(shù)據(jù)，一般根據(jù)某個基因的表達(dá)量或突變有無對個體進(jìn)行分組。

讀入數(shù)據(jù)

library(survival) BRCA <- read.table('BRCA.tsv', sep="\t", header=T) head(BRCA) ID SampleType PAM50Call_RNAseq Days.survival pathologic_stage 1 TCGA-E9-A2JT-01 Tumor_type LumA 288 stage iia 2 TCGA-BH-A0W4-01 Tumor_type LumA 759 stage iia 3 TCGA-BH-A0B5-01 Tumor_type LumA 2136 stage iiia 4 TCGA-AC-A3TM-01 Tumor_type Unknown 762 stage iiia 5 TCGA-E9-A5FL-01 Tumor_type Unknown 24 stage iib 6 TCGA-AC-A3TN-01 Tumor_type Unknown 456 stage iibvital_status 1 0 2 0 3 0 4 0 5 0 6 0

簡單地看下每一列都有什么內(nèi)容，方便對數(shù)據(jù)整體有個了解，比如有無特殊值。

summary(BRCA) ID SampleType PAM50Call_RNAseq Days.survival TCGA-3C-AAAU-01: 1 Tumor_type:1090 Basal :138 Min. : 0.0 TCGA-3C-AALI-01: 1 Her2 : 65 1st Qu.: 450.2 TCGA-3C-AALJ-01: 1 LumA :415 Median : 848.0 TCGA-3C-AALK-01: 1 LumB :194 Mean :1247.0 TCGA-4H-AAAK-01: 1 Normal : 24 3rd Qu.:1682.8 TCGA-5L-AAT0-01: 1 Unknown:254 Max. :8605.0 (Other) :1084 pathologic_stage vital_status stage iia :359 Min. :0.0000 stage iib :259 1st Qu.:0.0000 stage iiia:156 Median :0.0000 stage i : 90 Mean :0.1394 stage ia : 85 3rd Qu.:0.0000 stage iiic: 67 Max. :1.0000 (Other) : 74

計算壽命表

# Days.survival：跟蹤到的存活時間 # vital_status: 跟蹤到的存活狀態(tài) # ~1表示不進(jìn)行分組 fit <- survfit(Surv(Days.survival, vital_status)~1, data=BRCA)# 獲得的survial列就是生存率 summary(fit)Call: survfit(formula = Surv(Days.survival, vital_status) ~ 1, data = BRCA)time n.risk n.event survival std.err lower 95% CI upper 95% CI116 1021 1 0.999 0.000979 0.997 1.000158 1017 1 0.998 0.001386 0.995 1.000160 1016 1 0.997 0.001697 0.994 1.000172 1010 1 0.996 0.001962 0.992 1.000174 1008 1 0.995 0.002195 0.991 0.999197 1003 1 0.994 0.002406 0.989 0.999224 993 1 0.993 0.002604 0.988 0.998227 990 1 0.992 0.002788 0.987 0.998239 987 1 0.991 0.002961 0.985 0.997255 981 1 0.990 0.003125 0.984 0.996266 978 1 0.989 0.003282 0.983 0.996295 965 1 0.988 0.003435 0.981 0.995302 962 1 0.987 0.003581 0.980 0.994304 958 1 0.986 0.003723 0.979 0.993320 948 1 0.985 0.003862 0.977 0.993322 946 1 0.984 0.003995 0.976 0.992

繪制生存曲線，橫軸表示生存時間，縱軸表示生存概率，為一條梯形下降的曲線。下降幅度越大，表示生存率越低或生存時間越短。

library(survminer) # conf.int：是否顯示置信區(qū)間 # risk.table: 對應(yīng)時間存活個體總結(jié)表格 ggsurvplot(fit, conf.int=T,risk.table=T)

PAM50是通過50個基因的表達(dá)量把乳腺癌分為四種類型 (Luminal A, Luminal B, HER2-enriched, and Basal-like)作為預(yù)后的標(biāo)志。根據(jù)PAM50屬性對病人進(jìn)行分組，評估比較兩組之間生存率的差別。

# 這三步不是必須的，只是為了方便，選擇其中的4個確定了的分組進(jìn)行分析 # 同時為了簡化圖例，給列重命名一下，使得列名不那么長 BRCA_PAM50 <- BRCA[grepl("Basal|Her2|LumA|LumB",BRCA$PAM50Call_RNAseq),] BRCA_PAM50 <- droplevels(BRCA_PAM50) colnames(BRCA_PAM50)[colnames(BRCA_PAM50)=="PAM50Call_RNAseq"] <- 'PAM50' # 按PAM50分組 fit <- survfit(Surv(Days.survival, vital_status)~PAM50, data=BRCA_PAM50) # 繪制曲線 ggsurvplot(fit, conf.int=F,risk.table=T, risk.table.col="strata", pval=T)

簡化Stage信息，先只查看大的階段

BRCA_PAM50$pathologic_stage <- gsub('(i+v*).*', "\\1", BRCA_PAM50$pathologic_stage) BRCA_PAM50$pathologic_stage <- as.factor(BRCA_PAM50$pathologic_stage) colnames(BRCA_PAM50)[colnames(BRCA_PAM50)=="pathologic_stage"] <- 'PS' fit <- survfit(Surv(Days.survival, vital_status)~PS, data=BRCA_PAM50) # 繪制曲線 ggsurvplot(fit, conf.int=F,risk.table=T, risk.table.col="strata", pval=T)

參考資料

• http://rpubs.com/xuefliang/153247
• http://www.sthda.com/english/wiki/survminer-r-package-survival-data-analysis-and-visualization

總結(jié)

以上是生活随笔為你收集整理的R语言学习 - 非参数法生存分析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。