本文介绍生存分析，其实，在R中，生存分析很简单，大家在网上能找到无数的文章。利用survival包就可以。就是按照下列公式就可以完成简单的生存分析。

fit <- survfit(Surv(生存时间, 生存状态) ~ 分组, data=数据框)

我们这里就结合基因的表达量，来进行分析。

首先加载我们的数据。

options(stringsAsFactors = F)
#加载表达数据
load("F:/TCGA/HTSeq-FPKM/Rdata/data/TCGA-COAD-Exp.Rdata")
#加载临床数据
load("K:/TCGA/clinicalData/tidyAllCancerData/TCGA-COAD -Clindata.Rdata")

表达数据和临床数据，我之前已经上传到网盘

之前处理后的数据进行简单的处理，其实就是去掉正常组织的样本，再把列名变成3个字段。因为原来表达矩阵中病人的barcode长，"TCGA-AA-3662-11A-01R-1723-07"，而临床数据中的只有前3段。具体关于barcode这个在前面有介绍，也可以参考文章TCGAbiolinks包介绍，里面有详细介绍。

##提取肿瘤样本的表达矩阵
exp <- transomeData[["allGenesExp"]][["Exp"]]
tumExp <- exp[,transomeData[["allGenesExp"]][["TumorBarcode"]]]
barcode <- gsub("(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3",colnames(tumExp))
colnames(tumExp) <- barcode 
##

同样我们也要处理一下临床数据，我们之前处理的临床数据是这样的：

我们这里也需要简单处理一下。

#提取临床数据
library(dplyr)
clin <- clinData
clin <- filter(clin,follow != "--")
clin <- clin[,c("sample","survivalState","follow")]
clin[1,3]
clin$follow <- as.numeric(clin$follow)/365

其实，就是只要生存时间和生存状态这2个数据，再删除缺失值。follow除了365，单位就是年啦。

然后我们将表达矩阵与临床数据融合，因为不是每个病人的数据都是一一对应的，简单说，就是病人有表达数据，但他的临床数据就不全，我们也删除了缺失值的病人的临床数据，所以我们只需要具有临床数据又有表达数据的病人的数据。也就是取一个交集。

##融合数据
interbarcode <- intersect(barcode,clin$sample)
gene <- "LLGL2"
geneExp <- tumExp[gene,interbarcode]
geneExp <- t(geneExp) %>% as.data.frame() %>% data.matrix() %>% as.data.frame() 
geneExp$sample <- rownames(geneExp)
mergdata <- merge(clin,geneExp,by = "sample")

这里的gene我只写了一个，所以融合后是下面这样的数据。你也可以把所有基因都融合进去，这里是案例，就演示了一个。

接下来就是添加一个分组信息。

mergdata$Group <- ifelse(mergdata[,gene] > median(mergdata[,gene]),"High","Low")

我们以表达值的中位数为分界线，高于中位值为高表达，低于或等于中位值为低表达。也可以用均值。

得到上面这样的数据后，我们就可以按照刚刚的公式进行生存分析了：

######################### 生存分析
library(survival)
library(survminer)
fit <- survfit(Surv(follow, survivalState) ~ Group, data=mergdata)

绘图的话就用ggsurvplot函数。

ggsurvplot(fit,data = mergdata)

好像不是很美观，我们可以调整一下参数，比如y轴下部分很空，我们可以调整一下y轴坐标。

这样看着就好很多啦。我们在进行其他参数的调整。

ggsurvplot(fit, main = "Survival curve",
           font.main = c(16, "bold", "darkblue"),
           font.tickslab = c(12, "plain", "darkgreen"),
           legend.title = gene,
           legend.labs = c("High", "Low"),
           ylim = c(0.75,1),
           # Add p-value and tervals
           pval = TRUE,
           palette = c("#FF0000", "#2200FF"),
           ggtheme = theme_bw() 
)

如果我们要一次批量分析很多基因的高低表达与生存的关系，写一个循环，批量绘图了。

尽管本文是介绍基因表达量的生存分析，但其他的也是一样，就看你怎么分组，比如我们前面介绍SNP的数据处理后，能否做某基因突变与野生型的生存分析呢？其实都是一样的道理，其他的也是一样。照葫芦画瓢而已，大家自己去试试。