TCGA数据库：生存分析
本文介绍生存分析，其实，在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的数据处理后，能否做某基因突变与野生型的生存分析呢？其实都是一样的道理，其他的也是一样。照葫芦画瓢而已，大家自己去试试。
腾讯云开发者
