生物信息学【3】：相关理论方法

由于最近在做BRCA-lncRNA相关的生信课题研究，在看到相关论文中的一些模型和方法，整理一下，供自己和大家一起学习。

理论方法：
- 1. 生存分析：log-rank检验在什么情况下失效?
- 2. DESeq2详细用法
- 3. TCGA+biomarker——风险因子关联图
- 4. ROC曲线，混淆矩阵，开集闭集等概念
- 5. 非负矩阵分解NMF
- 6. 绘制ceRNA网络结构
- 7. starBase 数据库的使用
- 8. STRING 蛋白质组研究
- 9. DNA 甲基化
- 10. WGCNA分析
- 11. 单因素Cox回归 -> 多因素Cox分析

理论方法：

1. 生存分析：log-rank检验在什么情况下失效?

（1）生存分析：log-rank检验在什么情况下失效?

这一篇博客，写了关于生存分析检验方法的差别。
介绍了：1) log-rank检验 (对数秩检验) 2）Wilcoxon检验 3）Peto检验、Tarone-Ware检验

作者的总结：
如果log-rank检验有意义而Wilcoxon检验无意义，表明可能远期差异较大，早期则不一定，有可能差异不大。
如果log-rank检验无意义而Wilcoxon检验有意义，表明早期生存差别较大，远期生存差异不大。

总的来说，如果研究某种疗法对生存常期是一种状态效果，并不会随着时间变换而衰弱，就可以运用log-rank检验（对数秩检验），如果是认为这种预期效果会随着时间变化而减弱，可以考虑到加入权重N的Wilcoxon检验。

2. DESeq2详细用法

（2）：DESeq2详细用法
这一篇博客，写了创建一个DESeq的对象dds，两种数据转换的方法：vst，rlog，【数据集小于30个样品可以用rlog，数据集大于30个样品用vst，因为rlog速度慢。】，这两种方法的目的：得到一个近似为同方差的值矩阵（沿均值范围具有恒定的方差）。**

3. TCGA+biomarker——风险因子关联图

（3）：TCGA+biomarker——风险因子关联图
这一篇博客是关于预后分析的关联图。在我们得到了风险评分后，如何通过图像来展现出高低风险组样本的分类情况，这时可以通过风险因子关联图来展现出高低风险的差异，其中风险因子关联图包括三个部分。
这里这部分，加一点我的代码：


# 由多因素COX得到的5个基因，我们要进行风险分析，
# 来重新定义BRCA，将他们分组。
rm(list = ls())
library(survival)
library(glmnet)
library(ggplot2)
library(survminer)setwd("D:\\AProject\\Study_code\\Code_026_Risk_model_lncRNA\\Code_026_Risk_vst_v2\\Risk")# 读取表达数据
lncRNA_data <- read.csv("tcga_lncRNA_significance_in_immunecell_pathway.csv", row.names = 1)
lncRNA_exp <- lncRNA_data[, 1:7]#####(1)建立多因素cox回归的数据#######
fml <- as.formula(Surv(lncRNA_data$OS.time, lncRNA_data$status)~.)
mycox <- coxph(fml, data = lncRNA_exp)
summary(mycox)# 生存风险得分risk_level和评分risk_score
risk_score <- predict(mycox, type="risk", newdata = lncRNA_exp)
risk_level <- as.data.frame(ifelse(risk_score > median(risk_score), "High", "Low"))
colnames(risk_level) <- "risk_level"
risk_score <- as.data.frame(risk_score)
colnames(risk_score) <- "risk_score"
# 生存最后的数据
dat <- cbind(lncRNA_data, risk_score, risk_level)
write.csv(dat, "risk_score.csv")######(2)生存分析#####
colnames(dat)
# 影响因素分析
# survival包中的Surv函数可以创建一个生存对象
# gender: 0(wumen)、1(men)
fit <- survfit(Surv(OS.time1, status1) ~ risk_level, data = dat)
#survival包中的survfit函数用Kaplan-Meier法进行生存曲线的拟合
sur <- ggsurvplot(fit,pval = TRUE, # conf.int = TRUE,risk.table = TRUE, # Add risk tablerisk.table.col = "strata", # Change risk table color by groupslinetype = "strata", # Change line type by groupssurv.median.line = "hv", # Specify median survivalggtheme = theme_bw(), # Change ggplot2 themepalette = c("#E7B800", "#2E9FDF")
)# 存储图片
pdf("survival_sig_lncRNA_level1.pdf")#生成文件
sur
dev.off()#####(3) 绘制风险因子关联图#####
phe <- dat[order(dat$risk_score), ]
fp_dat <- data.frame(patientid = 1:nrow(phe),fp = phe$risk_score)
#添加风险分组，以风险评分的中位值将患者分为两组，大于中位值的 患者为高风险组，小于或等于中位值的患者为低风险组
fp_dat$riskgroup <- ifelse(fp_dat$fp>= median(fp_dat$fp),'high','low')###第一个图
library(ggplot2)
p1 = ggplot(fp_dat,aes(x = patientid,y = fp))+ geom_point(aes(color = riskgroup))+scale_colour_manual(values = c("red","green"))+theme_bw()+labs(x="Patient ID(increasing risk score)",y="Risk score")+geom_hline(yintercept=median(fp_dat$fp),colour="black", linetype="dotted",size=0.8)+geom_vline(xintercept=sum(fp_dat$riskgroup=="low"),colour="black", linetype="dotted",size=0.8)
p1# 生成sur_dat绘制图二
sur_dat <- data.frame(patientid = 1:nrow(risk_score), time = phe[,'OS.time1'],event = phe[, 'status1'])
sur_dat$event <- ifelse(sur_dat$event==0, 'alive', 'death')
sur_dat$event <- factor(sur_dat$event, levels = c("death","alive"))#### 第二个图
p2 <- ggplot(sur_dat, aes(x=patientid,y=time)) + geom_point(aes(col=event)) + theme_bw()+scale_colour_manual(values = c("red","green"))+labs(x = "Patient ID(increasing risk score)", y = "Survival time(year)")+geom_vline(xintercept=sum(fp_dat$riskgroup=="low"), colour="black", linetype="dotted",size=0.8)
p2#### 第三个图
library(pheatmap)
heatmap.data <- phe[, 1:7]  # 基因数据框heatmap.data[is.na(heatmap.data)] <- 0
heatmap.data <- as.matrix(heatmap.data)
heatmap.data.scale <- scale(heatmap.data)
heatmap.data.scale.new <- ifelse(abs(heatmap.data.scale) > 1, sign(heatmap.data.scale)*1, heatmap.data.scale)# 绘制热图
library(ComplexHeatmap)
library(qdapTools)  #need df2matrix function
library(RColorBrewer)  #colorRamp2
library(circlize)  #colorRamp2
library(tidyverse)  #select# 转置，得到行为基因，列为样本的矩阵
mat = t(heatmap.data.scale.new)
dim(mat)# 在热图上添加组标签
Groups <- as.character(phe$risk_level)
# 添加标签
annotation_col = HeatmapAnnotation(Groups = Groups,col = list(Groups = c("Low" = "#556B2F", "High" = "royalblue")),annotation_name_side = "right",  # 设置注释的名字在右边na_col = "#808080",  #设置空白值的颜色为灰色simple_anno_size = unit(5, "mm")  # 设置行宽度
)p3 = Heatmap(mat, col = colorRamp2(c(-1, 0, 1), c("CornflowerBlue", "#D3D3D3", "FireBrick")),column_title_gp = gpar(fontsize = 8, fontface = "bold"),top_annotation = annotation_col,# column_split = 2,# row_split = 2,# row_labels = FALSE,# row_names_gp = gpar(fontsize = 20),  # 设置行字体的大小# column_title = "Kmeans groups with heatmap",  # 标签名称cluster_columns = FALSE,  # 列不进行聚类# cluster_columns = hclust(dist(t(mat1))),  #列是样本，样本需要聚类# clustering_distance_columns = "euclidean",  #列是样本，样本需要聚类# clustering_method_columns = "complete",cluster_rows = hclust(dist(mat)),clustering_distance_rows = "euclidean",  #行是基因，基因需要聚类clustering_method_rows = "complete",show_row_dend = FALSE,  #是否显示树状图show_column_dend = FALSE,show_row_names = T,show_column_names = F,# column_names_gp = gpar(fontsize = 1),# row_names_gp = gpar(fontsize = 1),heatmap_legend_param = list(title = "Legend",title_position ="topcenter",title_gp = gpar(fontsize = 10, fontface = "bold"),labels_gp = gpar(fontsize = 10))
)
p3### 拼图实现三图联动
library(ggplotify)
plots = list(A = p1,B = p2,C = as.ggplot(as.grob(p3)))
library(gridExtra)
lay1 = rbind(c(rep(1,7)),c(rep(2,7)),c(rep(3,7))) #布局矩阵
riskdistru <- grid.arrange(grobs = plots, layout_matrix = lay1, heigths = c(2,2,3),weights=c(10,10,10))
ggsave("Risk factor association diagram.pdf", plot = riskdistru, dpi=600, width = 6.73, height = 6.7)# ggsave("风险因子关联图.tiff", plot = riskdistru, dpi=600, width = 6.73,height = 6.7)####(4) ROC曲线图#####
library(timeROC)
library(survival)
ROC <- timeROC(T= dat$OS.time/12, delta = dat$status,marker=dat$risk_score, cause=1,weighting="marginal",times=c(3,5,2),ROC=TRUE)# 显示全部细胞和基质占比
pdf(file = "ROC.pdf",width=5,height=5, useDingbats = FALSE)
plot(ROC,time=3,col="blue", title=FALSE, lwd=3)
plot(ROC,time=5,col="red", add=TRUE, title=FALSE,lwd=3)
legend("bottomright",title = 'AUC',c(paste("3-year: ",round(ROC$AUC[1],1)),paste("5-year: ",round(ROC$AUC[2],1))),col=c("blue","red"),lwd=2)dev.off()# ggsave("TCGA生存+timeROC.pdf", plot = ROC,dpi=600, width = 10, height = 4)

4. ROC曲线，混淆矩阵，开集闭集等概念

（4）：ROC、AUC、混淆矩阵

roc_curve()，ROC曲线，混淆矩阵，开集闭集等概念
混淆矩阵、召回率、准确率、ROC曲线、AUC
分类模型评估之ROC-AUC曲线和PRC曲线
混淆矩阵、ROC、AUC

特别是看一下第四个博客，博客中有详细讲解了ROC曲线中点是如何绘制出来的。
可以看到，ROC中的点是从原点（0，0）移动到（1，1）,依据真正类率(true positive rate ,TPR)和假正类率(true positive rate ,FPR)的值作为坐标值，同时，在此之前是按照得分排名的，个人认为：每一个Score下的矩阵表示形式为一个：混淆矩阵，通过计算TPR，FPR来最后得到该概率下的点坐标。

最理想的状况下是先上移动到（0，1），再右移动到（1，1）；理解就是，我们将样本分成的高低风险组与原有的分组之间是完全匹配的。
　　　　　　　　　　　

图片来自博主：chuanbanjun

5. 非负矩阵分解NMF

文章链接: 代谢相关的大肠癌分子分类.
链接: 非负矩阵分解NMF.

可以通过非负矩阵的分解将原始V矩阵（n×m）分解成W（n×r）与矩阵H（r×m）的乘积，矩阵H也叫系数矩阵，通常系数矩阵能够代替原始矩阵，实现降维。【其中矩阵V的每一个列向量为矩阵W所有列向量的线性组合的值，系数为矩阵H 的某一列】

6. 绘制ceRNA网络结构

如何使用Cytoscape绘制出ceRNA网络图

其他学习链接：
(1)ceRNA预测工具 —— GDCRNATools
http://www.360doc.com/content/18/0119/10/45852776_723289997.shtml
(2) Cytoscape史上最全攻略
http://www.360doc.com/content/19/0409/20/49059453_827533578.shtml

7. starBase 数据库的使用

官网：http://starbase.sysu.edu.cn/
（1）https://www.bilibili.com/video/BV1qh411o7UG?from=search&seid=2485808999820599866

我们在下载mRNA-Target靶基因数据时，有以下的参数需要知道：

colnames(lnc_mi)##  [1] "miRNAid"      "miRNAname"    "geneID"       "geneName"     "geneType"
##  [6] "chromosome"   "start"        "end"          "strand"       "clipExpNum"
## [11] "degraExpNum"  "RBP"          "merClass"     "miRseq"       "align"
## [16] "targetSeq"    "pancancerNum"

对starbase数据库提供的数据列名的解释里，有两个比较重要的：

clipExpNum ：The number of CLIP-seq experiments ;
表示CLIP测序实验支持的个数，而数据库默认的时低严格度，有一个测序结果就可以。

pancancerNum : Number of Cancer types (Pan-Cancer) that miRNA-target show anti-correlation relationships (pearson correlation: r<0, p-value<0.05).
所以pancerNum是几：就意味着这对miRNA-lncRNA在多少种癌症中表达量负相关，省掉了很多计算。

program number : 表示数据库预测的个数，如果选择三个表示至少保持三个数据库的预测结果。

数据库的tutorial里面还提到了一个比较严格的筛选标准：
CLIP evidence (>=5), degradome evidence (>=1), Pan-Cancer (>=10),
program number (>=5) and predicted program (None).

8. STRING 蛋白质组研究

官网：https://string-db.org/
（1）https://www.biomart.cn/experiment/793/2713465.html
（2）https://www.omicsclass.com/article/1126
（3）http://www.360doc.com/content/19/0419/13/52645714_829884871.shtml

9. DNA 甲基化

什么是 DNA 甲基化?

DNA甲基化是在DNA甲基化转移酶(Dnmt)的作用下将甲基选择性地添加到胞嘧啶上形成5-胞嘧啶的过程，刚被发现时被定义为第五种碱基，实际上它是一种重要的表观遗传学标记，在调控基因表达、维持染色质结构、基因印记、X染色体失活以及胚胎发育等生物学过程中发挥着重大的作用。

甲基化与肿瘤:
这些因素是导致肿瘤发展、转移、恶化最终导致患者死亡的重要原因。
(1) 甲基化的CpG岛二核苷酸中的胞嘧啶以较高的频率脱氨基变成胸腺嘧啶，造成基因突变；
(2) 抑癌基因和DNA修复基因由于超甲基化而沉默；
(3) 癌基因甲基化水平降低而活化；
(4) 基因组总体甲基化水平降低使转座子、重复序列活化导致染色体稳定性下降。

甲基化与肿瘤侵袭和转移:
利用DNA甲基转移酶（methyltransferase）抑制剂5-氮杂胞苷（5-aza-cytidine）处理MCF-7乳腺癌细胞，使其维持低甲基化状态，结果显示与EMT过程相关的促细胞侵袭基因（pro-invasive EMT-associated gene）表达上调，细胞的侵袭能力和转移能力增强。

甲基化与肿瘤治疗:
肿瘤预防和治疗的一个手段是通过去甲基化恢复某些关键的抑癌基因或DNA修复基因的活性，目前研究最多的是DNMTs抑制剂，它通过抑制DNMT活性以逆转异常的DNA甲基化。

10. WGCNA分析

（1）基本概念：
WGCNA其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。

（2）基本原理：
从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。

第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。因此，WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂，使得网络中的基因之间的连接服从无尺度网络分布(scale-freenetworks)，这种算法更具生物学意义。

第二步通过基因之间的相关系数构建分层聚类树，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。

（3）基本流程：

构建基因共表达网络：使用加权的表达相关性。
识别基因集：基于加权相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。
如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。
研究模型之间的关系，从系统层面查看不同模型的互作网络。
从关键模型中选择感兴趣的驱动基因，或根据模型中已知基因的功能推测未知基因的功能。
导出TOM矩阵，绘制相关性图。

具体参考下面两个博客：
（1）http://blog.sciencenet.cn/blog-118204-1111379.html
（2）http://www.bio-info-trainee.com/2535.html

11. 单因素Cox回归 -> 多因素Cox分析

1. 为什么使用单因素到多因素Cox回归分析？

在我看来是在考虑生存状态的情况下，判断一系列临床特征是否具有显著意义。当然这里的显著意义也就是说"回归"。因为，我们在使用多因素Cox分析时会得到显著性样本的系数，可以当作成对应的回归函数的系数。如下图系数：

通过上面的系数，我们可以建立函数：

当然这里输入的不是临床特征，而是差异基因lncRNAs，同样我们也是可以当作成一种特征，目的是分析一个整体（这些差异基因的联合），得到系数coef后，最后评出风险评分Risk score。

代码部分，具体可以参考这篇博客：
单因素Cox回归 -> 多因素Cox分析

2. 单因素和多因素结果分析

一般情况下，单因素和多因素分析相辅相成，单因素分析可以初步探索预测变量与响应变量的关系，并且当样本量不是很大的时候可以通过单因素分析删除部分无关的预测变量；而多因素分析可以进一步排出其它混杂因素的影响，从而确定预测变量与响应变量的相关性；此外，将单因素和多因素的结果进行比较更容易发现问题，如果单因素和多因素结果一致，结论就比较稳定且容易解释，但在一些情况下，结果可能存在差异。举个最简单的例子：如果其它变量（如a、b、c）与研究变量（如lx）之间是负向关系，而研究变量（如lx）与因变量（如y）之间是正向关系，且其它变量（如a、b、c）与因变量（如y）之间也是正向关系。这时候就很容易出现“单因素分析没有统计学意义，而多因素分析中有统计学意义”。

单因素和多因素分析的系统比较（转载自医咖会）

情况A： 单因素分析和多因素分析的结果都显示无统计学显著性，两者结果一致，均为阴性结果

在这种情况下，结果还是相对比较好解释的，一般基本上可以认为该因素对于结局事件来说，不是一个有意义的影响因素。

但是事情也并非这么简单，如果该因素作为一个混杂因素，在多因素分析中只是用来起到调整混杂作用的目的，那么虽然它在单因素和多因素分析中都是阴性结果，可能也不会太引起研究人员的重视；但是如果该因素是研究中所重点关注的一个因素，例如暴露/处理因素，此时单因素和多因素分析都得出阴性结果的话，就会让人感觉比较沮丧，不过也更值得我们好好去思考一下阴性结果背后的意义。

到底是该暴露/处理因素对结局事件真的没有影响，还是说因为其他原因而导致它的实际效应没有被显示出来？到底是研究设计的问题，还是指标定义的问题，亦或是统计方法的问题呢？都需要我们认真去查找一下原因，可以参考前期推送的有关介绍“阴性结果”的系列文章，或许可以帮助你寻找一下产生阴性结果的原因，开拓一下分析思路。

情况B： 单因素分析结果显示无统计学显著性，但多因素分析结果显示有统计学显著性

这种情况可能并不常见，主要是因为在单因素分析中没有统计学显著性的因素，按照一般的做法就不会再将此变量纳入到多因素分析中了，但其实上述做法小咖并不推荐，它是存在一定缺陷的。

我们在前期介绍《如何理解回归模型中的“调整”和“独立作用”》的内容中讲到，在单因素分析中，由于自变量之间存在一定的相互关联，自变量对因变量的影响反映的不仅仅单纯是它本身的作用，而是包含了该变量自身作用以及其他变量的混杂作用之后，呈现出来的一个综合的结果。而在多因素分析中，通过构建回归模型，调整了其他混杂因素的影响，从而才使该因素对因变量的真实效应显示出来。

因此不难理解，当某因素在单因素分析结果中无统计学显著性，而多因素分析结果有统计学显著性时，此时可能的原因是，该因素与其他混杂因素之间可能存在一定的关联，在单因素分析时，该因素的真实效应被其他混杂因素的作用所掩盖，通过多因素分析消除其他因素的影响后，才发现原来该因素对于结局事件来说是具有独立作用的。

举一个例子，例如某因素A是一个危险因素，而因素B是一个保护因素，由于具有因素A的个体，大部分人同时也具有因素B，因此在单因素分析中，因素A的作用并没有显现出来，这是因为因素A的危险作用被因素B的保护作用所掩盖了，无法体现因素A的实际效应。而通过多因素分析，将因素B的保护作用进行调整，从而暴露出因素A真实的危险作用。

情况C： 单因素分析结果显示有统计学显著性，但多因素分析结果显示无统计学显著性

想必大家都会经常遇见到这种情况，单因素分析时该因素有统计学显著性，然后就很兴奋地把它扔进多因素分析中，结果多因素分析结果却显示没有统计学显著性，感觉前功尽弃，很让人头痛，不知道是出了什么问题，到底该怎么办了。

如果你对情况B产生的原因已经理解，那么情况C也是同样的道理。在单因素分析中，自变量与因变量之间可能出现一定的假关联或者是间接的关联，例如某因素A对结局事件并无影响，而因素B对于结局事件是一个影响因素，但是由于因素A只是单纯的和因素B有强烈的相关性，两者存在共线性的现象，那么在单因素分析中，就可能出现因素A也存在显著差异的结果，从而导致因素A被误认为是一个影响因素而纳入到多因素分析中。

而在多因素分析中通过调整因素B的影响，因素A与因变量的“假关联”就消失了，此时可以认为因素A实际上对于结局事件并非是一个影响因素。

情况D： 单因素分析和多因素分析的结果都显示统计学显著性，两者结果一致，均为阳性结果

这种情况应该是大家最愿意看到的情况吧，往往单因素和多因素分析都出现阳性结果，以此结果作出的结论还算是比较稳定可靠，可以放心地写文章投稿了，但前提是单因素和多因素分析的阳性结果的方向是一致的，比如单因素分析显示病例组某因素的水平显著高于对照组，多因素分析也显示该因素为危险因素，两者的结果都倾向于该因素对结局事件具有危险作用。

不过偶尔也会遇见这样的情况，虽然单因素和多因素分析都得出阳性结果，但是有时单因素分析显示为危险因素，而多因素分析显示为保护因素，或者单因素分析显示为保护因素，而多因素分析显示为危险因素，两者的结果是相互矛盾的。

出现这样的情况，其实和上述的情况B和C是同样的道理，这是在统计分析中经常出现的一个陷阱，统计学上称之为“辛普森悖论”(Simpson’s Paradox)，是由英国统计学家E.H.Simpson于1951年提出。简单理解就是，在某个条件下的两组数据，分别讨论时都会满足某种性质，可是一旦将两组数据合并考虑，却可能导致相反的结论。

自变量筛选常用方法

参考：单因素和多因素分析