技術標籤：生物資訊學R語言生物資訊學統計學

由於最近在做BRCA-lncRNA相關的生信課題研究，在看到相關論文中的一些模型和方法，整理一下，供自己和大家一起學習。

目錄

• • 理論方法：
  • • 1. 生存分析：log-rank檢驗在什麼情況下失效?
    • 2. DESeq2詳細用法
    • 3. TCGA+biomarker——風險因子關聯圖
    • 4. ROC曲線，混淆矩陣，開集閉集等概念

理論方法：

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 table
                  risk.table.col = "strata", # Change risk table color by groups
                  linetype = "strata", # Change line type by groups
                  surv.median.line = "hv", # Specify median survival
                  ggtheme = theme_bw(), # Change ggplot2 theme
                  palette = 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、混淆矩陣

1. roc_curve()，ROC曲線，混淆矩陣，開集閉集等概念
2. 混淆矩陣、召回率、準確率、ROC曲線、AUC
3. 分類模型評估之ROC-AUC曲線和PRC曲線
4. 混淆矩陣、ROC、AUC

特別是看一下第四個部落格，部落格中有詳細講解了ROC曲線中點是如何繪製出來的。
可以看到，ROC中的點是從原點（0，0）移動到（1，1）,依據真正類率(true positive rate ,TPR)和假正類率(true positive rate ,FPR)的值作為座標值，同時，在此之前是按照得分排名的，個人認為：每一個Score下的矩陣表示形式為一個：混淆矩陣，通過計算TPR，FPR來最後得到該概率下的點座標。

最理想的狀況下是先上移動到（0，1），再右移動到（1，1）；理解就是，我們將樣本分成的高低風險組與原有的分組之間是完全匹配的。
　　　　　　　　　　　

圖片來自博主：chuanbanjun