詳細(xì)情況請(qǐng)前往
數(shù)據(jù)分析：完整的成體系的生存分析的結(jié)果匯總
介紹
本章節(jié)全面匯總了生存分析的相關(guān)各類分析方法攀芯，形成了一個(gè)系統(tǒng)化的生存分析教程两入。通過(guò)這個(gè)教程，讀者可以深入了解生存分析所涵蓋的多種分析技術(shù)敲才。組合下列結(jié)果的可視化可見(jiàn): 科研繪圖系列：組合多個(gè)文章圖
數(shù)據(jù)下載
本教程所需要的所有數(shù)據(jù)下載鏈接見(jiàn)：
- 百度網(wǎng)盤鏈接: https://pan.baidu.com/s/1KZLVvyrU1Yj78KL0CYYygQ
- 提取碼:
image.png
多變量生存分析
多變量生存分析模塊是一個(gè)統(tǒng)計(jì)分析工具裹纳，它提供了兩個(gè)主要步驟來(lái)識(shí)別和評(píng)估影響生存時(shí)間的變量：
- 單變量生存分析篩選：首先，該模塊基于單變量生存分析來(lái)篩選與生存時(shí)間相關(guān)的特征紧武。這一步驟涉及對(duì)每個(gè)變量單獨(dú)進(jìn)行生存分析剃氧，以確定它們與生存時(shí)間的關(guān)聯(lián)程度。通過(guò)這種方式阻星，可以識(shí)別出在單變量分析中顯著相關(guān)的變量朋鞍，為后續(xù)的多變量分析打下基礎(chǔ)已添。
- 多變量生存分析：其次，將單變量分析中顯著相關(guān)的特征組合在一起滥酥，進(jìn)行多變量生存分析更舞。這一步驟旨在評(píng)估多個(gè)變量同時(shí)對(duì)生存時(shí)間的影響，以及它們之間的相互作用坎吻。多變量生存分析可以提供更全面的風(fēng)險(xiǎn)評(píng)估缆蝉，因?yàn)樗紤]了多個(gè)變量的聯(lián)合效應(yīng)，從而允許更精確地預(yù)測(cè)生存時(shí)間和識(shí)別重要的預(yù)后因素瘦真。
library(tidyverse)
library(survminer)
library(survival)
library(reshape2)
library(cowplot)
library(caret)
library(pROC)
# Load data
############
load("./data/MHC_immunotherapy.RData")
ICB_data <- ICB_study[ICB_study$study == "Liu_2019", ]
# Remove NA
############
ICB_data <- na.omit(ICB_data)
# Binarize
#############
ICB_data$highTMB <- ICB_data$TMB > quantile(ICB_data$TMB, 0.5, na.rm = T)
ICB_data$strongMHC1 <- ICB_data$MGBS1 > quantile(ICB_data$MGBS1, 0.5, na.rm = T)
ICB_data$strongMHC2 <- ICB_data$MGBS2 > quantile(ICB_data$MGBS2, 0.5, na.rm = T)
ICB_data$strongMHCnorm <- ICB_data$MGBSnorm > quantile(ICB_data$MGBSnorm, 0.5, na.rm = T)
ICB_data$highPD1 <- ICB_data$PD1_noBatch > quantile(ICB_data$PD1_noBatch, 0.5, na.rm = T)
ICB_data$highCYT <- ICB_data$CYT_noBatch > quantile(ICB_data$CYT_noBatch, 0.5, na.rm = T)
ICB_data$highMHC1 <- ICB_data$MHC1_noBatch > quantile(ICB_data$MHC1_noBatch, 0.5, na.rm = T)
ICB_data$highMHC2 <- ICB_data$MHC2_noBatch > quantile(ICB_data$MHC2_noBatch, 0.5, na.rm = T)
ICB_data$highPurity <- ICB_data$purity > quantile(ICB_data$purity, 0.5, na.rm = T)
ICB_data$highPloidy <- ICB_data$ploidy > quantile(ICB_data$ploidy, 0.5, na.rm = T)
ICB_data$highHeterogeneity <- ICB_data$heterogeneity > quantile(ICB_data$heterogeneity, 0.5, na.rm = T)
ICB_data$highLDH <- ICB_data$LDH == 1
# Save
######
dir.create("./results/data/", recursive = T)
save.image(file = "./results/data/manuscript_multivariate.RData")
saveRDS(list(Cox_model_ls = Cox_model_ls, LR_model_ls = LR_model_ls), file = "./results/data/manuscript_multivariate.rds")
相關(guān)分析
library(corrplot)
library(RColorBrewer)
library(cowplot)
load("./data/MHC_immunotherapy.RData")
ICB_data <- ICB_study[ICB_study$study == "Liu_2019", ]
head(ICB_data)
# Create correlation plot
###########################
vars <- c("TMB", "neoAgB1", "neoAgB2", "neoAgBnorm", "pMHC1", "pMHC2",
"pMHCnorm", "MGBS1", "MGBS2", "MGBSnorm", "MHC1_noBatch",
"MHC2_noBatch", "heterogeneity", "ploidy", "purity",
"MHC1_zyg", "MHC2_zyg", "PD1_noBatch", "B2M_noBatch", "CYT_noBatch")
cor_df <- ICB_data[, vars]
colnames(cor_df) <- c("TMB", "abs. MHC-I neoAgB", "abs. MHC-II neoAgB",
"abs. diff. neoAgB", "norm. MHC-I neoAgB", "norm. MHC-II neoAgB",
"norm. diff. neoAgB", "MGBS-I", "MGBS-II", "MGBS-d", "MHC-I expr.",
"MHC-II expr.", "heterogeneity", "ploidy", "purity", "MHC-I zygosity",
"MHC-II zygosity", "PD-L1 (CD274) expr.", "B2M expr.", "CYT")
# Save
#########
save.image(file = "./results/data/manuscript_correlation_analysis.RData")
saveRDS(p, file = "./results/data/manuscript_correlation_analysis.rds")
比例分析
library(survminer)
library(survival)
library(reshape2)
# Load data
############
load("./data/MHC_immunotherapy.RData")
ICB_data <- ICB_study[ICB_study$study == "Liu_2019", ]
# MHC binders
#############
ICB_data$highTMB <- ICB_data$TMB > quantile(ICB_data$TMB, 0.5, na.rm = T)
ICB_data$strongMHC1 <- ICB_data$MGBS1 > quantile(ICB_data$MGBS1, 0.5, na.rm = T)
ICB_data$strongMHC2 <- ICB_data$MGBS2 > quantile(ICB_data$MGBS2, 0.5, na.rm = T)
ICB_data$strongMHCnorm <- ICB_data$MGBSnorm > quantile(ICB_data$MGBSnorm, 0.5, na.rm = T)
# Save
######
save.image(file = "./results/data/manuscript_RECIST.RData")
saveRDS(p_BR_ls, file = "./results/data/manuscript_RECIST.rds")
生存分析
生存分析是一種統(tǒng)計(jì)方法刊头，用于描述、測(cè)量和分析事件的特征诸尽，尋找其發(fā)生的原因原杂，并對(duì)生存以及事件發(fā)生的時(shí)間進(jìn)行預(yù)測(cè)。生存分析的目的在于研究哪些因素影響了事件的發(fā)生速度及生存時(shí)間的長(zhǎng)短您机，它探索的因變量是一個(gè)包含了刪失或者截除數(shù)據(jù)的事件時(shí)間變量穿肄。生存分析涉及的事件可以是患者的存活時(shí)間、疾病的復(fù)發(fā)际看、企業(yè)的存活時(shí)間等咸产。
生存分析的目的包括但不限于以下幾點(diǎn)：
- 描述生存時(shí)間：描述生存時(shí)間的分布情況，包括生存時(shí)間的中位數(shù)仿村、平均值等統(tǒng)計(jì)量锐朴。
- 識(shí)別風(fēng)險(xiǎn)因素：識(shí)別和評(píng)估影響生存時(shí)間的各種風(fēng)險(xiǎn)因素兴喂。
- 預(yù)測(cè)生存概率：預(yù)測(cè)個(gè)體在特定時(shí)間點(diǎn)之前生存或事件發(fā)生的概率蔼囊。
- 比較生存情況：比較不同群體或不同條件下的生存情況，例如比較不同治療方案的效果衣迷。
- 評(píng)估干預(yù)效果：評(píng)估醫(yī)療干預(yù)畏鼓、政策或其他措施對(duì)生存時(shí)間的影響。
生存分析統(tǒng)計(jì)方法涉及統(tǒng)計(jì)描述壶谒、差異性分析和回歸分析三大類云矫。其中，統(tǒng)計(jì)描述主要有Kaplan-Meier估計(jì)法和Life table估計(jì)法汗菜；差異性分析主要有對(duì)數(shù)秩檢驗(yàn)和Wilcoxon test让禀；而回歸分析主要有Cox比例和非比例風(fēng)險(xiǎn)回歸模型、參數(shù)回歸模型陨界。
library(survminer)
library(survival)
library(reshape2)
library(tidyverse)
# Load data
############
load("./data/MHC_immunotherapy.RData")
# Save
#########
save.image(file = "./results/data/manuscript_validation.RData")
saveRDS(list(p_val_ls = p_val_ls, HR_df_all = HR_df_all), file = "./results/data/manuscript_validation.rds")
生存分析2
生存分析是一種統(tǒng)計(jì)方法巡揍，用于描述、測(cè)量和分析事件的特征菌瘪，尋找其發(fā)生的原因腮敌，并對(duì)生存以及事件發(fā)生的時(shí)間進(jìn)行預(yù)測(cè)。生存分析的目的在于研究哪些因素影響了事件的發(fā)生速度及生存時(shí)間的長(zhǎng)短，它探索的因變量是一個(gè)包含了刪失或者截除數(shù)據(jù)的事件時(shí)間變量糜工。生存分析涉及的事件可以是患者的存活時(shí)間弊添、疾病的復(fù)發(fā)、企業(yè)的存活時(shí)間等捌木。
生存分析的目的包括但不限于以下幾點(diǎn)：
- 描述生存時(shí)間：描述生存時(shí)間的分布情況油坝，包括生存時(shí)間的中位數(shù)、平均值等統(tǒng)計(jì)量钮莲。
- 識(shí)別風(fēng)險(xiǎn)因素：識(shí)別和評(píng)估影響生存時(shí)間的各種風(fēng)險(xiǎn)因素免钻。
- 預(yù)測(cè)生存概率：預(yù)測(cè)個(gè)體在特定時(shí)間點(diǎn)之前生存或事件發(fā)生的概率。
- 比較生存情況：比較不同群體或不同條件下的生存情況崔拥，例如比較不同治療方案的效果极舔。
- 評(píng)估干預(yù)效果：評(píng)估醫(yī)療干預(yù)、政策或其他措施對(duì)生存時(shí)間的影響链瓦。
生存分析統(tǒng)計(jì)方法涉及統(tǒng)計(jì)描述拆魏、差異性分析和回歸分析三大類。其中慈俯，統(tǒng)計(jì)描述主要有Kaplan-Meier估計(jì)法和Life table估計(jì)法渤刃；差異性分析主要有對(duì)數(shù)秩檢驗(yàn)和Wilcoxon test；而回歸分析主要有Cox比例和非比例風(fēng)險(xiǎn)回歸模型贴膘、參數(shù)回歸模型卖子。
library(tidyverse)
library(scales)
library(survival)
# Load expression data and metadata
## Load ICB metadata in an encapsulated object "all_data"
all_data <- new.env()
load("./data/MHC_immunotherapy.RData", all_data)
metadata <- all_data$ICB_study
# Save
#########
save.image(file = "./results/data/manuscript_validation_compare_variables.RData")
saveRDS(p_ls, file = "./results/data/manuscript_validation_compare_variables.rds")
logreg分析 & ROC分析
logreg分析和ROC分析是兩種不同的統(tǒng)計(jì)方法，它們?cè)跈C(jī)器學(xué)習(xí)和數(shù)據(jù)分析中有著各自的應(yīng)用和目的刑峡。
logreg分析（Logistic Regression Analysis）
logreg分析洋闽，即邏輯回歸分析，是一種用于二分類問(wèn)題的統(tǒng)計(jì)方法突梦。它通過(guò)估計(jì)特定變量（或變量組合）對(duì)結(jié)果發(fā)生概率的影響來(lái)工作诫舅。邏輯回歸使用最大似然估計(jì)（MLE）方法來(lái)估計(jì)模型參數(shù)，即找到能夠使得觀測(cè)數(shù)據(jù)出現(xiàn)概率（似然函數(shù)）最大的參數(shù)值宫患。
ROC分析（Receiver Operating Characteristic Analysis）
ROC分析是一種用于評(píng)估分類模型性能的方法刊懈，特別適用于比較不同分類器的性能。ROC分析通過(guò)繪制真正例率（True Positive Rate, TPR）和假正例率（False Positive Rate, FPR）的對(duì)比圖來(lái)展示分類器的性能娃闲。TPR是正確分類為正例的比例虚汛，而FPR是錯(cuò)誤分類為正例的比例。ROC曲線下的面積（AUC）是衡量分類器性能的一個(gè)重要指標(biāo)皇帮，AUC值越高卷哩，表明分類器的區(qū)分能力越強(qiáng)
library(tidyverse)
library(survminer)
library(survival)
library(reshape2)
library(cowplot)
library(caret)
library(pROC)
# Load data
load("./data/MHC_immunotherapy.RData")
p_ls <- list()
# Get model from Liu - preIpi
ICB_data <- ICB_study[!is.na(ICB_study$study) & ICB_study$study == "Liu_2019", ]
ICB_data <- ICB_data[ICB_data$preIpi == T, ]
ICB_data$MGBS2_z <- (mean(ICB_data$MGBS2, na.rm = T) - ICB_data$MGBS2) / sd(ICB_data$MGBS2, na.rm = T) # Z-score normalization
ICB_data$MHC2_z <- (mean(ICB_data$MHC2, na.rm = T) - ICB_data$MHC2) / sd(ICB_data$MHC2, na.rm = T) # Z-score normalization
model <- glm(R ~ MHC2_z + MGBS2_z, data = ICB_data, family = binomial)
p_ls$model <- model
# summary(model)
# Save
#########
save.image(file = "./results/data/manuscript_validation_multivariate.RData")
saveRDS(p_ls, file = "./results/data/manuscript_validation_multivariate.rds")
生存分析3
生存分析是一種統(tǒng)計(jì)方法，用于描述玲献、測(cè)量和分析事件的特征殉疼，尋找其發(fā)生的原因梯浪，并對(duì)生存以及事件發(fā)生的時(shí)間進(jìn)行預(yù)測(cè)。生存分析的目的在于研究哪些因素影響了事件的發(fā)生速度及生存時(shí)間的長(zhǎng)短瓢娜，它探索的因變量是一個(gè)包含了刪失或者截除數(shù)據(jù)的事件時(shí)間變量挂洛。生存分析涉及的事件可以是患者的存活時(shí)間、疾病的復(fù)發(fā)眠砾、企業(yè)的存活時(shí)間等虏劲。
生存分析的目的包括但不限于以下幾點(diǎn)：
- 描述生存時(shí)間：描述生存時(shí)間的分布情況，包括生存時(shí)間的中位數(shù)褒颈、平均值等統(tǒng)計(jì)量柒巫。
- 識(shí)別風(fēng)險(xiǎn)因素：識(shí)別和評(píng)估影響生存時(shí)間的各種風(fēng)險(xiǎn)因素。
- 預(yù)測(cè)生存概率：預(yù)測(cè)個(gè)體在特定時(shí)間點(diǎn)之前生存或事件發(fā)生的概率谷丸。
- 比較生存情況：比較不同群體或不同條件下的生存情況堡掏，例如比較不同治療方案的效果。
- 評(píng)估干預(yù)效果：評(píng)估醫(yī)療干預(yù)刨疼、政策或其他措施對(duì)生存時(shí)間的影響泉唁。
生存分析統(tǒng)計(jì)方法涉及統(tǒng)計(jì)描述、差異性分析和回歸分析三大類揩慕。其中亭畜，統(tǒng)計(jì)描述主要有Kaplan-Meier估計(jì)法和Life table估計(jì)法；差異性分析主要有對(duì)數(shù)秩檢驗(yàn)和Wilcoxon test迎卤；而回歸分析主要有Cox比例和非比例風(fēng)險(xiǎn)回歸模型拴鸵、參數(shù)回歸模型。
library(survminer)
library(survival)
library(reshape2)
library(tidyverse)
# Load data
############
load("./data/MHC_immunotherapy.RData")
ICB_data <- ICB_study[ICB_study$study == "Liu_2019", ]
ICB_data$highTMB <- ICB_data$TMB > quantile(ICB_data$TMB, 0.5, na.rm = T)
# Save
#########
save.image(file = "./results/data/manuscript_surv_neoAgB.RData")
saveRDS(list(p_TMB = p_TMB, p_pMHC1 = p_pMHC1, p_pMHC2 = p_pMHC2,
p_pMHCnorm = p_pMHCnorm, p_neoAgB1 = p_neoAgB1,
p_neoAgB2 = p_neoAgB2, p_neoAgBnorm = p_neoAgBnorm,
fp_pMHC = fp_pMHC, fp_neoAgB = fp_neoAgB),
file = "./results/data/manuscript_surv_neoAgB.rds")
生存分析4
生存分析是一種統(tǒng)計(jì)方法蜗搔，用于描述劲藐、測(cè)量和分析事件的特征，尋找其發(fā)生的原因碍扔，并對(duì)生存以及事件發(fā)生的時(shí)間進(jìn)行預(yù)測(cè)瘩燥。生存分析的目的在于研究哪些因素影響了事件的發(fā)生速度及生存時(shí)間的長(zhǎng)短秕重，它探索的因變量是一個(gè)包含了刪失或者截除數(shù)據(jù)的事件時(shí)間變量不同。生存分析涉及的事件可以是患者的存活時(shí)間、疾病的復(fù)發(fā)溶耘、企業(yè)的存活時(shí)間等二拐。
生存分析的目的包括但不限于以下幾點(diǎn)：
- 描述生存時(shí)間：描述生存時(shí)間的分布情況，包括生存時(shí)間的中位數(shù)凳兵、平均值等統(tǒng)計(jì)量百新。
- 識(shí)別風(fēng)險(xiǎn)因素：識(shí)別和評(píng)估影響生存時(shí)間的各種風(fēng)險(xiǎn)因素。
- 預(yù)測(cè)生存概率：預(yù)測(cè)個(gè)體在特定時(shí)間點(diǎn)之前生存或事件發(fā)生的概率庐扫。
- 比較生存情況：比較不同群體或不同條件下的生存情況饭望，例如比較不同治療方案的效果仗哨。
- 評(píng)估干預(yù)效果：評(píng)估醫(yī)療干預(yù)、政策或其他措施對(duì)生存時(shí)間的影響铅辞。
生存分析統(tǒng)計(jì)方法涉及統(tǒng)計(jì)描述厌漂、差異性分析和回歸分析三大類。其中斟珊，統(tǒng)計(jì)描述主要有Kaplan-Meier估計(jì)法和Life table估計(jì)法苇倡；差異性分析主要有對(duì)數(shù)秩檢驗(yàn)和Wilcoxon test；而回歸分析主要有Cox比例和非比例風(fēng)險(xiǎn)回歸模型卒密、參數(shù)回歸模型奏纪。
library(survminer)
library(survival)
library(reshape2)
library(tidyverse)
# Load data
############
load("./data/MHC_immunotherapy.RData")
ICB_data <- ICB_study[ICB_study$study == "Liu_2019", ]
ICB_data$highTMB <- ICB_data$TMB > quantile(ICB_data$TMB, 0.5, na.rm = T)
# Save
#########
save.image(file = "./results/data/manuscript_surv.RData")
saveRDS(p_surv_ls, file = "./results/data/manuscript_surv.rds")
差異分析
limma差異分析是一種用于基因表達(dá)數(shù)據(jù)的統(tǒng)計(jì)分析方法枯夜，它主要用于識(shí)別差異表達(dá)的基因。以下是對(duì)limma差異分析的詳細(xì)解釋：
- limma包簡(jiǎn)介： limma（Linear Models for Microarray Data）是一個(gè)R/Bioconductor軟件包综慎，它提供了一套集成的解決方案，用于分析基因表達(dá)實(shí)驗(yàn)數(shù)據(jù)勤庐。limma包含了豐富的功能寥粹，用于處理復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)，并通過(guò)信息借用來(lái)克服樣本量小的問(wèn)題埃元。
- 主要特點(diǎn)：
- 線性模型：limma使用線性模型來(lái)分析基因表達(dá)數(shù)據(jù)涝涤，這包括簡(jiǎn)單復(fù)制設(shè)計(jì)、多組實(shí)驗(yàn)岛杀、直接設(shè)計(jì)阔拳、因子設(shè)計(jì)和時(shí)間序列實(shí)驗(yàn)等。
- 經(jīng)驗(yàn)貝葉斯統(tǒng)計(jì)：limma解釋了經(jīng)驗(yàn)貝葉斯測(cè)試統(tǒng)計(jì)量类嗤，這些統(tǒng)計(jì)量用于在基因表達(dá)分析中獲得后驗(yàn)方差估計(jì)器糊肠。
- 質(zhì)量權(quán)重：limma允許使用質(zhì)量權(quán)重、自適應(yīng)背景校正和控制點(diǎn)與線性模型結(jié)合使用遗锣，以提高小樣本實(shí)驗(yàn)中基因和基因集水平的推斷能力货裹。
- 應(yīng)用范圍：
- 最初，limma主要用于微陣列數(shù)據(jù)的分析精偿，但隨著技術(shù)的發(fā)展弧圆，limma的能力已經(jīng)顯著擴(kuò)展，現(xiàn)在可以執(zhí)行RNA sequencing（RNA-seq）數(shù)據(jù)的差異表達(dá)和差異剪接分析笔咽。
- 這使得之前僅限于微陣列數(shù)據(jù)的下游分析工具現(xiàn)在也可以用于RNA-seq數(shù)據(jù)搔预。
- 分析流程：
- limma分析的主要步驟包括數(shù)據(jù)預(yù)處理、線性模型擬合叶组、經(jīng)驗(yàn)貝葉斯方法應(yīng)用以及結(jié)果的診斷和解釋拯田。
- limma提供了一系列的函數(shù)，用于存儲(chǔ)數(shù)據(jù)或結(jié)果的不同類別甩十，以及在線文檔頁(yè)面船庇，用于每個(gè)單獨(dú)的函數(shù)和每個(gè)主要步驟吭产。
- 與其他方法的比較：
- 在RNA-seq數(shù)據(jù)的差異分析中，limma鸭轮、EdgeR和DESeq2是三種有效的工具垮刹。雖然它們的分析協(xié)議在某些步驟上有所不同（例如，limma使用線性模型進(jìn)行統(tǒng)計(jì)分析张弛，而EdgeR和DESeq2使用負(fù)二項(xiàng)分布）荒典，但它們的結(jié)果部分重疊，每種方法都有其自身的優(yōu)勢(shì)吞鸭。
library(tidyverse)
library(edgeR)
library(limma)
library(fgsea)
load("./data/MHC_immunotherapy.RData")
# Added: mapping table as ENSG -> HGNC must occur here (created by create_ENSG_HGNC_map_table.R)
ENSG_HGNC_map_table <- ENSG_HGNC_map_table %>%
select(gene_id = ensembl_gene_id, gene = external_gene_name, gene_biotype)
# Save
dir.create("./results/temp/", recursive = TRUE)
saveRDS(DGE_ls, "./results/temp/DGE_ls.rds")
功能分析
功能富集分析GSEA（Gene Set Enrichment Analysis）是一種用于解釋基因表達(dá)數(shù)據(jù)的計(jì)算方法寺董，它能夠確定特定基因集是否在某種生物學(xué)狀態(tài)下顯著富集。GSEA的核心思想是刻剥，相比單個(gè)基因的分析遮咖，一群基因（基因集）的協(xié)同變化可能更具有生物學(xué)意義。
GSEA分析的主要步驟包括：
- 計(jì)算富集分?jǐn)?shù)：評(píng)估每個(gè)基因集在排名列表中的表現(xiàn)造虏，富集分?jǐn)?shù)（Enrichment Score, ES）是一個(gè)介于0到1之間的值御吞，用來(lái)衡量基因集的富集程度。
- 估計(jì)富集分?jǐn)?shù)顯著性水平：通過(guò)計(jì)算基因集的p值來(lái)評(píng)估其統(tǒng)計(jì)顯著性漓藕，這通常涉及到使用排列測(cè)試來(lái)確定基因集富集分?jǐn)?shù)的隨機(jī)分布陶珠。
- 矯正多重假設(shè)驗(yàn)證：由于GSEA會(huì)同時(shí)測(cè)試多個(gè)基因集，因此需要對(duì)p值進(jìn)行校正以控制假陽(yáng)性率享钞，常用的方法包括Bonferroni校正和FDR（False Discovery Rate）校正揍诽。
GSEA分析的目的在于：
- 識(shí)別生物過(guò)程中的關(guān)鍵基因集：GSEA能夠識(shí)別在特定條件下顯著富集的基因集，這些基因集可能涉及特定的生物學(xué)過(guò)程或疾病相關(guān)的通路栗竖。
- 揭示基因表達(dá)的總體趨勢(shì)：GSEA不僅關(guān)注差異表達(dá)的基因暑脆，而是利用所有基因的表達(dá)數(shù)據(jù)，揭示整個(gè)基因集的表達(dá)趨勢(shì)狐肢，從而判斷特定通路是被激活還是抑制添吗。
- 提高數(shù)據(jù)分析的可解釋性：通過(guò)將基因集與已知的生物學(xué)知識(shí)相聯(lián)系，GSEA提高了基因表達(dá)數(shù)據(jù)的解釋性份名，使得研究者能夠更好地理解數(shù)據(jù)背后的生物學(xué)意義碟联。
library(tidyverse)
library(ggrepel)
library(fgsea)
library(cowplot)
load("./data/MHC_immunotherapy.RData")
DGE_ls <- readRDS("./results/temp/DGE_ls.rds")
# Save
#########
p_DGE_ls <- list(p_fGSEA = fGSEA_plot, p_RS_ls = p_RS_ls, p_GSEA_preIpi = p_GSEA_preIpi)
save.image(file = "./results/data/manuscript_DGE_analysis.RData")
saveRDS(p_DGE_ls, file = "./results/data/manuscript_DGE_analysis.rds")
介紹
組合多個(gè)文章圖，圖所需要的結(jié)果文件見(jiàn)：數(shù)據(jù)分析：完整的成體系的生存分析的結(jié)果匯總
數(shù)據(jù)下載
數(shù)據(jù)下載鏈接見(jiàn)：
- 百度網(wǎng)盤鏈接: https://pan.baidu.com/s/1KZLVvyrU1Yj78KL0CYYygQ
- 提取碼:
圖所需要的結(jié)果文件見(jiàn)：數(shù)據(jù)分析：完整的成體系的生存分析的結(jié)果匯總
圖1
p_ls <- readRDS("./results/data/manuscript_surv.rds")
p_BR_ls <- readRDS("./results/data/manuscript_RECIST.rds")
# Create fig
p1 <- p_ls$p_TMB +
ggtitle("Tumor Mutation Burden")
p2 <- p_ls$p_strongMHC1 +
ggtitle("MGBS-I")
p3 <- p_ls$p_strongMHC2 +
ggtitle("MGBS-II")
p4 <- p_ls$p_strongMHCnorm +
ggtitle("MGBS-d")
dir.create("./results/figs/", recursive = T)
ggplot2::ggsave(paste0("./results/figs/manuscript_fig1.pdf"), p, width = 178, height = 265 / 1.5, units = "mm")
image.png
圖2
p_ls <- readRDS("./results/data/manuscript_surv.rds")
p1 <- plot_grid(
NA, p_ls$p_ipi_vs_noIpi,
ncol = 2
)
ggplot2::ggsave(paste0("./results/figs/manuscript_fig2.pdf"), p, width = 178, height = 0.5 * 265, units = "mm")
image.png
圖3
p_ls <- readRDS("./results/data/manuscript_multivariate.rds")
# ROC & surv
p_ROC_surv <- plot_grid(
p_ls$LR_model_ls$"TRUE"$p_ROC, p_ls$LR_model_ls$"TRUE"$p_surv, NA,
ncol = 1,
labels = c("a", "c"),
label_size = 8
)
ggplot2::ggsave("./results/figs/manuscript_fig3.pdf", p, width = 178, height = 265 / 2, units = "mm")
image.png
圖4
p_ls <- readRDS("./results/data/manuscript_DGE_analysis.rds")
ggplot2::ggsave("./results/figs/manuscript_fig4.pdf", p, width = 3 / 4 * 178, height = 265 / 3, units = "mm")
image.png
附圖2
p_ls <- readRDS("./results/data/manuscript_surv.rds")
p <- plot_grid(
p_ls$p_th_ls$TMBclass$`0.75` + ggtitle("Tumor Mutation Burden") + theme(legend.background = element_blank()),
p_ls$p_th_ls$TMBclass$`0.9` + ggtitle("Tumor Mutation Burden") + theme(legend.background = element_blank()),
p_ls$p_th_ls$TMBclass$`0.95` + ggtitle("Tumor Mutation Burden") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC1class$`0.75` + ggtitle("MGBS-I") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC1class$`0.9` + ggtitle("MGBS-I") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC1class$`0.95` + ggtitle("MGBS-I") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC2class$`0.75` + ggtitle("MGBS-II") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC2class$`0.9` + ggtitle("MGBS-II") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHC2class$`0.95` + ggtitle("MGBS-II") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHCnormclass$`0.75` + ggtitle("MGBS-d") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHCnormclass$`0.9` + ggtitle("MGBS-d") + theme(legend.background = element_blank()),
p_ls$p_th_ls$MHCnormclass$`0.95` + ggtitle("MGBS-d") + theme(legend.background = element_blank()),
ncol = 3
)
ggplot2::ggsave(paste0("./results/figs/manuscript_figS2_surv_th.pdf"), p, width = 178, height = 265, units = "mm")
image.png
附圖3
p <- readRDS("./results/data/manuscript_correlation_analysis.rds")
# Create figure
ggplot2::ggsave(paste0("./results/figs/manuscript_figS3_correlation.pdf"), p, width = 178, height = 265, units = "mm")
image.png
附圖4
p_ls <- readRDS("./results/data/manuscript_validation.rds")
# Forest plot
# Adjust sizes facet
gt <- ggplot_gtable(ggplot_build(p_ls$p_val_ls$fp))
# gtable_show_layout(gt)
gt$heights[10] <- 1 / 2 * gt$heights[10]
# fp<- grid.draw(gt)
# Forest plot preIpi
fp_preIpi <- p_ls$p_val_ls$fp_preIpi
# Create figure
ggplot2::ggsave(paste0("./results/figs/manuscript_figS4_validation.pdf"), p, width = 178, height = 265, units = "mm")
image.png
附圖5
p_ls <- readRDS("./results/data/manuscript_surv_neoAgB.rds")
p <- plot_grid(
p_ls$fp_neoAgB + ggtitle("Absolute Neoantigen Burden") + theme(title = element_text(size = 8)), p_ls$fp_pMHC + ggtitle("Normalized Neoantigen Burden") + theme(title = element_text(size = 8)),
p_ls$p_neoAgB1 + ggtitle("MHC-I NeoAgB"), p_ls$p_pMHC1 + ggtitle("norm. MHC-I NeoAgB"),
p_ls$p_neoAgB2 + ggtitle("MHC-II NeoAgB"), p_ls$p_pMHC2 + ggtitle("norm. MHC-II NeoAgB"),
p_ls$p_neoAgBnorm + ggtitle("diff. NeoAgB"), p_ls$p_pMHCnorm + ggtitle("norm. diff. NeoAgB"),
ncol = 2
)
# Create figures
ggplot2::ggsave(paste0("./results/figs/manuscript_figS5_neoAgB.pdf"), p_neoAgB, width = 178, height = 265, units = "mm")
image.png
附圖6-7
p_ls <- readRDS("./results/data/manuscript_multivariate.rds")
# 1. COX
t_preIpi <- p_ls$Cox_model_ls$"TRUE"$t_uni
t_ipiNaive <- p_ls$Cox_model_ls$"FALSE"$t_uni
p1 <- plot_grid(
t_ipiNaive, t_preIpi,
ncol = 2
)
p_cox <- plot_grid(
p1, p_ls$Cox_model_ls$fp_multi + scale_y_continuous(limits = c(0.15, 12), trans = "log10"), NA,
ncol = 1,
rel_heights = c(1, 1 / 2, 1)
)
# Create figures
ggplot2::ggsave(paste0("./results/figs/manuscript_figS6_LogReg.pdf"), p_logReg, width = 178, height = 265, units = "mm")
ggplot2::ggsave(paste0("./results/figs/manuscript_figS7_Cox.pdf"), p_cox, width = 178, height = 265, units = "mm")
image.png
附圖8
p_ls <- readRDS("./results/data/manuscript_validation_multivariate.rds")
p_compare_ls <- readRDS("./results/data/manuscript_validation_compare_variables.rds")
# Create figure
ggplot2::ggsave(paste0("./results/figs/manuscript_figS8_validation_multivariate.pdf"), p, width = 178, height = 265, units = "mm")
image.png
附圖9
p_ls <- readRDS("./results/data/manuscript_DGE_analysis.rds")
ggplot2::ggsave(paste0("./results/figs/manuscript_figS9_GSEA.pdf"), p, width = 178, height = 265, units = "mm")
image.png
系統(tǒng)信息
- 數(shù)據(jù)分析
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Asia/Shanghai
tzcode source: internal
attached base packages:
[1] grid stats graphics grDevices utils datasets methods base
other attached packages:
[1] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
[7] tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0 gridExtra_2.3 ggplot2_3.5.1 cowplot_1.1.3
loaded via a namespace (and not attached):
[1] vctrs_0.6.5 cli_3.6.3 rlang_1.1.4 stringi_1.8.3 generics_0.1.3
[6] glue_1.7.0 colorspace_2.1-0 hms_1.1.3 scales_1.3.0 fansi_1.0.6
[11] munsell_0.5.0 tzdb_0.4.0 lifecycle_1.0.4 compiler_4.3.3 timechange_0.3.0
[16] pkgconfig_2.0.3 rstudioapi_0.16.0 R6_2.5.1 tidyselect_1.2.1 utf8_1.2.4
[21] pillar_1.9.0 magrittr_2.0.3 tools_4.3.3 withr_3.0.0 gtable_0.3.4
- 畫圖
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Asia/Shanghai
tzcode source: internal
attached base packages:
[1] grid stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggrepel_0.9.5 fgsea_1.28.0 edgeR_4.0.16 limma_3.58.1 scales_1.3.0 pROC_1.18.5
[7] caret_6.0-94 lattice_0.22-6 reshape2_1.4.4 survival_3.7-0 survminer_0.4.9 ggpubr_0.6.0
[13] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
[19] tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0 gridExtra_2.3 ggplot2_3.5.1 cowplot_1.1.3
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 timeDate_4032.109 digest_0.6.35 rpart_4.1.23
[5] timechange_0.3.0 lifecycle_1.0.4 statmod_1.5.0 magrittr_2.0.3
[9] compiler_4.3.3 rlang_1.1.4 tools_4.3.3 utf8_1.2.4
[13] data.table_1.15.2 knitr_1.45 ggsignif_0.6.4 plyr_1.8.9
[17] BiocParallel_1.36.0 abind_1.4-5 withr_3.0.0 nnet_7.3-19
[21] stats4_4.3.3 fansi_1.0.6 xtable_1.8-4 colorspace_2.1-0
[25] future_1.33.1 globals_0.16.3 iterators_1.0.14 MASS_7.3-60.0.1
[29] cli_3.6.3 generics_0.1.3 rstudioapi_0.16.0 future.apply_1.11.1
[33] km.ci_0.5-6 tzdb_0.4.0 splines_4.3.3 parallel_4.3.3
[37] survMisc_0.5.6 vctrs_0.6.5 hardhat_1.3.1 Matrix_1.6-5
[41] carData_3.0-5 car_3.1-2 hms_1.1.3 rstatix_0.7.2
[45] listenv_0.9.1 locfit_1.5-9.9 foreach_1.5.2 gower_1.0.1
[49] recipes_1.0.10 glue_1.7.0 parallelly_1.37.1 codetools_0.2-19
[53] stringi_1.8.3 gtable_0.3.4 munsell_0.5.0 pillar_1.9.0
[57] ipred_0.9-14 lava_1.8.0 R6_2.5.1 KMsurv_0.1-5
[61] backports_1.4.1 broom_1.0.5 class_7.3-22 fastmatch_1.1-4
[65] Rcpp_1.0.12 nlme_3.1-164 prodlim_2023.08.28 xfun_0.43
[69] zoo_1.8-12 pkgconfig_2.0.3 ModelMetrics_1.2.2.2
