人人都可以学会生存分析（学徒数据挖掘）

下面是《生信入门第6期》学员的分享

她上一个笔记是：学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？

这次带来的是2个基因组合分群的生存分析，如下：

学徒作业，完成生存分析；

• 数据集是：METABRIC dataset (Curtis et al., 2012).
• 生存结局事件是：Disease-specific (DSS) survival
• 表达量分组策略是：
  • BCL11A High + MBNL1 Low (n=148)
  • BCL11A Low + MBNL1 High (n=345)
  • Intermediate (n=1478).

首先看看Disease-specific (DSS) survival的含义

OS 总体生存期：Overall Survival

• 定义： 结局指标是死亡时间，这个死亡是任何原因导致的死亡都算进去，只关心是否死亡，不关心因为何种原因死亡。
• 优点 ：能比较方便的记录，因为患者死亡的日期确认没有困难。只要研究结果显示生存有提高，就可认为是是临床又获益。
• 缺点：随访的时间较长

PFS 无进展生存期：Progression Free Survival

• 定义："the length of time during and after the treatment of a disease, such as cancer, that a patient lives with the disease but it does not get worse" ；指的是疾病经过治疗后，没有进一步恶化的生存期，结局指标是 发生恶化或 死亡
• **优点：**增加了发生恶化这一结局指标节点，随访时间短一些，对应的改善是未恶化与未死亡，可以反映临床获益。如果 PFS提高了，可以认为临床有获益。
• **缺点：**因为增加了 发生恶化这一结局指标，我们就要问一个问题了，**何为发生恶化？有没有明确的标准？**相对于记录是否死亡，判断病人是否病情恶化的难度要大得多，因此这就要求对 **发生恶化的标准进行明确的定义。**发生恶化的定义通常涉及影像学资料（普通X线，CT扫描，MRI，PET扫描，超声）或其他方面：生化进展可以根据肿瘤标志物的增加。

DFS 无病生存期: Disease Free Survival

• 定义：经过治疗后未发现肿瘤，结局指标为 疾病复发或 死亡，同样不关心死亡原因。The measure of time after treatment during which no sign of cancer is found. 同样 Relapse Free survival的定义也类似。
• 优点：是临床获益的重要反映，随访时间可以缩短，因为增加了疾病复发这一节点。没有复发或没有死亡可以反映临床获益。
• **缺点：**同理，因为增加了 疾病复发这一节点，我们就要问了， **何为复发？如何明确有无疾病复发？**对于记录死亡，明确是否复发的难度要大得多。记录比较困难。

DSS 疾病特异性生存期: Disease Free Survival

• 定义：结局指标改变为 由特定疾病导致的死亡，这时候开始关心死亡的原因是否是由特定疾病导致的。如果不是特定疾病导致的则不计入结局指标。
• 优点：针对性的反应临床获益，DSS提升能够很好的反应特定疾病的临床获益，特定疾病导致的死亡减少或增加。
• 缺点：同样很明显，相比简单确认患者是否死亡，这时候我们需要明确 何为由疾病导致的死亡？有没有明确的标准？，这个问题需要一个专业的判断。患者的死因经常并不容易明确。

其次了解一下METABRIC dataset数据库的背景知识

网址：http://molonc.bccrc.ca/aparicio-lab/research/metabric/

网址：https://ega-archive.org/dacs/EGAC00001000484

国际乳腺癌协会的分子分类数据库(Molecular Taxonomy of Breast Cancer International Consortium, METABRIC) 是一个加拿大-英国联合项目，旨在根据有助于确定最佳治疗过程的分子特征将乳腺肿瘤进一步分类。我们迄今为止已经根据肿瘤的基因指纹将乳腺癌重新分类为10个全新的类别。这些基因可以对乳腺癌生物学提供迫切需要的洞察力，使医生能够预测肿瘤是否会对某种特定的治疗产生反应。肿瘤是否有可能扩散到身体的其他部位，或者治疗后是否有可能复发。

最后BCL11A和MBNL1

来自维基百科

• BCL11A（BAF染色质重塑复合物亚基BCL11A）是一种蛋白质编码基因。与BCL11A相关的疾病包括智力发育障碍伴持续性胎儿血红蛋白和胼胝体发育不全。与该基因有关的基因本体论（GO）注释包括蛋白质同源二聚活性和RNA聚合酶II近端启动子序列特异性DNA结合。所述BCL11A用于调节C2H2型锌指蛋白基因编码，可结合到DNA。BCL11A在几个造血谱系中高度表达，并在胎儿向成人红细胞生成过渡过程中，从γ-珠蛋白表达转变为β-珠蛋白发挥作用。此外，BCL11A在大脑中表达，在大脑中与CASK形成蛋白质复合物，以调节轴突的生长和分支。在新皮质中，BCL11A与TBR1调节区结合并抑制TBR1的表达。相应BCL11A小鼠基因是一个常见的部位逆转录病毒整合在髓细胞性白血病。在造血细胞分化过程中，该基因被下调。它可能与淋巴瘤的发病机制有关，因为与B细胞恶性肿瘤相关的易位也使它的表达失控。另外，已经发现BCL11A在抑制胎儿血红蛋白产生中起作用。
• MBNL1 基因编码一种RNA剪接蛋白肌肉盲样剪接调节剂1（MBNL1）。肌盲蛋白与扩增的dsCUG RNA特异结合，但不与正常大小的CUG重复序列结合，因此可能在强直性肌营养不良的病理生理学中发挥作用。缺乏这种基因的老鼠表现出肌肉异常和白内障。已经描述了几个选择性剪接转录本变体，但仅确定了其中一些的全长性质。不同的异构体被认为具有不同的结合特异性和/或剪接活性。MBNL1(肌肉盲样剪接调节因子1)是一个蛋白质编码基因。与MBNL1相关的疾病包括强直性肌营养不良症和强直性肌强直病。它的相关途径包括脂肪形成。与该基因相关的基因本体论(GO)注释包括双链RNA结合。MBNL1基因介导前mRNA选择性剪接调控。在特定的Pre-mRNA靶点上作为剪接的激活物或抑制物。抑制心肌肌钙蛋白-T(TNNT2)前mRNA外显子包涵体，诱导胰岛素受体(IR)前mRNA外显子包涵体。拮抗CELF蛋白的选择性剪接活性模式。通过与U2AF2的竞争来调节TNNT2外显子5的跳跃。抑制TNNT2前mRNA内含子4上剪接体A复合体的形成。在剪接体组装过程中，与TNNT2内含子4的多嘧啶束内的茎环结构结合。与5‘-YGCu(U/G)Y-3’一致序列结合。

生存分析

• 随便去xena中含有DSS临床信息的数据中看看，发现其实DSS.time和OS.time时间是一样的，且不同结局的患者时间也是不变化的。即在没有DSS.time的数据中，可以用OS.time来替代（我猜测）

• 代码源于生活
• 以下代码均来自生信技能树Jimmy老师~
• 改动了一点数据

rm(list = ls())
library(cgdsr)
library(DT)

mycgds = CGDS("http://www.cbioportal.org/")
#test(mycgds)
setVerbose(mycgds, TRUE)
all_TCGA_studies = getCancerStudies(mycgds)
DT::datatable(all_TCGA_studies)
all_gen_pro = getGeneticProfiles(mycgds,'brca_metabric')
all_tables <- getCaseLists(mycgds,'brca_metabric')

# 获得两种基因的表达水平
my_dataset <- 'brca_metabric_mrna'
my_table <- "brca_metabric_mrna"
exp <- getProfileData(mycgds, c("BCL11A","MBNL1"), my_dataset, my_table)
exp$patient=rownames(exp)

# 获得临床信息
cli_dat = getClinicalData(mycgds,my_table)
myclinicaldata = cli_dat
DT::datatable(myclinicaldata,
              extensions = 'FixedColumns',
              options = list(                    #dom = 't',
                scrollX = TRUE,
                fixedColumns = TRUE
              ))

colnames(myclinicaldata)
cli_dat=myclinicaldata[,c("OS_MONTHS","VITAL_STATUS")]
cli_dat$patient=rownames(myclinicaldata)

# 融合两种数据并微调
tmp=merge(exp,cli_dat,by="patient")
tmp=na.omit(tmp)
colnames(tmp)[4]="DSS.time"
colnames(tmp)[5]="DSS"
tmp$DSS=ifelse(tmp$DSS=="Died of Disease",1,0)
tmp$DSS.time=tmp$DSS.time*30

library(survminer)
# 找到合适的“高低表达量”的定义，高于cut，就是高表达
tmp.cut <- surv_cutpoint(
  tmp,
  time = "DSS.time",
  event = "DSS",
  variables = c("BCL11A", "MBNL1")
)
summary(tmp.cut)
plot(tmp.cut, "BCL11A", palette = "npg")
# 将数字的表达信息根据粗体转化为high or low的二元信息
tmp.cat <- surv_categorize(tmp.cut)
head(tmp.cat)

library(survival)
# 常规的画法
fit <- survfit(Surv(DSS.time, DSS) ~ BCL11A + MBNL1,data = tmp.cat)
ggsurvplot(
  fit,
  pval = TRUE,
  conf.int = FALSE,
  title = 'METABRIC',
  xlim = c(0,10000),
  break.time.by = 2000,
  ggtheme = theme_RTCGA(),
  risk.table.y.text.col = T,
  risk.table.y.text = FALSE
)

• 不知道是按照什么的标准来分高低组
• 常见的分组方式可以包括：中位数、上/下四分位数，上/ 下三分位数，Z~score等
• 下面是中位数

group <- ifelse(tmp$BCL11A>median(tmp$BCL11A)&tmp$MBNL1<median(tmp$MBNL1),'Bh-ML',
                ifelse(tmp$BCL11A<median(tmp$BCL11A)&tmp$MBNL1>median(tmp$MBNL1),'BL-Mh',
                       ifelse(tmp$BCL11A>median(tmp$BCL11A)&tmp$MBNL1>median(tmp$MBNL1),'intermediate','intermediate')))

• 归一化后分组

tmp[,'BCL11A'] = scale(as.numeric(tmp[,'BCL11A']))
tmp[,'MBNL1'] = scale(as.numeric(tmp[,'MBNL1']))


group <- ifelse(tmp$BCL11A>0&tmp$MBNL1< 0,'Bh-ML',
                ifelse(tmp$BCL11A< 0&tmp$MBNL1> 0,'BL-Mh',
                       ifelse((tmp$BCL11A> 0&tmp$MBNL1>0)&(tmp$BCL11A< 0&tmp$MBNL1< 0),'intermediate','intermediate')))
 table(group)

这个分组有点接近原文了，Samples (all subtypes) were divided into BCL11A High + MBNL1 Low (n=148), BCL11A Low + MBNL1，High (n=345) or Intermediate (n=1478).那就这样吧

• 画出来的图跟原文超像有木有。

最后做出来的图就是上面这样了。需要找准分组的cutoff，关乎着你是否能重现文中的图表。

当然对于我而言，任务重要的是探索的过程，当然这句话Jimmy大大也经常说，教给你的不是现成的知识，而是你搜索记住的信息。授人以鱼不如授人以渔，棒！！！