今天为大家详细复现一篇于今年8月发表在《 International Immunopharmacology》（IF=4.8，中科院二区）杂志上的纯生信文章，本文用到的生信分析数据都来自公共数据库，本文研究思路如下：

首先使用差异表达基因分析寻找两疾病（PD和pSS）各自的差异基因，并且在两疾病的差异表达基因之间取交集；使用加权基因共表达网络(WGCNA)寻找与各自疾病高度相关的模块和基因，并且在两疾病之间取交集基因。对差异表达分析的交集结果和WGCNA的交集结果再取交集，获得了与两疾病高度相关且差异表达的共同基因。

随后，对共同基因进行了GO和KEGG富集分析。使用三种机器学习算法(LASSO、支持向量机递归特征消除(SVM-RFE)和随机森林(RF))进一步筛选，获得了两疾病的潜在的生物标志(CSF2RB, CXCR4, LYN)，为验证其诊断价值，分别在发现集和验证集上用ROC曲线分析了它们的诊断效力。

最后通过CIBERSORT进行免疫细胞浸润分析；NetworkAnalyst分析转录因子(TFs)-基因和miRNAs-基因调控网络；DSigDB预测相关的药物靶点，鉴定了PD和pSS之间的共同分子机制、转录因子(TFs), miRNAs 以及候选药物。

复写图片集锦可见：

数据来源和研究路线图

--- ·4个数据集· ---

文献：

研究路线：

插一段优质广告，收藏备用以备不时之需！CNSknowall公司位于上海市陆家嘴，汇集了几十位涵盖所有组学分析的顶尖生信硕博工程师，不仅拥有最好用、规模最大的在线数据分析云平台，同时也拥有强大的生信课题设计和定制分析能力，更拥有极致的性价比和服务到发表的完美售后：
0 准备工作

0.1 R包载入

没有的包可以通过：

多数情况下我们通过BiocManager::install()安装，提示没找到的，去搜索引擎检索一下它的来源和安装方式。

library(GEOquery) #读入GEO数据
library(oligo) #Affymatrix的原始数据（.CEL文件）读入方法1
library(affy) #Affymatrix的原始数据（.CEL文件）读入方法2
library(limma) #芯片数据差异表达分析
library(sva) #ComBat函数矫正批次效应
library(tidyverse) #载入该包时，ggplot2, dplyr, stringr等包也会被一同载入
# library(ggplot2) #可视化
# library(dplyr) #管道符及数据清洗
# library(stringr) #对字符串进行精细操作
library(FactoMineR) #主成分分析图需要加载这两个包
library(factoextra)
library(patchwork) #拼图
library(pheatmap) #绘制热图
library(grDevices) #用于图形设备和图形输出
library(ggplotify) #as.ggplot函数
library(VennDiagram) #韦恩图绘制
library(org.Hs.eg.db) #人类的org.db包
library(AnnotationDbi) #富集分析id转换
library(clusterProfiler) #富集分析
library(enrichplot) #富集分析可视化
library(WGCNA) #WGCNA加权基因共表达网络分析
library(glmnet) #lasso回归
library(e1071) #svm-ref
library(caret)
library(randomForest) #随机森林
library(rstatix) #用于统计
library(ggdist) #用于在可视化中标注置信区间
library(Hmisc) #Hmisc包与rms包，用于统计分析和回归建模
library(rms)
library(pROC) #计算ROC曲线
library(CIBERSORT) #免疫浸润
library(magrittr) #管道符%<>%
library(ggpubr) #另一个在gplot2图形中添加显著性的方法
library(Ipaper)
library(corrplot) #绘制相关性矩阵图

0.2 数据下载

可以用迅雷或者GEOquery包的getGEO函数 ，因为我开了迅雷会员，下的比函数快不少，所以复制连接打开迅雷下载。

# 清空工作环境
rm(list = ls())
options(stringsAsFactors = F) #创建数据框时，不自动将字符串转换为因子类型

# 创建画图和保存数据的文件夹
dir.create("./plots")
dir.create("./data")

# 读入GEO数据
GEOquery
# 列出文件夹下所有文件
dir()
# 牙周炎PD数据集
gse1 <- getGEO(filename = "./data/GSE16134_series_matrix.txt.gz", getGPL = F)
# 原发性干燥综合征pSS数据集
gse2 <- getGEO(filename = "./data/GSE40611_series_matrix.txt.gz", getGPL = F)

# 建议一下进来就读取临床信息，在数据质控过滤过程中临床信息很重要
pd1 <- pData(gse1)
pd2 <- pData(gse2)

0.3 原始数据处理

芯片厂商是affymatrix，使用affy包进行.CEL文件原始数据的读入和处理

library(affy)
# 读入.CEL文件
raw1 <- ReadAffy(celfile.path = "./GSE16134_RAW")
# rma一步进行背景矫正，标准化
exp1 <- rma(raw1)
# 从ExpressionSet对象中提取表达矩阵
exp1 <- exprs(exp1)

raw2 <- ReadAffy(celfile.path = "./GSE40611_RAW")
exp2 <- rma(raw2)
exp2 <- exprs(exp2)

save(exp1, exp2, file = "./data/step0_Affy_raw.Rdata")

1 表达矩阵质控过滤

对从.CEL原始数据得到的表达矩阵进行过滤，id转换

load("./data/step0_Affy_raw.Rdata")
range(exp1)
range(exp2)

# 箱线图看一下各样本表达量分布方差齐性
boxplot(exp1)
boxplot(exp2)

# pSS数据集可以做下批次效应矫正
# library(sva)
# exp2 <- ComBat(exp2, batch = pd2$`batch:ch1`)
# boxplot(exp2)

#做一下归一化，消除不同样本间的技术性偏差，使得各样本的数据更加可比
exp1 <- normalizeBetweenArrays(exp1)
exp2 <- normalizeBetweenArrays(exp2)

# 探针id转换成gene symbol
# 芯片厂商是昂飞：Affymetrix Human Genome U133 Plus 2.0 Array，比较常见
# 从网页https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570下载GPL570-55999.txt表格
gpl <- read.csv("GPL570-55999.txt", sep = "\t", header = T, comment.char = "#")
id2s <- data.frame(ID = gpl$ID, SYMBOL = gpl$Gene.Symbol)
id2s[1:4,]
# 可以看到存在一个探针对应多个symbol的情况，有两种处理方法，一是保留第一个对应的，二是删除所有一个探针对应多个基因的行。这里我们采取第二种方法
id2s <- id2s[!str_detect(id2s$SYMBOL, "///"), ]
# 接下来转化表达矩阵行名，因为存在多个探针对应一个基因的情况，我们按照文献计算多个探针的均值

# 写一个id转换的函数
IDtrans <- function(id2s, exp){
  
  # 让表达矩阵与id2s的探针ID顺序一致
  exp <- exp[match(id2s$ID, rownames(exp)),]
  rownames(exp) <- id2s$SYMBOL
  exp <- as.data.frame(exp)
  
  exp$RowName <- id2s$SYMBOL
  
  # 去除空列或者na列
  exp <- exp %>%
  filter(!is.na(RowName) & RowName != "")

  # 使用 dplyr 进行分组和计算平均值
  result <- exp %>%
    group_by(RowName) %>%
      summarise(across(everything(), mean))

  result <- as.data.frame(result)
  # 将行名设置为分组名称
  rownames(result) <- result$RowName
  result <- result[,-1]

  return(result)
}

# 然后分别对exp1和exp2执行，注意，因为这两个表达矩阵的GPL平台都是GPL570，才能共用一套探针id转换规则（整理得到的id2s数据框）
exp1 <- IDtrans(id2s, exp1)
exp2 <- IDtrans(id2s, exp2)

dim(exp1); dim(exp2)

save(exp1, exp2, file = "./data/Step1_transid_exp.RData")

2 临床信息读取整理

将临床信息与表达矩阵中的样本一一对应

# 读取临床信息
pd1 <- pData(gse1)
pd2 <- pData(gse2)

# 先看一下表达矩阵和临床信息有无格式的差别
colnames(exp1)
rownames(pd1)

# 点号（.）在正则表达式中是一个特殊字符，表示任意单个字符，因此需要使用双反斜杠（\\.）进行转义。
colnames(exp1) <- str_split(colnames(exp1), "\\.", simplify = T)[,1]

colnames(exp2)
rownames(pd2)

colnames(exp2) <- str_split(colnames(exp2), "_", simplify = T)[,1]

# 看下临床信息是否与表达矩阵一一对应
identical(rownames(pd1), colnames(exp1))
identical(rownames(pd2), colnames(exp2))

# 提取其中的分组信息

table(pd1$`gingival tissue phenotype:ch1`)
# Diseased (BoP + interproximal PPD of>3 mm + concomitant CAL of>2 mm) 
# GPT翻译：
#（患病）：牙龈出血（BoP）+ 龈间袋探诊深度（PPD）> 3 毫米 + 同时存在的临床附着丧失（CAL）> 2 毫米
# Healthy (no BoP + PPD of<5 mm + CAL of<3 2 mm
#（健康）：无牙龈出血（BoP）+ 探诊深度（PPD）< 5 毫米 + 临床附着丧失（CAL）< 2 毫米
group1 <- ifelse(str_detect(pd1$`gingival tissue phenotype:ch1`, "Diseased"), "Diseased", "Healthy")

table(pd2$source_name_ch1)
# parotid tissue from control (Control)
# 来自对照组的腮腺组织
# parotid tissue from primary Sjogren syndrome patient (pSS)
# 来自原发性干燥综合征患者的腮腺组织
# parotid tissue from sicca patient (Sicca)
# 来自干燥症患者的腮腺组织
group2 <- pd2$`disease status:ch1`

# 筛选：这里只保留pss和control
filt <- group2 != "Sicca"
group2 <- group2[filt]
exp2 <- exp2[,filt]
length(group2); ncol(exp2)

# 最后将group类型改为factor，并且设定levels是实验组在前，对照组在后，以方便后续进行差异表达分析（DEG）
table(group1)
group1 <- factor(group1, levels = c("Diseased", "Healthy"))
table(group2)
group2 <- factor(group2, levels = c("pSS", "Control"))

save(exp1, group1, pd1,
     exp2, group2, pd2,
     file = "./data/step2_exp&group.RData")

3 质控-PCA主成分分析
    主要看下各分组的样本是否能聚在一起，这张图在文献中并没有展现，是我们常规进行质控需要看的图。

#load("./data/step2_exp&group.RData")

# 定义一个用于画PCA图的函数qc_PCA
# 输入数据形式要求exp是列为样本，行为基因的表达矩阵；group是包含实验-对照分组信息的向量
qc_PCA <- function(exp, group){
  cg = names(tail(sort(apply(exp,1,sd)),1000))
  # apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
  
  # 下面是画PCA的必须操作，需要看说明书。
  dat = t(exp[cg,])#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
  dat = as.data.frame(dat)#将matrix转换为data.frame 
  
  dat.pca <- PCA(dat , graph = FALSE)
  # 现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  p <- fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind =  group, # color by groups 
               addEllipses = T,
               legend.title = "Groups"
  )
  return(p)
}

# 调用函数绘制PCA图
p1 <- qc_PCA(as.matrix(exp1), group1) + ggtitle("PD_PCA")
p2 <- qc_PCA(exp2, group2) + ggtitle("pSS_PCA")

# 拼图
p1 + p2

# 保存图片
ggsave('./plots/QC_all_samples_PCA.png', p1 + p2, height = 5, width = 8)

# 我们通常使用limma包进行芯片数据的差异表达分析
# 定义一个进行差异表达分析的函数doDEG
# 这里vs参数是设置实验组和对照组的，若缺省，则默认使用group因子的第一个level作为实验，第二个level作为对照
doDEG <- function(dat, group, logFC_t = NULL, padj_t = 0.05, vs = NULL){
  design <- model.matrix(~0 + group)
  colnames(design)=levels(group)
  rownames(design)=colnames(dat)
  if(is.null(vs)){
    vs <- paste0(levels(group)[1], "-", levels(group)[2])
  }

  contrast.matrix<-makeContrasts(contrasts = vs,levels = design)

  fit <- lmFit(dat,design)
  fit2 <- contrasts.fit(fit, contrast.matrix)
  fit2 <- eBayes(fit2)
  deg = topTable(fit2, coef=1, n=Inf)
  deg = na.omit(deg)
  if(is.null(logFC_t)){
    logFC_t <- with(deg,mean(abs(deg$logFC)) + 2*sd(abs(deg$logFC)) )
  }
  print(paste0("The threshold of logFoldChange is ", logFC_t))
  deg$v= -log10(deg$adj.P.Val)
  deg$change <- ifelse(deg$adj.P.Val > padj_t,'stable',
                       ifelse(deg$logFC > logFC_t,'up',
                              ifelse(deg$logFC < -logFC_t,'down','stable'))
  )
  deg$name=rownames(deg)
  return(deg)
}

# 使用0.5作为logFoldChange的阈值，0.05作为p.adjust的阈值进行差异表达分析
DEG1 <- doDEG(exp1, group1, logFC_t = 0.5)
table(DEG1$change)

DEG2 <- doDEG(exp2, group2, logFC_t = 0.5)
table(DEG2$change)
save(DEG1, DEG2, file = "./data/step3_deg.Rdata")

火山图绘制代码：

文章原图Fig.2A, 2B（上）vs CNSknowall R语言复现（下）

解析：复现结果中，差异表达基因的数量和文献关注的三个基因的logFC和校正后p值在火山图可视化结果中与原文几乎一致。

图2 C+D+E：根据PD和pss的所有差异表达基因，分别绘制了聚类热图(图2C和D)。利用维恩图，PD和pss之间共得到173个共同DEGs(图2E)。

复现代码：

热图绘制代码：

#定义绘制热图的函数
DEG_heatmap <- function(dat, group, target_gene) 
{
  n = dat
  ac = data.frame(group = group)
  rownames(ac) = colnames(n)
  color <- c("#1F77B4FF", "#FF7F0EFF", "#2CA02CFF", "#D62728FF", 
             "#FF9896FF", "#9467BDFF", "#C5B0D5FF", "#8C564BFF", "#98DF8AFF", 
             "#C49C94FF", "#E377C2FF", "#FFBB78FF", "#F7B6D2FF", "#7F7F7FFF", 
             "#C7C7C7FF")[1:length(levels(group))]
  names(color) <- names(table(group))

  p <- pheatmap(n,
                color = colorRampPalette(c("#00008B", "#555597", "#FFFFB0", 
                                           "#B1553A", "#8B0000"))(100),
                labels_row = target_gene, 
                show_colnames = F,
                show_rownames = T, 
                cluster_cols = F,
                annotation_colors = list(group = color), 
                annotation_col = ac)
  return(p)
}

# 这里是为了显示target_gene进行的
# 选出差异表达基因绘制热图
plot_exp1 <- plot_exp1[DEG1[DEG1$change!="stable", "name"],]
# 显示target_gene作为行名
rownames_custom1 <- ifelse(rownames(plot_exp1) %in% target_gene,
                          rownames(plot_exp1), "")
# 选出差异表达基因绘制热图
plot_exp2 <- plot_exp2[DEG2[DEG2$change!="stable", "name"],]
# 显示targetgene作为行名
rownames_custom2 <- ifelse(rownames(plot_exp2) %in% target_gene, rownames(plot_exp2), "")

# 调用函数绘制热图
p1 <- DEG_heatmap(plot_exp1, plot_group1, rownames_custom1)
p2 <- DEG_heatmap(plot_exp2, plot_group2, rownames_custom2)

ggsave(as.ggplot(p1) + as.ggplot(p2),
       filename = "./plots/Fig2_DEG_heatmap.png",
       width = 12,       height = 6)

韦恩图绘制代码：

PD_up <- DEG1[DEG1$change == "up","name"]
PD_down <- DEG1[DEG1$change == "down","name"]
pSS_up <- DEG2[DEG2$change == "up","name"]
pSS_down <- DEG2[DEG2$change == "down","name"]

#R Markdown 有时会在处理图形设备时出现问题，特别是在 png 和 dev.off() 之间。
#复制粘贴在左下窗口跑就行

png(filename = "./plots/Fig2_venn.png", width = 480, height = 480)
venn.plot <- venn.diagram(
  x = list(Set1 = PD_up, Set2 = pSS_down, Set3 = pSS_up, Set4 = PD_down),
  category.names = c("PD_up", "pSS_down", "pSS_up", "PD_down"),
  filename = NULL,
  output = TRUE,
  cex = 3,    # 调整交集数字的大小
  fill = c("#D62728FF", "#2CA02CFF", "#FF7F0EFF", "#1F77B4FF"),
  cat.cex = 1.5,  # 调整集合名称的大小
  cat.col = c("#D62728FF", "#2CA02CFF", "#FF7F0EFF", "#1F77B4FF")
)
grid.draw(venn.plot)dev.off()

复现结果：

图2 C+D：

文章原图Fig.2A, 2B（上）vs CNSknowall R语言复现（下）

解析：复现结果中，三个关键基因都在实验组高表达，对照组中低表达

图2 E：

解析：两疾病差异表达的上下调基因取交集结果，CNSknowall的复现结果与原文几乎一致

5 富集分析

复现代码：

#PD和pSS共有的DEG
commonDEGs <- intersect(DEG1[DEG1$change!="stable", "name"],
                        DEG2[DEG2$change!="stable", "name"])

#使用AnnotationDbi包将symbol转成entrezid用于富集分析
library(org.Hs.eg.db)
library(AnnotationDbi)
keytypes(org.Hs.eg.db)
k <- AnnotationDbi::keys(org.Hs.eg.db, keytype = "ENTREZID") 
e2s <- AnnotationDbi::select(org.Hs.eg.db, keys = k,
                             columns="SYMBOL", keytype = "ENTREZID")
ENTREZID_commonDEGs <- e2s[match(commonDEGs, e2s$SYMBOL),"ENTREZID"] %>% na.omit %>% as.numeric

library(clusterProfiler)
library(enrichplot)
#写一个循环，分别对go数据库的BP, MF, CC三个ontology做ORA富集分析
ont <- c("BP", "MF", "CC")
for(i in ont){
  go <- enrichGO(ENTREZID_commonDEGs,
               OrgDb = "org.Hs.eg.db",
               ont = i,
               qvalueCutoff = 0.2,
               pvalueCutoff = 0.05)
  p <- enrichplot::dotplot(go, showCategory = 10) +
    labs(title = i) +
    theme(plot.title = element_text(size = 20))
  ggsave(p, filename = paste0("./plots/Fig3_ORA_GO_", i, ".png"))
}

复现结果：

图3 A+B+C：

文章原图Fig.2A, 2B（左）vs CNSknowall R语言复现（右）

解析：复现结果与原文献略有不一致，推测原因为：

    ①上游数据处理时有差异导致找到的common_DEGs有一定差异

    ②GO数据库版本变化导致通路中基因种类数量可能存在差异。

    但是富集的结果基本都与炎症、细胞因子和免疫相关过程相关，符合文献的生物学结论。

复现代码：

#对KEGG数据库进行ORA富集分析
kegg <- enrichKEGG(ENTREZID_commonDEGs,
                   organism = "hsa")
p <- enrichplot::dotplot(kegg, showCategory = 10) +
    labs(title = "KEGG enrichment") +
    theme(plot.title = element_text(size = 20))
ggsave(p, filename = paste0("./plots/Fig3_ORA_KEGG.png"))

复现结果：

图3 D：

文章原图Fig.2A, 2B（左）vs CNSknowall R语言复现（右）

图4 A+B：PD (GSE16134)和pss (GSE40611)的最佳软功率B分别为14和28(图4A和B)。

复现代码：

load("./data/step2_exp&group.RData")
library(WGCNA)

#该函数输入为选择的方差最大的前多少基因和表达矩阵
Pick_SFT <- function(top_sd_gene_number = 4000,
                     exp,
                     name,
                     RsquaredCut = 0.85){
  #这里我们按文献取方差前10000基因做WGCNA
  apply(exp, 1, sd) %>% sort(.,decreasing = T) %>% head(., top_sd_gene_number) %>% names %>% exp1[.,] -> WGCNA_exp
  
  t_exp <- t(WGCNA_exp)
  
  #seq生成12，14，16，18，20的等差数列
  powers = c(c(1:10), seq(from = 12, to=30, by=2))
  # Call the network topology analysis function
  #选择合适的软阈值
  sft = pickSoftThreshold(t_exp, 
                          RsquaredCut = RsquaredCut,
                          powerVector = powers, 
                          verbose = 5)
  #po是软阈值
  po <- sft$powerEstimate
  assign(paste0("po_", name), po, envir = .GlobalEnv)
  print(po)
  
  png(paste0("./plots/Fig4_pick_SFT_", name, "RsquaredCut_", RsquaredCut,".png"), width = 900, height = 500)
  
  # Plot the results:
  par(mfrow = c(1,2));
  cex1 = 0.8;
  
  # Scale-free topology fit index as a function of the soft-thresholding power
  # 无标度拓扑拟合指数R^2
  plot(sft$fitIndices[,1], 
       -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       xlab="Soft Threshold (power)",
       ylab="Scale Free Topology Model Fit,signed R^2",
       type="n",
       main = paste("Scale independence"));
  text(sft$fitIndices[,1], 
       -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       labels=powers,cex=cex1,col="red");
  
  # this line corresponds to using an R^2 cut-off of h
  abline(h=0.90,col="red")
  
  # Mean connectivity as a function of the soft-thresholding power
  plot(sft$fitIndices[,1], 
       sft$fitIndices[,5],
       xlab="Soft Threshold (power)",
       ylab="Mean Connectivity", 
       type="n",main = paste("Mean connectivity"))
  text(sft$fitIndices[,1], 
       sft$fitIndices[,5], 
       labels=powers, cex=cex1,col="red")
  
  dev.off()
}

Pick_SFT(top_sd_gene_number = 10000, exp1, name = "exp1", RsquaredCut = 0.9)
Pick_SFT(top_sd_gene_number = 10000, exp2, name = "exp2", RsquaredCut = 0.9)

复现结果：

图4 A+B：

解析：两张图的横坐标都是软阈值β，绘制其的目的是指导选择合适的软阈值来构建网络。

    左图的纵坐标是无标度拓扑拟合指数，其值越接近1越像无标度拓扑网络，生物的基因网络不可能完全满足无标度拓扑网络，因此我们可以选取R^2＞0.85（这里我们用0.90，这是人为划定的阈值，可以根据自己的数据进行调整）之后的所有β值进行网络构建。

    右图的纵坐标是基因间的连通性，选取的β值越大就越有可能得到基因之间高度相关的模块，同时也可能丢掉与模块相关性略差一点的基因。（简而言之：β值过小会导致与模块无关的基因被纳入模块，β值过大会导致与模块相关的基因被剔出模块）

    这里我们的复现方式选择使用较大的β值分步构建加权基因共表达网络，使用动态树切割得到模块，之后对较小的模块进行合并。

6.1 PD组的WGCNA分析

图4 C+E：在PD中，共鉴定出10个模块，其中绿松石模块与PD的正相关性最强(相关系数= 0.66,p = 2e-39)，包含1192个基因(图4C和E)。

图4 D+F：在pss中鉴定出6个模块，其中青色模块(相关系数= 0.76,p = 1e-07)显示出最高的正相关性，包含211个基因(图4D和F)。

图4 G+H：此外，散点图显示了PD的绿松石模块和pss的青色模块中模块隶属度与基因显著性之间的相关性(图4G和H)。

复现代码：

# 重新选之前pickSoftThreshold用到的基因：
apply(exp1, 1, sd) %>% sort(.,decreasing = T) %>% head(., 10000) %>% names %>% exp1[.,] -> WGCNA_exp1
t_exp1 <- t(WGCNA_exp1)

allowWGCNAThreads()
#po_exp1 = 22
#共表达相似性和邻接性
adjacency = adjacency(t_exp1, power = po_exp1)
#TOM矩阵
TOM = TOMsimilarity(adjacency)
save(TOM, file = paste0("./data/step6_TOM_exp1_po", po_exp1, ".RData"))
load(paste0("./data/step6_TOM_exp1_po", po_exp1, ".RData"))
dissTOM = 1 - TOM

#根据TOM矩阵构建层次聚类
geneTree = hclust(as.dist(dissTOM), method = "average")
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)

# 模块识别使用动态树切割，这里敏感度deepSplit我设置为2
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 4, pamRespectsDendro = FALSE, minClusterSize = 50) 
# 查看模块数
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
plotDendroAndColors(geneTree,
                    dynamicColors,
                    "Dynamic Tree Cut",
                    dendroLabels = FALSE,
                    hang = 0.03,
                    addGuide = TRUE,
                    guideHang = 0.05,
                    main = "Gene dendrogram and module colors")

# Calculate eigengenes
MEList = moduleEigengenes(t_exp1, colors = dynamicColors)
MEs = MEList$eigengenes
#  Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)
#  Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")
# 画图
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

# 在树状图中加入切割线
abline(h=0.18, col = "red")
# 调用自动归并函数
merge = mergeCloseModules(t_exp1, dynamicColors, cutHeight = 0.18, verbose = 3)
# 合并模块的颜色
mergedColors = merge$colors
table(mergedColors)
# 合并后新模块
mergedMEs = merge$newMEs

png("./plots/Fig4_dendrogram_and_module_colors_exp1.png")
plotDendroAndColors(geneTree,
                    cbind(dynamicColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE,
                    hang = 0.03,
                    addGuide = TRUE,
                    guideHang = 0.05)
dev.off()

# 对模块颜色重命名
moduleColors = mergedColors
# 将颜色转换为数据标签
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder)-1
MEs = mergedMEs
# 保存模块颜色和标签，以供后续部分使用
save(MEs, moduleLabels, moduleColors, geneTree, file = "./data/step6_networkConstruction_stepByStep_exp1.RData")

复现代码：

图4 C：

文章原图（左） vs CNSknowall R语言复现（右）

解析：若下一步中未找到与疾病高度相关的模块，则可以适当减小切割动态树和合并模块的阈值 "cutHeight" 参数，更细致的划分模块从而得到更小的模块。若下一步中找到多个与疾病高度相关的模块，则可以通过层次聚类树查看这些模块能否合并，增加切割动态树和合并模块的阈值 "cutHeight" 参数，合并相似的模块从而得到更大的模块。

复现代码：

# Modules-traits relationships----
# Define numbers of genes and samples
nGenes = ncol(t_exp1)
nSamples = nrow(t_exp1)

# 这里使用merge过后的模块识别eigengenes

#使用分组信息
design=model.matrix(~0+ group1)
colnames(design)=levels(group1)
moduleTraitCor = cor(MEs, design, use = "p")
head(moduleTraitCor)

moduleTraitPvalue = corPvalueStudent(moduleTraitCor, 
                                     nSamples)
head(moduleTraitPvalue)

#可视化

# Will display correlations and their p-values
textMatrix =  paste(signif(moduleTraitCor, 2), 
                    "\n(",signif(moduleTraitPvalue, 1),
                    ")", 
                    sep = "");
dim(textMatrix) = dim(moduleTraitCor)
textMatrix[1:6,1:2]

design1 <- as.data.frame(design)
# Display the correlation values within a heatmap plot
png("./plots/Fig4_module_traits_relationship_exp1.png", width = 480, height = 720)
par(mar = c(5, 10, 4, 2));
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design1),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))
dev.off()


#2.Module membership----
# Define variable weight containing the weight column of datTrait
trait = as.data.frame(design1$Diseased);
names(trait) = "trait"

# names (colors) of the modules
modNames = substring(names(MEs), 3)

geneModuleMembership = as.data.frame(cor(t_exp1, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")

geneTraitSignificance = as.data.frame(cor(t_exp1, trait, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))

names(geneTraitSignificance) = paste("GS.", names(trait), sep="")
names(GSPvalue) = paste("p.GS.", names(trait), sep="")

#选择与PD性状相关性最高的模块
module = "turquoise"
column = match(module, modNames);
moduleGenes = mergedColors == module;

png("./plots/Fig4_GS_MM_exp1.png")
par(mar = c(4,4,4,4));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for body weight",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()

#提取该模块中的基因
PD_module <- colnames(t_exp1)[mergedColors == "turquoise"]
PD_module
save(PD_module, file = "./data/step6_WGCNA_PD_module.RData")

复现结果：

图4 E+G：

文章原图（左） vs CNSknowall R语言复现（右）

解析：热图结果展现出一个与PD性状高度正相关的模块（文献中是turquoise模块，复现结果也是turquoise模块），且复现结果的相关性与文献差距很小。（复现0.64-文献0.66）

    散点图结果展现了MM（基因与特定模块的关联性）与 GS （基因与疾病性状的关联程度）之间的相关性。

    可以看出复现结果中的模块，其MM和GS的相关性更强（复现0.82>文献0.72）。

6.2 pSS组的WGCNA分析

复现代码：

# 重新选之前pickSoftThreshold用到的基因：
apply(exp2, 1, sd) %>% sort(.,decreasing = T) %>% head(., 10000) %>% names %>% exp2[.,] -> WGCNA_exp2
t_exp2 <- t(WGCNA_exp2)

allowWGCNAThreads()
#po_exp2 = 14
#共表达相似性和邻接性
adjacency = adjacency(t_exp2, power = po_exp2)
#TOM矩阵
TOM = TOMsimilarity(adjacency)
save(TOM, file = paste0("./data/step6_TOM_exp2_po", po_exp2, ".RData"))
load(paste0("./data/step6_TOM_exp2_po", po_exp2, ".RData"))
dissTOM = 1 - TOM

#根据TOM矩阵构建层次聚类
geneTree = hclust(as.dist(dissTOM), method = "average")
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)

# 模块识别使用动态树切割，这里敏感度deepSplit我设置为2
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = 50) 
# 查看模块数
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
plotDendroAndColors(geneTree,
                    dynamicColors,
                    "Dynamic Tree Cut",
                    dendroLabels = FALSE,
                    hang = 0.03,
                    addGuide = TRUE,
                    guideHang = 0.05,
                    main = "Gene dendrogram and module colors")

# Calculate eigengenes
MEList = moduleEigengenes(t_exp2, colors = dynamicColors)
MEs = MEList$eigengenes
#  Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)
#  Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")
# 画图
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

# 在树状图中加入切割线
abline(h=0.25, col = "red")
# 调用自动归并函数
merge = mergeCloseModules(t_exp2, dynamicColors, cutHeight = 0.25, verbose = 3)
# 合并模块的颜色
mergedColors = merge$colors
table(mergedColors)
# 合并后新模块
mergedMEs = merge$newMEs

png("./plots/Fig4_dendrogram_and_module_colors_exp2.png")
plotDendroAndColors(geneTree,
                    cbind(dynamicColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE,
                    hang = 0.03,
                    addGuide = TRUE,
                    guideHang = 0.05)
dev.off()

# 对模块颜色重命名
moduleColors = mergedColors
# 将颜色转换为数据标签
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder)-1
MEs = mergedMEs
# 保存模块颜色和标签，以供后续部分使用
save(MEs, moduleLabels, moduleColors, geneTree, file = "./data/step6_networkConstruction_stepByStep_exp2.RData")

复现结果：

图4 D：

文章原图（左） vs CNSknowall R语言复现（右）

# Modules-traits relationships----
# Define numbers of genes and samples
nGenes = ncol(t_exp2)
nSamples = nrow(t_exp2)

# 这里使用merge过后的模块识别eigengenes

#使用分组信息
design=model.matrix(~0+ group2)
colnames(design)=levels(group2)
moduleTraitCor = cor(MEs, design, use = "p")
head(moduleTraitCor)

moduleTraitPvalue = corPvalueStudent(moduleTraitCor, 
                                     nSamples)

head(moduleTraitPvalue)

#可视化

# Will display correlations and their p-values
textMatrix =  paste(signif(moduleTraitCor, 2), 
                    "\n(",signif(moduleTraitPvalue, 1),
                    ")", 
                    sep = "");
dim(textMatrix) = dim(moduleTraitCor)
textMatrix[1:6,1:2]

design2 <- as.data.frame(design)
# Display the correlation values within a heatmap plot
png("./plots/Fig4_module_traits_relationship_exp2.png", width = 480, height = 720)
par(mar = c(5, 10, 4, 2));
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design2),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))
dev.off()


#2.Module membership----
# Define variable weight containing the weight column of datTrait
trait = as.data.frame(design2$pSS);
names(trait) = "trait"

# names (colors) of the modules
modNames = substring(names(MEs), 3)

geneModuleMembership = as.data.frame(cor(t_exp2, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")

geneTraitSignificance = as.data.frame(cor(t_exp2, trait, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))

names(geneTraitSignificance) = paste("GS.", names(trait), sep="")
names(GSPvalue) = paste("p.GS.", names(trait), sep="")

#选择与pSS性状相关性最高的模块
module = "darkorange"
column = match(module, modNames);
moduleGenes = mergedColors == module;

png("./plots/Fig4_GS_MM_exp2.png")
par(mar = c(4,4,4,4));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for body weight",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()

#提取该模块中的基因
pSS_module <- colnames(t_exp2)[mergedColors == "darkorange"]
pSS_module
save(pSS_module, file = "./data/step6_WGCNA_pSS_module.RData")

复现结果：

图4 F+H：

文章原图（左） vs CNSknowall R语言复现（右）

解析：与图4E+G一致，热图结果展现出一个与pSS性状高度正相关的模块（文献中是cyan模块，复现结果是darkorange模块），且复现结果的相关性与文献差距很小（复现0.73-文献0.76）。

    可以看出复现结果中的模块，其MM和GS的相关性更强（复现0.66>文献0.53）。

6.3 韦恩图取交集

图5 A：然后，我们利用维恩图得到了与PD和pss正相关的模块中34个重叠基因(图5A)。

复现代码：

library(VennDiagram)

png(filename = "./plots/Fig5_venn_WCGNA.png", width = 480, height = 480)
venn.plot <- venn.diagram(
  x = list(set1 = PD_module, set2 = pSS_module),
  category.names = c("PD_module", "pSS_module"),
  cat.cex = 1.5,    # 调整标签字体大小
  margin = 0.1,     # 调整页边距
  cex = 1.5,        # 调整区域内文本字体大小
  scaled = F,
  filename = NULL,
  output = TRUE,
  fill = c("#2CA02CFF", "#D62728FF"),
  cat.col = c("#2CA02CFF", "#D62728FF")
)
grid.draw(venn.plot)
dev.off()
common_WGCNA_genes <- intersect(PD_module, pSS_module)


load("./data/step3_deg.Rdata")
common_DEG_genes <- intersect(DEG1$name[DEG1$change != "stable"], DEG2$name[DEG2$change != "stable"])
png(filename = "./plots/Fig5_venn_W&D.png", width = 480, height = 480)
venn.plot <- venn.diagram(
  x = list(set1 = common_DEG_genes, set2 = common_WGCNA_genes),
  category.names = c("common_DEG_genes", "common_WGCNA_genes"),
  cat.cex = 1.5,    # 调整标签字体大小
  margin = 0.15,     # 调整页边距
  cex = 1.5,        # 调整区域内文本字体大小
  scaled = F,
  filename = NULL,
  output = TRUE,
  fill = c("#2CA02CFF", "#D62728FF"),
  cat.col = c("#2CA02CFF", "#D62728FF")
)
grid.draw(venn.plot)
dev.off()

#取交集后的基因用于机器学习
gene_for_ML <- intersect(common_DEG_genes, common_WGCNA_genes)
gene_for_ML

文章原图（左） vs CNSknowall R语言复现（右）

解析：复现结果找到的交集基因略多于原文献的结果。但是随后的生物学结论与原文献高度一致。

6.4 对交集基因进行富集分析

    进行了GO和KEGG途径富集分析，以揭示疾病相关常见基因的生物学功能和途径。共鉴定出82个GO术语，包括76个生物过程(BP)和6个分子功能(MF)，以及4个KEGG途径(详细结果见补充表S2)。

图5 B：在BP中，大多数基因主要参与免疫应答激活细胞表面受体信号通路、免疫应答调节细胞表面受体信号通路以及对脂多糖的应答(图5B)。

图5 C：在MF中，大多数基因主要参与光蛋白结合、磷酸酪氨酸残基结合和蛋白磷酸化氨基酸结合(图5C)。

图5 D：KEGG通路富集分析显示，大多数基因在白细胞跨内皮迁移、趋化因子信号通路和B细胞受体信号通路中高度增强(图SD)。

    这些结果表明，PD和psS主要与炎症和免疫反应有关，与常见DEGS的富集分析一致。

复现代码：

#使用AnnotationDbi包将symbol转成entrezid用于富集分析
library(org.Hs.eg.db)
library(AnnotationDbi)
keytypes(org.Hs.eg.db)
k <- AnnotationDbi::keys(org.Hs.eg.db, keytype = "ENTREZID") 
e2s <- AnnotationDbi::select(org.Hs.eg.db, keys = k,
                             columns="SYMBOL", keytype = "ENTREZID")
ENTREZID_commonWGCNA <- e2s[match(common_WGCNA_genes, e2s$SYMBOL),"ENTREZID"] %>% na.omit %>% as.numeric

library(clusterProfiler)
library(enrichplot)
#写一个循环，分别对go数据库的BP, MF, CC三个ontology做ORA富集分析
ont <- c("BP", "MF", "CC")
for(i in ont){
  go <- enrichGO(ENTREZID_commonWGCNA,
               OrgDb = "org.Hs.eg.db",
               ont = i,
               minGSSize = 0,
               pvalueCutoff = 0.05)
  p <- enrichplot::dotplot(go, showCategory = 10) +
    labs(title = i) +
    theme(plot.title = element_text(size = 20))
  ggsave(p,
         filename = paste0("./plots/Fig5_ORA_GO_", i, ".png"),
         width = 6,
         height = 6)
}

#对KEGG数据库进行ORA富集分析
kegg <- enrichKEGG(ENTREZID_commonWGCNA,
                   organism = "hsa",
                   pAdjustMethod = "none",
                   minGSSize = 0,
                   qvalueCutoff = 0.9)
p <- enrichplot::dotplot(kegg, showCategory = 12) +
    labs(title = "KEGG enrichment") +
    theme(plot.title = element_text(size = 20))
p
ggsave(p,
       filename = paste0("./plots/Fig5_ORA_KEGG.png"),
       width = 6, height = 6)

复现结果：

图5 B+C+D：

文章原图（左） vs CNSknowall R语言复现（右）

解析：由于前面找到的交集基因数量/种类有所不同，所以富集分析结果有一定差别，但是其主要通路是一致的。

7 机器学习

    为了更深入地研究PD和pss的共同发病机制，我们通过将173个DEGs和WGCNA筛选的34个基因取交集，鉴定了22个基因。随后，我们利用三种机器学习算法基于这22个共享基因筛选潜在的候选生物标志物。

7.1 PD数据集机器学习

图6 A-D：在PD数据集中，LASSO回归识别了7个基因(图6A)， SVM-RFE算法发现了10个均方根误差(RMSE)最低的基因(图6B)， 随机森林的重要性量表上的前10个基因(图6C)。

复现代码：

#筛选得到的基因提取表达矩阵
exprSet = exp1[gene_for_ML,]

#转置表达矩阵成行为样本
#x为转置后的矩阵
x = as.matrix(t(exprSet))
#y为生存或死亡事件
y = ifelse(group1 == "Healthy", 0, 1)

library(glmnet)
#调优参数
set.seed(21)
#cv.glmnet函数是glmnet包中的一个功能，主要用于进行交叉验证，以确定最优的正则化参数（lambda）用于广义线性模型（GLM）。
png("./plots/Fig6_lasso_1.png", width = 960, height = 480)
par(mfrow = c(1, 2))
#系数图
#glmnet 函数是R包glmnet中的一个核心函数，用于拟合带有L1（Lasso）或L2（Ridge）正则化的广义线性模型（GLM）。与cv.glmnet函数不同，glmnet不执行交叉验证，而是直接拟合模型，并为多个 λ 值提供解。
fit <- glmnet(x=x, y=y)
plot(fit,xvar = "lambda")

cv_fit <- cv.glmnet(x=x, y=y)
plot(cv_fit)
dev.off()

model_lasso_min <- glmnet(x=x, y=y, lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y, lambda=cv_fit$lambda.1se)

#提取model_lasso_1se构建模型所用的基因
lasso_1 <- model_lasso_1se$beta@Dimnames[[1]][model_lasso_1se$beta@i]
lasso_1

复现结果：

图6 A：

解析：Lasso回归选择正则化系数构建模型时：

若选择cv_fit$lambda.min（最小化均方误差的正则化系数），会得到均方误差（MSE）最小化的模型； 若选择cv_fit$lambda.1se（最小均方误差加上1个标准误差（SE）的正则化系数）。会得到最简单的模型，即在误差可接受范围内选择具有更强正则化（较大的alpha值）的一组回归系数。 相比最小误差标准，"1 SE" 标准会选择一个相对更大的正则化系数，使得模型更为稀疏（更多系数为零），从而减小方差，避免过拟合。

这里我们复现对两个正则化系数都构建了模型，1se模型筛选得到的基因：CXCR4, CSF2RB, LYN, HCLS1, KLHL6 （红色基因为文章鉴定到的三个biomarker）

复现代码：

library(e1071)
library(caret)

dat <- cbind(as.data.frame(t(exprSet)), group1)

# 设置种子
set.seed(1785)
# cv交叉验证次数10
control <- rfeControl(functions = caretFuncs, method = "cv", number = 10)
# 执行SVM-RFE算法
num <- ncol(dat)-1

if(file.exists("./data/step7_SVM-RFE_1.Rdata")) {
  load("./data/step7_SVM-RFE_1.Rdata")
  } else {
  results <- rfe(x = dat[, 1:num], #除去最后一列，其余列均为预测变量（也就是hubgene的表达量）
               y = dat$group1, # 分组信息
               sizes = c(1:num), 
               rfeControl = control,
               method = "svmRadial")
# save(results, file = "./data/step7_SVM-RFE_1.Rdata")
}


svmrfe_result <- data.frame(symbol = predictors(results)) 

plotdf <- data.frame(Variables = results$results$Variables, RMSE = sqrt(1 - results$results$Accuracy))
#RMSE最小点标注
min_y <- min(plotdf$RMSE)
min_x <- plotdf$Variables[which.min(plotdf$RMSE)]
#看图其实可以用N=4的点，因为此时的Accuracy有0.8968,比N=26时的0.8970相差不大
min_y <- 0.3213
min_x <- 4

# SVM-RFE结果简单可视化
p1 <- ggplot(data = plotdf, aes(x = Variables, y = RMSE))+
  geom_point(size = 2, fill = "#0097E0", color = "#0097E0")+
  geom_line(colour = "#0097E0")+
  theme_bw()+
  theme(plot.title = element_text(hjust =0.5,colour="black",face="bold",size=25),
        axis.text.x = element_text(size = 12), 
        axis.text.y = element_text(size = 12), 
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_blank(),
        legend.title = element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())+
  geom_point(data = data.frame(Variables = min_x, RMSE = min_y), 
             aes(x = Variables, y = RMSE), 
             color = "red", 
             size = 4, 
             shape = 21, 
             fill = "red") +  # 突出显示最低点
  geom_text(data = data.frame(Variables = min_x, RMSE = min_y), 
            aes(x = Variables, y = RMSE, label = paste0("N = ",min_x)), 
            vjust = -1.2, 
            hjust = 1.2, 
            color = "red", 
            size = 5)       # 标注 x 值
p1
ggsave("./plots/Fig6_svm_rfe_1.png", p1, width = 8, height = 6)

svm_rfe_1 <- head(results$optVariables, 4)
svm_rfe_1

复现结果：

图6 B：

文章原图（左） vs CNSknowall R语言复现（右）

解析：分类问题应该得到的是Accuracy，这里的文献中的RMSE不知道是怎么算得的。只能依葫芦画瓢画画看。我们就用Accuracy最高的点（这里使用次高的点，因为在保证accuracy仅次于最高点的同时使用了更少的特征）用到的基因，得到的4个基因为：TAGAP, CXCR4, CSF2RB, LYN（红色基因为文章鉴定到的三个biomarker）

复现代码：

library(randomForest)

x = t(exprSet)
y = group1
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
rf_output = randomForest(x=x, y=y, importance = TRUE, ntree = 500,
                         proximity=TRUE )

png("./plots/Fig6_RF_1.png", width = 1080, height = 480)
par(mfrow = c(1, 2))
plot(rf_output)
#top10的基因
varImpPlot(rf_output, type=2, n.var=10, scale=FALSE,
           main="MeanDecreaseGini",
           cex =2)
dev.off()

rf_importances = importance(rf_output, scale=FALSE)
head(rf_importances, 10)

#解释量top10的基因,和图上是一样的。
rf_1 = rev(rownames(tail(rf_importances[order(rf_importances[,4]),],10)))
rf_1

复现结果：

图6 C：

解析：选择了解释量前十的基因：CXCR4, TAGAP, LYN, MEF2C, CSF2RB, FCRLA, SRGN, HCLS1, FCRL5, ELTD1（红色基因为文章鉴定到的三个biomarker）

复现代码：

library(ggvenn)
ggvenn(list(LASSO = lasso_1, SVM_RFE = svm_rfe_1, Random_forest = rf_1),
  fill_color = c("green","red","lightblue"),   #填充圆的颜色，默认颜色会自动分配
  stroke_size = 0.5,   # 圆的边框线条粗细
  text_size = 4,
  set_name_size = 4,    # 集合名称的字体大小
  show_percentage = FALSE    #省略百分比
)
ggsave("./plots/Fig6_venn_exp1.png")

common_ML1 <- intersect(intersect(lasso_1, svm_rfe_1), rf_1)
common_ML1

复现结果：

图6 D：

文章原图（左） vs CNSknowall R语言复现（右）

解析：对PD数据集机器学习的三个结果取交集，文献最后找到的三个biomarker (CXCR4, LYN, CSF2RB) 全部都在三个集合的交集中。

7.2 pSS数据集机器学习

图6 E-H：同样，在pSS数据集中，LASSO回归发现了5个基因，SVM-RFE算法发现了8个基因，随机森林发现了10个基因(图6E-G)。在pss组中对鉴定到的基因取交集(图6H)。

复现代码：

#随机种子
set.seed(1)
#筛选得到的基因提取表达矩阵
exprSet = exp2[gene_for_ML,]

#转置表达矩阵成行为样本
#x为转置后的矩阵
x = as.matrix(t(exprSet))
#y为生存或死亡事件
y = ifelse(group2 == "Control", 0, 1)

library(glmnet)
#调优参数
#cv.glmnet函数是glmnet包中的一个功能，主要用于进行交叉验证，以确定最优的正则化参数（lambda）用于广义线性模型（GLM）。
png("./plots/Fig6_lasso_2.png", width = 960, height = 480)
par(mfrow = c(1, 2))
#系数图
#glmnet 函数是R包glmnet中的一个核心函数，用于拟合带有L1（Lasso）或L2（Ridge）正则化的广义线性模型（GLM）。与cv.glmnet函数不同，glmnet不执行交叉验证，而是直接拟合模型，并为多个 λ 值提供解。
fit <- glmnet(x=x, y=y)
plot(fit,xvar = "lambda")

cv_fit <- cv.glmnet(x=x, y=y ,nfolds = 10)
plot(cv_fit)
dev.off()

model_lasso_min <- glmnet(x=x, y=y, lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y, lambda=cv_fit$lambda.1se)

#提取model_lasso_1se构建模型所用的基因
lasso_2 <- model_lasso_min$beta@Dimnames[[1]][model_lasso_min$beta@i]
lasso_2

复现结果：

图6 E：

解析：得到用于构建1min模型的基因：CXCR4, CSF2RB, SRGN, FCRL5, RGS1（红色基因为文章鉴定到的2个biomarker，这里没找到LYN）

复现代码：

library(e1071)
library(caret)

dat <- cbind(as.data.frame(t(exprSet)), group2)

# 设置种子
set.seed(3)
# cv交叉验证次数10
control <- rfeControl(functions = caretFuncs, method = "cv", number = 10)
# 执行SVM-RFE算法
num <- ncol(dat)-1

if(file.exists("./data/step7_SVM-RFE_2.Rdata")) {
  load("./data/step7_SVM-RFE_2.Rdata")
  } else {
  results <- rfe(x = dat[, 1:num], #除去最后一列，其余列均为预测变量（也就是hubgene的表达量）
               y = dat$group2, # 分组信息
               sizes = c(1:num), 
               rfeControl = control,
               method = "svmRadial")
save(results, file = "./data/step7_SVM-RFE_2.Rdata")
}


svmrfe_result <- data.frame(symbol = predictors(results)) 

plotdf <- data.frame(Variables = results$results$Variables, RMSE = sqrt(1 - results$results$Accuracy))
#RMSE最小点标注
min_y <- min(plotdf$RMSE)
min_x <- plotdf$Variables[which.min(plotdf$RMSE)]
# SVM-RFE结果简单可视化
p2 <- ggplot(data = plotdf, aes(x = Variables, y = RMSE))+
  geom_point(size = 2, fill = "#0097E0", color = "#0097E0")+
  geom_line(colour = "#0097E0")+
  theme_bw()+
  theme(plot.title = element_text(hjust =0.5,colour="black",face="bold",size=25),
        axis.text.x = element_text(size = 12), 
        axis.text.y = element_text(size = 12), 
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_blank(),
        legend.title = element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())+
  geom_point(data = data.frame(Variables = min_x, RMSE = min_y), 
             aes(x = Variables, y = RMSE), 
             color = "red", 
             size = 4, 
             shape = 21, 
             fill = "red") +  # 突出显示最低点
  geom_text(data = data.frame(Variables = min_x, RMSE = min_y), 
            aes(x = Variables, y = RMSE, label = paste0("N = ",min_x)), 
            vjust = -1.2, 
            hjust = 1.2, 
            color = "red", 
            size = 5)       # 标注 x 值
p2
ggsave("./plots/Fig6_svm_rfe_2.png", p2, width = 8, height = 6)

svm_rfe_2 <- head(results$optVariables, 11)
svm_rfe_2

复现结果：

图6 F：

文章原图（左） vs CNSknowall R语言复现（右）

解析：我们使用Acurracy最高的点用到的基因，得到的4个基因为：CXCR4, CRLA, CSF2RB, MEF2C, FCGR2B, CTSS, LYN, EVI2B, PAG1, RCSD1, ZC3H12D（红色基因为文章鉴定到的三个biomarker）

复现代码：

library(randomForest)

x = t(exprSet)
y = group2
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
rf_output = randomForest(x=x, y=y, importance = TRUE, ntree = 500,
                         proximity=TRUE )

png("./plots/Fig6_RF_2.png", width = 1080, height = 480)
par(mfrow = c(1, 2))
plot(rf_output)
#top10的基因
varImpPlot(rf_output, type=2, n.var=10, scale=FALSE, 
           main="MeanDecreaseGini",
           cex = 2)
dev.off()

rf_importances = importance(rf_output, scale=FALSE)
head(rf_importances, 10)

#解释量top10的基因,和图上是一样的。
rf_2 = rev(rownames(tail(rf_importances[order(rf_importances[,4]),],10)))
rf_2

复现结果：

图6 G：

解析：选择了解释量前十的基因：CXCR4, FCRLA, CSF2RB, EVI2B, MEF2C, ELTD1, CTA-250D10.23, LYN, CTSS, PAG1（红色基因为文章鉴定到的三个biomarker）

复现代码：

ggvenn(list(LASSO = lasso_2, SVM_RFE = svm_rfe_2, Random_forest = rf_2),
  fill_color = c("green","red","lightblue"),   #填充圆的颜色，默认颜色会自动分配
  stroke_size = 0.5,   # 圆的边框线条粗细
  text_size = 4,
  set_name_size = 4,    # 集合名称的字体大小
  show_percentage = FALSE    #省略百分比
)
ggsave("./plots/Fig6_venn_exp2.png")

common_ML2 <- intersect(intersect(lasso_2, svm_rfe_2), rf_2)
common_ML2

复现结果：

图6 H：

文章原图（左） vs CNSknowall R语言复现（右）

ggvenn(list(LASSO = lasso_2, SVM_RFE = svm_rfe_2, Random_forest = rf_2),
  fill_color = c("green","red","lightblue"),   #填充圆的颜色，默认颜色会自动分配
  stroke_size = 0.5,   # 圆的边框线条粗细
  text_size = 4,
  set_name_size = 4,    # 集合名称的字体大小
  show_percentage = FALSE    #省略百分比
)
ggsave("./plots/Fig6_venn_exp2.png")

common_ML2 <- intersect(intersect(lasso_2, svm_rfe_2), rf_2)
common_ML2

解析：对pSS数据集机器学习的三个结果取交集，由于lasso回归没找到LYN，因此这里三种机器学习方法的交集基因只有CXCR4和CSF2RB

7.3 韦恩图取交集

图6 I：最后，通过维恩图确定CSF2RB、CXCR4和LYN基因为候选生物标志物(图6I)。

复现代码：

ggvenn(list(PD = common_ML1, pSS = common_ML2),
  fill_color = c("green","red"),   #填充圆的颜色，默认颜色会自动分配
  stroke_size = 0.5,   # 圆的边框线条粗细
  text_size = 4,
  set_name_size = 4,    # 集合名称的字体大小
  show_percentage = FALSE    #省略百分比
)
ggsave("./plots/Fig6_venn_biomarker.png")

复现结果：

图6 I：

文章原图（左） vs CNSknowall R语言复现（右）

解析：对两个数据集三种ML方法共同识别到的基因取交集

8 箱线图，列线图，ROC曲线    

    确定CSF2RB、CXCR4和LYN是否可以作为分析PD和psS的诊断标志物。

8.1 测试集（GSE16134，GSE40611）

图7 A+B：如图7A和B所示，上述三个基因在疾病缓解组的表达水平明显高于对照组。

图7 C+D：接下来，根据三个基因构建诺莫图，每个基因的相对表达量对应一个评分，将每个基因的得分相加计算总分(图7C和D)。

图7 E+F：ROC曲线分析显示，各生物标志物具有较好的预测性能，如下图所示(图7E和F):对于PD (GSE16134)数据集，CSF2RB (AUC: 0.914, 95% CI: 0.875-0.953)， CXCR4 (AUC: 0.917, 95% CI: 0.875-0.959)， LYN (AUC: 0.909, 95% CI 0.867-0.971)， Nomoscore(AUC: 0.935, 95% AUC: 0.899-0.971)。对于pSS (GSE40611)数据集 ，CSF2RB (AUC: 0.958, 95% CI: 0.897-1)， CXCR4 (AUC: 0.99, 95% CI: 0.97-1)， LYN (AUC: 0.954, 95% CI 0.891-1)，Nomoscore (AUC: 0.993, 95% CI 0.978-1)。

复现代码：

library(ggplot2)
library(ggpubr)
library(rstatix)
library(tidyverse)
library(ggdist)

#之前写的画箱线图的函数：
# 唯一参数：plotdf 用于可视化的数据框
#要求第一列是实验分组（factor）
#第二列是基因表达量（numeric）
#第三列是基因id（factor）

drawBoxplot <- function(plotdf, filename){
    if (!ncol(plotdf)==3) {
      stop("plotdf must be a dataframe with 3 columns.")
    }
    if (!class(plotdf$col1)=="factor"&class(plotdf$col3)=="factor") {
      stop("The first column and third column of plotdf must be factor.")
    }
    if (!class(plotdf$col2)=="numeric") {
      stop("The second column of plotdf must be numeric.")
    }
  color<-c("#1F77B4FF","#FF7F0EFF","#FFBB78FF","#2CA02CFF","#98DF8AFF","#D62728FF","#FF9896FF","#9467BDFF","#C5B0D5FF","#8C564BFF","#C49C94FF","#E377C2FF","#F7B6D2FF","#7F7F7FFF","#C7C7C7FF")
  
  #显著性
  stat.test <- plotdf %>%
    group_by(col3) %>%
    pairwise_wilcox_test(
      col2 ~ col1, paired = F,
      p.adjust.method = "bonferroni"
    )
  stat.test <- stat.test %>% add_xy_position(x = "col3",step.increase = 0.08)
  
  #散点图抖动宽度
  plotdf <- plotdf %>%
  group_by(col3) %>%
  mutate(jittered_x = as.numeric(col3) + 
                      (as.numeric(col1) - 1.5) * 0.4
                      )  # 自定义抖动宽度
  
  b<-ggplot(data=plotdf,aes(x = col3,y = col2))+
    geom_boxplot(mapping=aes(colour = col1 ), #箱线图
                 alpha = 0.5,
                 size = 0.7)+
     geom_jitter(aes(x = jittered_x, y = col2, color = col1),
                 height = 0,
                 width = 0.15,
                 alpha = 0.6) +
    # stat_dots(mapping=aes(colour = col1, fill = col1),
    #           position = position_dodge(0.75),
    #           scale = 1, alpha = 1,
    #           dotsize = 1,
    #           side = "left") +
    scale_fill_manual((limits=c(colnames(plotdf))),
                      values=color[1:length(plotdf$col1[!duplicated(plotdf$col1)])])+
    scale_color_manual(name = "group", (limits=c(colnames(plotdf))),
                      values=color[1:length(plotdf$col1[!duplicated(plotdf$col1)])])+
    theme_classic(base_line_size = 1)+
    theme(plot.title = element_text(size = 15,
                                    colour = "black",
                                    hjust = 0.5),
          axis.title.x = element_text(size = 15,
                                      color = "black",
                                      face = "bold",
                                      vjust = 1.9,
                                      hjust = 0.5,
                                      angle = 90),
          axis.title.y = element_text(size = 15,
                                      color = "black",
                                      face = "bold",
                                      vjust = 1.9,
                                      hjust = 0.5,
                                      angle = 90),
          legend.title = element_text(color="black", # 修改图例的标题
                                      size=15,
                                      face="bold"),
          legend.text = element_text(color="black", # 设置图例标签文字
                                     size = 10,
                                     face = "bold"),
          axis.text.x = element_text(size = 13,  # 修改X轴上字体大小，
                                     color = "black", # 颜色
                                     face = "bold", #  face取值：plain普通，bold加粗，italic斜体，bold.italic斜体加粗
                                     vjust = 0.5, # 位置
                                     hjust = 0.5,
                                     angle = 0), #角度
          axis.text.y = element_text(size = 13,
                                     color = "black",
                                     face = "bold",
                                     vjust = 0.5,
                                     hjust = 0.5,
                                     angle = 0)
    )+
    stat_pvalue_manual(
      stat.test, label = "p.adj.signif",
      linetype = 1,
      step.increase=0
    )+
    ylim(floor(min(plotdf$col2)),
         ceiling(max(plotdf$col2)))+
    labs(y = "Gene expression")
  ggsave(filename, b)
  return(b)
}


plotdf1 <- data.frame(col1 = factor(rep(group1, times = 3),
                                    levels = c("Healthy", "Diseased")),
                      col2 = c(t(exp1["CSF2RB",]),
                              t(exp1["CXCR4",]),
                              t(exp1["LYN",])),
                      col3 = factor(rep(c("CSF2RB", "CXCR4", "LYN"),
                                       each = length(group1))))

drawBoxplot(plotdf1, "./plots/Fig7_boxplot_exp1.png")

plotdf2 <- data.frame(col1 = factor(rep(group2, times = 3),
                                    levels = c("Control", "pSS")),
                      col2 = c(t(exp2["CSF2RB",]),
                              t(exp2["CXCR4",]),
                              t(exp2["LYN",])),
                      col3 = factor(rep(c("CSF2RB", "CXCR4", "LYN"),
                                       each = length(group2))))
drawBoxplot(plotdf2,
            "./plots/Fig7_boxplot_exp2.png")

复现结果：

图7 A+B：

解析：这里定义了一个绘制箱线图的函数drawBoxplot，是我常用的一个可视化风格。相较于对照组，可以看到这三个biomarker的表达量在PD组和pSS组都有显著升高。

复现代码：

library(Hmisc) ##加载Hmisc包与rms包
library(rms)

data = data.frame(group = 2 - as.numeric(group1),
                  CSF2RB = t(exp1["CSF2RB",]),
                  CXCR4 = t(exp1["CXCR4",]),
                  LYN = t(exp1["LYN",]))

dd <- datadist(data)
options(datadist = "dd")

#构建回归方程，使用lrm()函数构建二元LR，这里面的方程式就是结局~变量A+变量B+... 
f_lrm <- lrm(group ~ CSF2RB + CXCR4 + LYN, data = data)

#若报错：
# Error in reformulate(attr(termobj, "term.labels")[-dropx], response = if (keep.response) termobj[[2L]], : 'termlabels'必需是长度至少为一的字节矢量
#应该是rms包调用stats包的函数有问题，更新R或者给rms降版本可能可以解决。我是更新了整个R

#查看回归方程的结果 
summary(f_lrm)

nomogram <- nomogram(f_lrm,fun=function(x)1/(1+exp(-x)),   ##回归方程
                     fun.at = c(0.01,0.05,0.2,0.5,0.9,0.99), 
                     funlabel = "Risk",    ##风险轴刻度，结局是我的结局
                     conf.int = F,  #每个得分的置信度区间
                     abbrev = F)   #是否用简称代表因子变量
png(file = "./plots/Fig7_nomogram.png", width = 960, height = 480)
par(mgp=c(1.6,0.6,0), mar=c(5,5,3,1))
plot(nomogram, cex.axis=1.5)    #输出图片，右边plots可见列线图模型
dev.off()

Nomoscore <- predict(f_lrm, newdata = data[,2:4], type = "fitted")
head(Nomoscore)

#ROC曲线
roc_data = data.frame(group = 2 - as.numeric(group1),
                      CSF2RB = t(exp1["CSF2RB",]),
                      CXCR4 = t(exp1["CXCR4",]),
                      LYN = t(exp1["LYN",]),
                      Nomoscore = Nomoscore)
head(roc_data)

library(pROC)

#要求control组是0，所以前面的roc_data中group需要2 - group2
res <- roc(group ~ CSF2RB + CXCR4 + LYN + Nomoscore,
           data = roc_data,
           aur = TRUE,
           ci = TRUE, # 95%CI
           # percent=TRUE, ##是否需要以百分比显示
           direction="<" #设置分组方向
           )

library(ggplot2)
#CSF2RB
ci.auc(res$CSF2RB)
p1 <- ggroc(res$CSF2RB,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CSF2RB')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CSF2RB$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CSF2RB)[1],3), 
                          "-", 
                          round(ci.auc(res$CSF2RB)[3],3)))

#CXCR4
ci.auc(res$CXCR4)
p2 <- ggroc(res$CXCR4,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CXCR4')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CXCR4$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CXCR4)[1],3), 
                          "-", 
                          round(ci.auc(res$CXCR4)[3],3)))

#LYN
ci.auc(res$LYN)
p3 <- ggroc(res$LYN,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('LYN')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$LYN$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$LYN)[1],3), 
                          "-", 
                          round(ci.auc(res$LYN)[3],3)))

#Nomoscore
ci.auc(res$Nomoscore)
p4 <- ggroc(res$Nomoscore,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('Nomoscore')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$Nomoscore$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$Nomoscore)[1],3), 
                          "-", 
                          round(ci.auc(res$Nomoscore)[3],3)))

p <- cowplot::plot_grid(p1, p2, p3, p4, labels = "AUTO", nrow = 1)

ggsave("./plots/Fig7_ROC_exp1.png", p, width = 20, height = 5)

复现结果：

图7 C：

图7 E：

解析：诺莫图展示了各风险因素对PD临床结局事件的贡献大小，原文献中CXCR4贡献较大，CSF2RB和LYN贡献较小。复现图中三者贡献的得分比较相近。

    ROC曲线表明使用这三个基因单独预测临床结局或者使用诺莫图得分预测临床结局都具有较高的预测性能（AUC值均接近0.9）。且我们复现结果中四张图的AUC值均优于原文献结果。

复现代码：

#列线图
data = data.frame(group = 2 - as.numeric(group2),
                  CSF2RB = t(exp2["CSF2RB",]),
                  CXCR4 = t(exp2["CXCR4",]),
                  LYN = t(exp2["LYN",]))

dd <- datadist(data)
options(datadist = "dd")

#构建回归方程，使用lrm()函数构建二元LR，这里面的方程式就是结局~变量A+变量B+... 
f_lrm <- lrm(group ~ CSF2RB + CXCR4 + LYN, data = data)

#若报错：
# Error in reformulate(attr(termobj, "term.labels")[-dropx], response = if (keep.response) termobj[[2L]], : 'termlabels'必需是长度至少为一的字节矢量
#应该是rms包调用stats包的函数有问题，更新R或者给rms降版本可能可以解决。我是更新了整个R

#查看回归方程的结果 
summary(f_lrm)

nomogram <- nomogram(f_lrm,fun=function(x)1/(1+exp(-x)),   ##回归方程
                     fun.at = c(0.01,0.05,0.2,0.5,0.9,0.99), 
                     funlabel = "Risk",    ##风险轴刻度，结局是我的结局
                     conf.int = F,  #每个得分的置信度区间
                     abbrev = F)   #是否用简称代表因子变量
png(file = "./plots/Fig7_nomogram2.png", width = 960, height = 480)
par(mgp=c(1.6,0.6,0), mar=c(5,5,3,1))
plot(nomogram, cex.axis=2)    #输出图片，右边plots可见列线图模型
dev.off()


Nomoscore <- predict(f_lrm, newdata = data[,2:4], type = "fitted")
head(Nomoscore)

#ROC曲线
roc_data = data.frame(group = 2 - as.numeric(group2),
                      CSF2RB = t(exp2["CSF2RB",]),
                      CXCR4 = t(exp2["CXCR4",]),
                      LYN = t(exp2["LYN",]),
                      Nomoscore = Nomoscore)
head(roc_data)

library(pROC)

#要求control组是0，所以前面的roc_data中group需要2 - group2
res <- roc(group ~ CSF2RB + CXCR4 + LYN + Nomoscore,
           data = roc_data,
           aur = TRUE,
           ci = TRUE, # 95%CI
           # percent=TRUE, ##是否需要以百分比显示
           direction="<" #设置分组方向
           )

library(ggplot2)
#CSF2RB
ci.auc(res$CSF2RB)
p1 <- ggroc(res$CSF2RB,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CSF2RB')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CSF2RB$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CSF2RB)[1],3), 
                          "-", 
                          round(ci.auc(res$CSF2RB)[3],3)))

#CXCR4
ci.auc(res$CXCR4)
p2 <- ggroc(res$CXCR4,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CXCR4')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CXCR4$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CXCR4)[1],3), 
                          "-", 
                          round(ci.auc(res$CXCR4)[3],3)))

#LYN
ci.auc(res$LYN)
p3 <- ggroc(res$LYN,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('LYN')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$LYN$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$LYN)[1],3), 
                          "-", 
                          round(ci.auc(res$LYN)[3],3)))

#Nomoscore
ci.auc(res$Nomoscore)
p4 <- ggroc(res$Nomoscore,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('Nomoscore')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$Nomoscore$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$Nomoscore)[1],3), 
                          "-", 
                          round(ci.auc(res$Nomoscore)[3],3)))

p <- cowplot::plot_grid(p1, p2, p3, p4, labels = "AUTO", nrow = 1)

ggsave("./plots/Fig7_ROC_exp2.png", p, width = 20, height = 5)

复现结果：

图7 D+F：

解析：诺莫图展示了各风险因素对pSS临床结局事件的贡献大小，复现图与原文献CXCR4和和LYN4对结局贡献程度大小有差异，原因可能是上游进行数据清洗质控的步骤与原文献存在差异，因此复现结果中LYN表达量与原文献产生差异。

    ROC曲线表明使用这三个基因单独预测临床结局或者使用诺莫图得分预测临床结局都具有较高的预测性能（AUC值均大于0.9）。且我们复现结果中四张图使用三个基因单独预测临床结局的AUC值优于原文献结果，使用诺莫图得分预测临床结局的AUC值与文献结果基本一致。

9 使用验证数据集进行验证

    为了进一步评估候选生物标志物的准确性，我们在两个外部验证集(GSE10334用于PD和GSE23117用于pss)中对它们进行了验证。与既往结果一致，验证数据集的疾病组中3个biomarker的表达均显著上调(图8A和B)，且经诺莫图和ROC曲线分析证实，候选基因均具有良好的诊断价值(图8C-F)。

    综上所述，上述研究结果表明CSF2B、CXCR4和LYN可能作为PD和pss的有希望的诊断生物标志物。

数据下载：

dir.create("./validation")
dir.create("./validation/data")
dir.create("./validation/plots")
dir.create("./validation/GSE10334_RAW")
dir.create("./validation/GSE23117_RAW")

options(stringsAsFactors = F)

#读入GEO数据
library(GEOquery)
#列出文件夹下所有文件
dir()
# 牙周炎PD数据集
v1 <- getGEO(filename = "./validation/GSE10334_series_matrix.txt.gz", getGPL = F)
# 原发性干燥综合征pSS数据集
v2 <- getGEO(filename = "./validation/GSE23117_series_matrix.txt.gz", getGPL = F)

#建议一下进来就读取临床信息，在数据质控过滤过程中临床信息很重要
pd_v1 <- pData(v1)
pd_v2 <- pData(v2)

#表达矩阵(from series_matrix)
se_v1 <- exprs(v1)
se_v2 <- exprs(v2)
boxplot(se_v1)
boxplot(se_v2)
se_v2 <- log2(se_v2+1)


#表达矩阵(原始数据)
library(affy)
raw_v1 <- ReadAffy(celfile.path = "./validation/GSE10334_RAW/")
exp_v1 <- rma(raw_v1)
exp_v1 <- exprs(exp_v1)

raw_v2 <- ReadAffy(celfile.path = "./validation/GSE23117_RAW/")
exp_v2 <- rma(raw_v2)
exp_v2 <- exprs(exp_v2)

save(exp_v1, exp_v2, file = "./validation/data/Affy_raw.Rdata")

表达矩阵质控(同一基因对应多个探针取多行均值)和id转换

load("./validation/data/Affy_raw.Rdata")
range(exp_v1); dim(exp_v1)
range(exp_v2); dim(exp_v2)

#箱线图看一下各样本表达量分布方差齐性
boxplot(exp_v1)
boxplot(exp_v2)

# #pSS数据集做下批次效应矫正
# batch2 <- str_split(str_split(pd_v2$title, ",", simplify = T)[,4], " ", simplify = T)[,3]
# batch2[3:6] <- 2
# 
# library(sva)
# exp_v2 <- ComBat(exp_v2, batch = batch2)
# boxplot(exp_v2)

#*(选做)做一下归一化
# exp_v1 <- normalizeBetweenArrays(exp_v1)
# exp_v2 <- normalizeBetweenArrays(exp_v2)

#探针id转换成gene symbol
#芯片厂商是昂飞：Affymetrix Human Genome U133 Plus 2.0 Array，比较常见
#从网页https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570下载GPL570-55999.txt表格
gpl <- read.csv("GPL570-55999.txt", sep = "\t", header = T, comment.char = "#")
id2s <- data.frame(ID = gpl$ID, SYMBOL = gpl$Gene.Symbol)
id2s[1:4,]
#可以看到存在一个探针对应多个symbol的情况，有两种处理方法，一是保留第一个，二是删除所有一个探针对应多个基因的行。这里我们采取第二种方法
id2s <- id2s[!str_detect(id2s$SYMBOL, "///"), ]
#接下来转化表达矩阵行名，因为存在多个探针对应一个基因的情况，我们按照文献计算多个探针的均值

#写一个id转换的函数
IDtrans <- function(id2s, exp){
  
  #让表达矩阵与id2s的探针ID顺序一致
  exp <- exp[match(id2s$ID, rownames(exp)),]
  rownames(exp) <- id2s$SYMBOL
  exp <- as.data.frame(exp)
  
  exp$RowName <- id2s$SYMBOL
  
  #去除空列或者na列
  exp <- exp %>%
  filter(!is.na(RowName) & RowName != "")

  # 使用 dplyr 进行分组和计算平均值
  result <- exp %>%
    group_by(RowName) %>%
      summarise(across(everything(), mean))

  result <- as.data.frame(result)
  # 将行名设置为分组名称
  rownames(result) <- result$RowName
  result <- result[,-1]

  return(result)
}

#然后分别对exp1和exp2执行，注意，因为这两个表达矩阵的GPL平台都是GPL570，才能共用一套探针id转换规则（整理得到的id2s数据库）
exp_v1 <- IDtrans(id2s, exp_v1)
exp_v2 <- IDtrans(id2s, exp_v2)

dim(exp_v1); dim(exp_v2)

save(exp_v1, exp_v2, file = "./validation/data/transid_exp.RData")

整理临床信息顺序至与表达矩阵一致（原始数据）

#读取临床信息
pd_v1 <- pData(v1)
pd_v2 <- pData(v2)

colnames(exp_v1)
rownames(pd_v1)

#点号（.）在正则表达式中是一个特殊字符，表示任意单个字符，因此需要使用双反斜杠（\\.）进行转义。
colnames(exp_v1) <- str_split(colnames(exp_v1), "\\.", simplify = T)[,1]

colnames(exp_v2)
rownames(pd_v2)

colnames(exp_v2) <- str_split(colnames(exp_v2), "\\.", simplify = T)[,1]


#看下临床信息是否与表达矩阵一一对应
identical(rownames(pd_v1), colnames(exp_v1))
identical(rownames(pd_v2), colnames(exp_v2))
#第二个数据集false了


#提取其中的分组信息
table(pd_v1$`Gingival tissue phenotype:ch1`)
#Diseased (BoP + interproximal PPD of>3 mm + concomitant CAL of>2 mm)
#GPT翻译：（患病）：牙龈出血（BoP）+ 龈间袋探诊深度（PPD）> 3 毫米 + 同时存在的临床附着丧失（CAL）> 2 毫米
#Healthy (no BoP + PPD of<5 mm + CAL of<3 2 mm
#（健康）：无牙龈出血（BoP）+ 探诊深度（PPD）< 5 毫米 + 临床附着丧失（CAL）< 2 毫米
group_v1 <- ifelse(str_detect(pd_v1$`Gingival tissue phenotype:ch1`, "Diseased"), "PD", "Healthy")
group_v2 <- ifelse(str_detect(pd_v2$title, "non"), "Control", "pSS")

#最后将group类型改为factor，并且设定levels是实验组在前，对照组在后
table(group_v1)
group_v1 <- factor(group_v1, levels = c("PD", "Healthy"))
table(group_v2)
group_v2 <- factor(group_v2, levels = c("pSS", "Control"))

save(exp_v1, group_v1, pd_v1,
     exp_v2, group_v2, pd_v2,
     file = "./validation/data/step0.RData")

质控-PCA主成分分析 主要看下各分组的样本是否能聚在一起

#load("./validation/data/step0.RData")

qc_PCA <- function(exp, group){
  cg = names(tail(sort(apply(exp,1,sd)),1000))
  #apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
  
  ## 下面是画PCA的必须操作，需要看说明书。
  dat = t(exp[cg,])#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
  dat = as.data.frame(dat)#将matrix转换为data.frame 
  
  dat.pca <- PCA(dat , graph = FALSE)
  #现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  p <- fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind =  group, # color by groups 
               addEllipses = T,
               legend.title = "Groups"
  )
  return(p)
}

p1 <- qc_PCA(exp_v1, group_v1) + ggtitle("validation_PD_PCA")
p2 <- qc_PCA(exp_v2, group_v2) + ggtitle("validation_pSS_PCA")

#拼图
p1 + p2

ggsave('./validation/plots/QC_all_samples_PCA.png', p1 + p2, height = 5, width = 8)

箱线图绘制

library(ggplot2)
library(ggpubr)
library(rstatix)
library(tidyverse)
library(ggdist)

#之前写的画箱线图的函数：
# 唯一参数：plotdf 用于可视化的数据框
#要求第一列是实验分组（factor）
#第二列是基因表达量（numeric）
#第三列是基因id（factor）

drawBoxplot <- function(plotdf, filename){
    if (!ncol(plotdf)==3) {
      stop("plotdf must be a dataframe with 3 columns.")
    }
    if (!class(plotdf$col1)=="factor"&class(plotdf$col3)=="factor") {
      stop("The first column and third column of plotdf must be factor.")
    }
    if (!class(plotdf$col2)=="numeric") {
      stop("The second column of plotdf must be numeric.")
    }
  color<-c("#1F77B4FF","#FF7F0EFF","#FFBB78FF","#2CA02CFF","#98DF8AFF","#D62728FF","#FF9896FF","#9467BDFF","#C5B0D5FF","#8C564BFF","#C49C94FF","#E377C2FF","#F7B6D2FF","#7F7F7FFF","#C7C7C7FF")
  
  #显著性
  stat.test <- plotdf %>%
    group_by(col3) %>%
    pairwise_wilcox_test(
      col2 ~ col1, paired = F,
      p.adjust.method = "bonferroni"
    )
  stat.test <- stat.test %>% add_xy_position(x = "col3",step.increase = 0.08)
  
  #散点图抖动宽度
  plotdf <- plotdf %>%
  group_by(col3) %>%
  mutate(jittered_x = as.numeric(col3) + 
                      (as.numeric(col1) - 1.5) * 0.4
                      )  # 自定义抖动宽度
  
  b<-ggplot(data=plotdf,aes(x = col3,y = col2))+
    geom_boxplot(mapping=aes(colour = col1 ), #箱线图
                 alpha = 0.5,
                 size = 0.7)+
     geom_jitter(aes(x = jittered_x, y = col2, color = col1),
                 height = 0,
                 width = 0.15,
                 alpha = 0.6) +
    # stat_dots(mapping=aes(colour = col1, fill = col1),
    #           position = position_dodge(0.75),
    #           scale = 1, alpha = 1,
    #           dotsize = 1,
    #           side = "left") +
    scale_fill_manual((limits=c(colnames(plotdf))),
                      values=color[1:length(plotdf$col1[!duplicated(plotdf$col1)])])+
    scale_color_manual(name = "group", (limits=c(colnames(plotdf))),
                      values=color[1:length(plotdf$col1[!duplicated(plotdf$col1)])])+
    theme_classic(base_line_size = 1)+
    theme(plot.title = element_text(size = 15,
                                    colour = "black",
                                    hjust = 0.5),
          axis.title.x = element_text(size = 15,
                                      color = "black",
                                      face = "bold",
                                      vjust = 1.9,
                                      hjust = 0.5,
                                      angle = 90),
          axis.title.y = element_text(size = 15,
                                      color = "black",
                                      face = "bold",
                                      vjust = 1.9,
                                      hjust = 0.5,
                                      angle = 90),
          legend.title = element_text(color="black", # 修改图例的标题
                                      size=15,
                                      face="bold"),
          legend.text = element_text(color="black", # 设置图例标签文字
                                     size = 10,
                                     face = "bold"),
          axis.text.x = element_text(size = 13,  # 修改X轴上字体大小，
                                     color = "black", # 颜色
                                     face = "bold", #  face取值：plain普通，bold加粗，italic斜体，bold.italic斜体加粗
                                     vjust = 0.5, # 位置
                                     hjust = 0.5,
                                     angle = 0), #角度
          axis.text.y = element_text(size = 13,
                                     color = "black",
                                     face = "bold",
                                     vjust = 0.5,
                                     hjust = 0.5,
                                     angle = 0)
    )+
    stat_pvalue_manual(
      stat.test, label = "p.adj.signif",
      linetype = 1,
      step.increase=0
    )+
    ylim(floor(min(plotdf$col2)),
         ceiling(max(plotdf$col2)) + 1)+
    labs(y = "Gene expression")
  ggsave(filename, b)
  return(b)
}


plotdf1 <- data.frame(col1 = factor(rep(group_v1, times = 3),
                                    levels = c("Healthy", "PD")),
                      col2 = c(t(exp_v1["CSF2RB",]),
                              t(exp_v1["CXCR4",]),
                              t(exp_v1["LYN",])),
                      col3 = factor(rep(c("CSF2RB", "CXCR4", "LYN"),
                                       each = length(group_v1))))

drawBoxplot(plotdf1, "./plots/Fig8_boxplot_expV1.png")

plotdf2 <- data.frame(col1 = factor(rep(group_v2, times = 3),
                                    levels = c("Control", "pSS")),
                      col2 = c(t(exp_v2["CSF2RB",]),
                              t(exp_v2["CXCR4",]),
                              t(exp_v2["LYN",])),
                      col3 = factor(rep(c("CSF2RB", "CXCR4", "LYN"),
                                       each = length(group_v2))))
drawBoxplot(plotdf2,
            "./plots/Fig8_boxplot_expV2.png")

复现结果：图8 A+B：

解析：在验证集中，相较于对照组，可以看到三个biomarker的表达量在PD组都显著升高。在pSS组中CSF2RB和LYN表达量显著升高，CXCR4有表达升高趋势，但不显著，推测原因为验证集样本数量太少，组内方差太高。复现代码：

library(Hmisc) ##加载Hmisc包与rms包
library(rms)

data = data.frame(group = 2 - as.numeric(group_v1),
                  CSF2RB = t(exp_v1["CSF2RB",]),
                  CXCR4 = t(exp_v1["CXCR4",]),
                  LYN = t(exp_v1["LYN",]))

dd <- datadist(data)
options(datadist = "dd")

#构建回归方程，使用lrm()函数构建二元LR，这里面的方程式就是结局~变量A+变量B+... 
f_lrm <- lrm(group ~ CSF2RB + CXCR4 + LYN, data = data)

#若报错：
# Error in reformulate(attr(termobj, "term.labels")[-dropx], response = if (keep.response) termobj[[2L]], : 'termlabels'必需是长度至少为一的字节矢量
#应该是rms包调用stats包的函数有问题，更新R或者给rms降版本可能可以解决。我是更新了整个R

#查看回归方程的结果 
summary(f_lrm)

nomogram <- nomogram(f_lrm,fun=function(x)1/(1+exp(-x)),   ##回归方程
                     fun.at = c(0.01,0.05,0.2,0.5,0.9,0.99), 
                     funlabel = "Risk",    ##风险轴刻度，结局是我的结局
                     conf.int = F,  #每个得分的置信度区间
                     abbrev = F)   #是否用简称代表因子变量
png(file = "./plots/Fig8_nomogram1.png", width = 960, height = 480)
par(mgp=c(1.6,0.6,0), mar=c(5,5,3,1))
plot(nomogram, cex.axis=2)    #输出图片，右边plots可见列线图模型
dev.off()

Nomoscore <- predict(f_lrm, newdata = data[,2:4], type = "fitted")
head(Nomoscore)

#ROC曲线
roc_data = data.frame(group = 2 - as.numeric(group_v1),
                      CSF2RB = t(exp_v1["CSF2RB",]),
                      CXCR4 = t(exp_v1["CXCR4",]),
                      LYN = t(exp_v1["LYN",]),
                      Nomoscore = Nomoscore)
head(roc_data)

library(pROC)

#要求control组是0，所以前面的roc_data中group需要2 - group2
res <- roc(group ~ CSF2RB + CXCR4 + LYN + Nomoscore,
           data = roc_data,
           aur = TRUE,
           ci = TRUE, # 95%CI
           # percent=TRUE, ##是否需要以百分比显示
           direction="<" #设置分组方向
           )

library(ggplot2)
#CSF2RB
ci.auc(res$CSF2RB)
p1 <- ggroc(res$CSF2RB,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CSF2RB')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CSF2RB$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CSF2RB)[1],3), 
                          "-", 
                          round(ci.auc(res$CSF2RB)[3],3)))

#CXCR4
ci.auc(res$CXCR4)
p2 <- ggroc(res$CXCR4,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CXCR4')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CXCR4$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CXCR4)[1],3), 
                          "-", 
                          round(ci.auc(res$CXCR4)[3],3)))

#LYN
ci.auc(res$LYN)
p3 <- ggroc(res$LYN,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('LYN')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$LYN$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$LYN)[1],3), 
                          "-", 
                          round(ci.auc(res$LYN)[3],3)))

#Nomoscore
ci.auc(res$Nomoscore)
p4 <- ggroc(res$Nomoscore,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('Nomoscore')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$Nomoscore$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$Nomoscore)[1],3), 
                          "-", 
                          round(ci.auc(res$Nomoscore)[3],3)))

p <- cowplot::plot_grid(p1, p2, p3, p4, labels = "AUTO", nrow = 1)

ggsave("./plots/Fig8_ROC_exp1.png", p, width = 20, height = 5)

复现结果：图8 C+E：

解析：诺莫图展示了各风险因素对PD验证集临床结局事件的贡献大小，原文献中CXCR4贡献较大，CSF2RB和LYN贡献较小。复现图与原文献结果相近。

    ROC曲线表明使用这三个基因单独预测临床结局或者使用诺莫图得分预测临床结局都具有较高的预测性能（AUC值均接近0.9）。我们的复现结果与原文献较一致。

复现代码：

#列线图
data = data.frame(group = 2 - as.numeric(group_v2),
                  CSF2RB = t(exp_v2["CSF2RB",]),
                  CXCR4 = t(exp_v2["CXCR4",]),
                  LYN = t(exp_v2["LYN",]))

dd <- datadist(data)
options(datadist = "dd")

#构建回归方程，使用lrm()函数构建二元LR，这里面的方程式就是结局~变量A+变量B+... 
f_lrm <- lrm(group ~ CSF2RB + CXCR4 + LYN, data = data)

#若报错：
# Error in reformulate(attr(termobj, "term.labels")[-dropx], response = if (keep.response) termobj[[2L]], : 'termlabels'必需是长度至少为一的字节矢量
#应该是rms包调用stats包的函数有问题，更新R或者给rms降版本可能可以解决。我是更新了整个R

#查看回归方程的结果 
summary(f_lrm)

nomogram <- nomogram(f_lrm,fun=function(x)1/(1+exp(-x)),   ##回归方程
                     fun.at = c(0.01,0.05,0.2,0.5,0.9,0.99), 
                     funlabel = "Risk",    ##风险轴刻度，结局是我的结局
                     conf.int = F,  #每个得分的置信度区间
                     abbrev = F)   #是否用简称代表因子变量
png(file = "./plots/Fig8_nomogram2.png", width = 960, height = 480)
par(mgp=c(1.6,0.6,0), mar=c(5,5,3,1))
plot(nomogram, cex.axis=2)    #输出图片，右边plots可见列线图模型
dev.off()


Nomoscore <- predict(f_lrm, newdata = data[,2:4], type = "fitted")
head(Nomoscore)

#ROC曲线
roc_data = data.frame(group = 2 - as.numeric(group_v2),
                      CSF2RB = t(exp_v2["CSF2RB",]),
                      CXCR4 = t(exp_v2["CXCR4",]),
                      LYN = t(exp_v2["LYN",]),
                      Nomoscore = Nomoscore)
head(roc_data)

library(pROC)

#要求control组是0，所以前面的roc_data中group需要2 - group2
res <- roc(group ~ CSF2RB + CXCR4 + LYN + Nomoscore,
           data = roc_data,
           aur = TRUE,
           ci = TRUE, # 95%CI
           # percent=TRUE, ##是否需要以百分比显示
           direction="<" #设置分组方向
           )

library(ggplot2)
#CSF2RB
ci.auc(res$CSF2RB)
p1 <- ggroc(res$CSF2RB,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CSF2RB')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CSF2RB$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CSF2RB)[1],3), 
                          "-", 
                          round(ci.auc(res$CSF2RB)[3],3)))

#CXCR4
ci.auc(res$CXCR4)
p2 <- ggroc(res$CXCR4,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('CXCR4')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$CXCR4$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$CXCR4)[1],3), 
                          "-", 
                          round(ci.auc(res$CXCR4)[3],3)))

#LYN
ci.auc(res$LYN)
p3 <- ggroc(res$LYN,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('LYN')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$LYN$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$LYN)[1],3), 
                          "-", 
                          round(ci.auc(res$LYN)[3],3)))

#Nomoscore
ci.auc(res$Nomoscore)
p4 <- ggroc(res$Nomoscore,
            legacy.axes = TRUE, # 将X轴改为0-1，（默认是1-0）
            colour = "red"
            )+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
               color="darkgrey", linetype=4)+
  theme_bw()+# 设置背景
  ggtitle('Nomoscore')+
  theme(
    plot.title = element_text(
      size = 16,       # 设置标题大小
      face = "bold",   # 设置标题字体样式，例如加粗
      hjust = 0.5,     # 设置标题水平对齐，0是左对齐，0.5是居中，1是右对齐
      vjust = 1        # 设置标题垂直对齐，通常设置为1
    )
  )+
  annotate("text", x = 0.80, y = 0.20, size = 7,
           label = paste("AUC:", round(res$Nomoscore$auc,3)))+
  annotate("text", x = 0.80, y = 0.10, size = 7,
           label = paste0("CI: ", 
                          round(ci.auc(res$Nomoscore)[1],3), 
                          "-", 
                          round(ci.auc(res$Nomoscore)[3],3)))

p <- cowplot::plot_grid(p1, p2, p3, p4, labels = "AUTO", nrow = 1)

ggsave("./plots/Fig8_ROC_exp2.png", p, width = 20, height = 5)

复现结果：图8 D：图8 F：

解析：

    诺莫图展示了各风险因素对pSS验证集临床结局事件的贡献大小，原文献中LYB贡献较大，CXCR4贡献程度中等，CSF2RB贡献较小。复现图中CXCR4和CSF2RB的贡献程度与原文献相反，可能是上游质控产生的差异。

    ROC曲线表明使用这三个基因单独预测临床结局或者使用诺莫图得分预测临床结局都具有较高的预测性能（AUC值均大于0.9）。且我们复现结果中除CXCR4预测的AUC值偏低，其余AUC值与原文献比较接近。

10 免疫浸润

    鉴于PD和psS均为免疫反应较强的自身免疫性疾病，采用CIBERSORT分析免疫浸润细胞比例。

    鉴于PD和psS均为免疫反应较强的自身免疫性疾病，采用CIBERSORT分析免疫浸润细胞比例。

10.1 PD组免疫浸润

图9 A+C：图9A和B使用条形图展示了PD组每个样本中22种免疫细胞的分布情况。

    与对照样本相比，PD组织中plasma cells, activated memory CD4 + T cells, γδT cells, M0 macrophages, and neutrophils水平较高，而memory B cells, CD8 + T cells, follicular helper T cells, regulatory T cells (Tregs), activated NK cells, M1 macrophages, resting dendritic cells, activated dendritic cells, and resting mast cells （图9 C）水平较低。

复现代码：

堆叠柱状图：

library(CIBERSORT)
library(magrittr)
#读取包自带的LM22文件（免疫细胞特征基因文件）
sig_matrix <- system.file("extdata", "LM22.txt", package = "CIBERSORT")

data(LM22)
LM22[1:4,1:4]

result1 <- cibersort(sig_matrix = LM22, mixture_file = exp1, perm = 1000, QN = F)
save(result1, file = "./data/step10_cibersort1.Rdata")

#1. 堆叠柱状图
cibersort.res <- result1 %>% as.data.frame()   
  colnames(cibersort.res) %<>% gsub("_CIBERSORT","",.)
  colnames(cibersort.res) %<>% gsub("_"," ",.)
  cibersort.res <- cibersort.res[,-c(23:25)]  

res1 <- t(cibersort.res) %>% as.data.frame()
res1$cell_type <- rownames(res1)

#因为我们之前·1却人类group1和exp1是一一对应的才能这么写
group = data.frame(group = group1, id = colnames(exp1))

#定义颜色
fixed_colors <- colors_dutch <- c(
  '#EA2027','#EC3D25','#EE5A24','#F79F1F','#FFC312','#E1D324',
  '#C4E538','#A3CB38','#009432','#006266',
  '#12CBC4','#1289A7','#0652DD','#1B1464',
  '#5758BB','#9980FA','#D980FA','#FDA7DF',
  '#ED4C67','#B53471','#833471','#6F1E51'
)


png('./plots/Fig9_ICI_barplot_exp1.png',w=12,h=8,units='in',res=600, family = 'Times')
res1 %>%
  gather(id, fraction, -cell_type) %>%
  merge(group,by='id')%>%
  # 绘制堆积条形图
  ggplot(aes(x=id, y=fraction, fill=cell_type)) +
  geom_bar(position = 'stack',stat = 'identity')+
  scale_y_continuous(expand = c(0,0))+
  theme_bw()+
  labs(x='',
       y='Relative Percent',
       fill='')+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = 'top',
        text = element_text(size=15)) +
  scale_fill_manual(values = fixed_colors)+
  facet_grid(~group,scales= "free",space= "free")
dev.off()

箱线图：

#2. 箱线图
# Wilcoxon秩和检验比较两组（PD vs对照样本）免疫细胞浸润的差异（p< 0.05） --------------------------------------------------------------------
tiics_result <- t(cibersort.res) %>% as.matrix()
rownames(group) <- group$id
colnames(group) <- c("group", "sample") 
#新建三个空的矩阵
pvalue = padj = log2FoldChange <- matrix(0, nrow(tiics_result), 1)

group_High <- group$sample[which(group$group == 'Diseased')]
# group_High<-as.character(group_High$sample)
group_Low <- group$sample[which(group$group == 'Healthy')]
# group_Low<-as.character(group_Low$sample)

#wilcox.test计算p值，以及log2FC
for (i in 1:nrow(tiics_result)){
  # print(tiics_result[i][group_High])
  pvalue[i, 1] = p.value = 
    wilcox.test(tiics_result[i, group_High], tiics_result[i, group_Low])$p.value
  
  log2FoldChange[i, 1] = mean(tiics_result[i, group_High]) - 
    mean(tiics_result[i, group_Low])
}

padj <- p.adjust(as.vector(pvalue), "fdr", n = length(pvalue))
rTable <- data.frame(log2FoldChange, 
                     pvalue, 
                     padj,
                     row.names = rownames(tiics_result))

control <- signif(apply(tiics_result[rownames(rTable), group_Low], 
                        1,
                        mean), 4)
case <- signif(apply(tiics_result[rownames(rTable), group_High], 
                     1, 
                     mean), 4)

rTable <- data.frame(control, 
                     case,
                     rTable[, c("padj", "pvalue", "log2FoldChange")])


rTable$immune_cell <- rownames(rTable)
rTable$sig <- ifelse(rTable$pvalue < 0.05,
                     ifelse(rTable$pvalue < 0.01, 
                            ifelse(rTable$pvalue < 0.001,
                                   ifelse(rTable$pvalue < 0.0001,
                                          paste(rTable$immune_cell, "****",  sep = ""),
                                          paste(rTable$immune_cell, "***", sep = "")),
                                   paste(rTable$immune_cell, "**", sep = "")),
                            paste(rTable$immune_cell, "*",  sep = "")), 
                     rTable$immune_cell)  

write.csv(rTable,
          file = "./data/tiics_wilcox_test_1.csv",
          row.names = F)

# diff_Table <- rTable[which(rTable$pvalue<0.05),]
#diff_Table <- rTable
# write.csv(diff_Table,
#           file = "diff_tiics_wilcox_test.csv",
#           row.names = F)


library(tidyr)
library(ggplot2)
library(ggpubr)
#devtools::install_github("kongdd/Ipaper")
library(Ipaper)

xCell2 <- data.frame(Immune_Cell = rownames(tiics_result), 
                     tiics_result, 
                     pvalue=rTable$pvalue)
# xCell2 <- xCell2[which(xCell2$pvalue<0.05),]

diff_tiics <- rownames(xCell2)
violin_dat <- gather(xCell2,
                     key = Group, #新的列名：变量名称
                     value = score, #新的列名：变量值
                     -c("Immune_Cell","pvalue")) #要聚合的列，除了Immune_Cell和pvalue

violin_dat$Group <- ifelse(gsub("\\.","-",violin_dat$Group) %in% group_Low,
                           "Control", "PD") 
head(violin_dat)
boxplot_diff_TIICs <- ggplot(violin_dat, aes(x=Immune_Cell, 
                                             y=score,
                                             fill=Group)) +
  geom_violin(trim=T,color="black") + #绘制小提琴图, “color=”设置小提琴图的轮廓线的颜色(#不要轮廓可以设为white以下设为背景为白色，其实表示不要轮廓线)#"trim"如果为TRUE(默认值),则将小提琴的尾部修剪到数据范围。如果为FALSE,不修剪尾部。
  stat_boxplot(geom="errorbar",
               width=0.1,
               position = position_dodge(0.9)) +
  geom_boxplot(width=0.7,
               position=position_dodge(0.9),
               outlier.shape = NA)+ #绘制箱线图，此处width=0.1控制小提琴图中箱线图的宽窄
  geom_point(aes(fill = Group),
             size = 0.05,
             position = position_dodge(0.9))+
  scale_fill_manual(values= c("#8DD3C7","#E31A1C"))+ #设置填充的颜色
  labs(title="", x="", y = "Score",size=20) +
  stat_compare_means(data = violin_dat,
                     mapping = aes(group = Group),
                     label ="p.signif",
                     hide.ns = F) +
  theme_bw()+#把背景设置为白底
  theme(plot.title = element_text(hjust =0.5,colour="black",face="bold",size=18), # 将图表标题居中
        axis.text.x=element_text(angle=45,hjust=1,colour="black",face="bold",size=10), #设置x轴刻度标签的字体显示倾斜角度为45度，并向下调整1(hjust = 1)，字体大小为14
        axis.text.y=element_text(hjust=0.5,colour="black",size=12), #设置y轴刻度标签的字体簇，字体大小，字体样式为plain
        axis.title.x=element_text(size=16,face="bold"),#设置x轴标题的字体属性
        axis.title.y=element_text(size=14,face="bold"), #设置y轴标题的字体属性
        legend.text=element_text(face="bold", hjust = 0.5,colour="black", size=11), #设置图例的子标题的字体属性
        legend.title=element_text(face="bold", colour="black", size=11),#设置图例的总标题的字体属性
        #legend.justification=c(-0.1,1.2), #可调整图例的位置。##(1,1)第一个1是调整图例在图的内外(左右移动)，第二个1是在图中上下移动。
        #legend.position=c(0, 1.04), #legend.position=c(0,1)左上角，(1,1)是在右上角。
        panel.grid.major = element_blank(), #不显示网格线
        panel.grid.minor = element_blank()) #不显示网格线
boxplot_diff_TIICs
ggsave('./plots/Fig9_boxplot_exp1.png', boxplot_diff_TIICs, w=10, h=6)

复现结果：

图9 A：

文章原图（上） vs CNSknowall R语言复现（下）

图9 C：

文章原图（上） vs CNSknowall R语言复现（下）

解析：复现结果中能观察到plasma cells, activated memory CD4 + T cells, γδT cells, M0 macrophages, and neutrophils 水平升高（红框框出），而memory B cells, CD8 + T cells, follicular helper T cells, regulatory T cells (Tregs), activated NK cells, M1 macrophages, resting dendritic cells, activated dendritic cells, and resting mast cells 水平降低（蓝框框出）。表明复现结果中免疫细胞的变化趋势与原文献结果几乎一致。

图9 E：随后，PD样本免疫细胞之间的相关性热图显示，activated NK cells与

resting mast cells呈正相关(r = 0.58)，而plasma cells与resting dendritic cells负相关(r = -0.64)(图9E)。

图9 F：见图9 F+H

复现代码：

免疫细胞间相关性热图：

#3 相关性矩阵
library(corrplot)

#创建空矩阵
cor_mat <- matrix(NA, nrow = 22, ncol = 22)
rownames(cor_mat) <- rownames(tiics_result)
colnames(cor_mat) <- rownames(tiics_result)
#计算相关性
for(i in 1:nrow(tiics_result)){
  for(j in 1:nrow(tiics_result)){
    cor_mat[i,j] <- stats::cor(tiics_result[i,],tiics_result[j,])
  }
}

png("./plots/Fig9_cor_matrix_1.png", width = 1200, height = 1200)

#要先画下面的一半
corrplot(cor_mat,
         mar = c(18, 6, 6, 12),
         method = c('number'), 
         type = c('lower'), 
         col = colorRampPalette(c("#00008B", "#555597", "#FFFFB0", 
                                           "#B1553A", "#8B0000"))(100),
         order = c('AOE'), 
         diag = TRUE, 
         number.cex = 1,
         tl.pos = 'tl', #不显示文本标签
         tl.col = 'black', #横纵坐标文本文字颜色
         tl.cex = 1, #横纵坐标文本文字大小
         tl.srt = 45 #横纵坐标文本旋转角度
)

corrplot(cor_mat, add = TRUE,
         method = c('pie'), 
         type = c('upper'), 
         col = colorRampPalette(c("#00008B", "#555597", "#FFFFB0", 
                                           "#B1553A", "#8B0000"))(100),
         outline = 'grey', 
         order = c('AOE'), 
         diag = FALSE,
         tl.pos = 'n' #不显示文本标签
)

dev.off()

免疫与三个biomarker基因间相关性热图：

#4. 热图
### 关键基因与差异细胞相关性----------------------------------------------------
library(pheatmap)

dat_allcell <- t(cibersort.res) %>% as.matrix()
dim(dat_allcell)

dat_diffcell <- dat_allcell %>% t %>% as.data.frame 

hubgene <- c("CSF2RB", "CXCR4", "LYN")
dat_hubgene <- exp1[hubgene, group$sample] %>% t %>% as.data.frame 


cor_r <- cor(dat_diffcell, dat_hubgene,method = "spearman") 
cor_p <- WGCNA::corPvalueStudent(cor_r,length(rownames(dat_hubgene)))


cor_r2 <- cor_r %>% as.data.frame %>% tibble::rownames_to_column(var = "cell") %>%
  tidyr::gather(., gene,Correlation,-cell) #转换数据长短 cor 用于gglot2可视化
cor_p2 <- cor_p %>% as.data.frame %>% tibble::rownames_to_column(var = "cell") %>% 
  tidyr::gather(., gene, Pvalue, -cell) #转换数据长短 p 用于ggplot2可视化


cor_dat <- cbind(cor_r2, cor_p2)[,c("cell", "gene", "Correlation", "Pvalue")]
colnames(cor_dat) <- c('cell','gene','correlation','p')
cor_dat$cell <- gsub('.',' ',cor_dat$cell, fix=T)
write.csv(cor_dat, "./data/correlation_cor_exp1.csv")

cor_dat$sig <- cut(cor_dat$p,
                   breaks = c(0, 0.001, 0.01, 0.05, 1),
                   labels = c("***","**","*",""),
                   right=FALSE)
#可视化

p <- ggplot(cor_dat,aes(x = gene, y = cell))+
  geom_tile(aes(fill = correlation), size = 0.5)+
  geom_text(aes(label = sig), color = "black", size = 3) +
  scale_x_discrete(expand=c(0,0)) +
  scale_y_discrete(expand=c(0,0),position="left")+
  theme(text=element_text(family="Roboto"),
        axis.text.x=element_text(color="black"),
        axis.text.y=element_text(color="black"),
        axis.ticks.x = element_blank(),
        axis.ticks.y=element_blank(),
        panel.border = element_rect(fill=NA,color="grey80",
        size=1, linetype="solid"))+
  scale_fill_gradient2(low = "#00008B", high = "#8B0000")

ggsave("./plots/Fig9_heatmap_exp1.png", p, height = 8, width = 5)

复现结果：

图9 E：

文章原图（上） vs CNSknowall R语言复现（下）

图9 F：

解析：图9E显示activated NK cells与resting mast cells呈正相关(原文献r = 0.58，复现结果r = 0.35存在一定差异，可能是上游质控处理有差异导致)，而plasma cells与resting dendritic cells负相关(原文献r = -0.64，复现结果r = -0.64)。图9F与图9H的解析放在一起，在推文后面。

10.2 pSS组免疫浸润

图9 B+D：图9A和B使用条形图展示了pSS组中每个样本中22种免疫细胞的分布情况。

    pSS组中初始naïve B cells, memory B cells, naïve CD4 + T cells, activated memory CD4+ T cells, γδT cells, M0macrophages, and M1 macrophages 水平升高，plasma cells, CD8 + T cells, Tregs, resting NK cells, monocytes, activated dendritic cells, 和 activated mast cells 水平降低（图9 D）。

    有趣的是，在PD和pSS样本中部分细胞水平的变化展现出类似的趋势，activated memory CD4 + T cells, γδT cells, M0 macrophages 增加；CD8 + T cells, Tregs, activated dendritic cells减少。

复现代码：

堆叠柱状图：

library(CIBERSORT)
library(magrittr)
#读取包自带的LM22文件（免疫细胞特征基因文件）
sig_matrix <- system.file("extdata", "LM22.txt", package = "CIBERSORT")

data(LM22)
LM22[1:4,1:4]

result2 <- cibersort(sig_matrix = LM22, mixture_file = exp2, perm = 1000, QN = F)
save(result2, file = "./data/step10_cibersort2.Rdata")

#1. 堆叠柱状图
cibersort.res <- result2 %>% as.data.frame()   
  colnames(cibersort.res) %<>% gsub("_CIBERSORT","",.)
  colnames(cibersort.res) %<>% gsub("_"," ",.)
  cibersort.res <- cibersort.res[,-c(23:25)]  

res2 <- t(cibersort.res) %>% as.data.frame()
res2$cell_type <- rownames(res2)

#因为我们之前·1却人类group1和exp1是一一对应的才能这么写
group = data.frame(group = group2, id = colnames(exp2))

#定义颜色
fixed_colors <- colors_dutch <- c(
  '#EA2027','#EC3D25','#EE5A24','#F79F1F','#FFC312','#E1D324',
  '#C4E538','#A3CB38','#009432','#006266',
  '#12CBC4','#1289A7','#0652DD','#1B1464',
  '#5758BB','#9980FA','#D980FA','#FDA7DF',
  '#ED4C67','#B53471','#833471','#6F1E51'
)


png('./plots/Fig9_ICI_barplot_exp2.png',w=12,h=8,units='in',res=600, family = 'Times')
res2 %>%
  gather(id, fraction, -cell_type) %>%
  merge(group,by='id')%>%
  # 绘制堆积条形图
  ggplot(aes(x=id, y=fraction, fill=cell_type)) +
  geom_bar(position = 'stack',stat = 'identity')+
  scale_y_continuous(expand = c(0,0))+
  theme_bw()+
  labs(x='',
       y='Relative Percent',
       fill='')+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = 'top',
        text = element_text(size=15)) +
  scale_fill_manual(values = fixed_colors)+
  facet_grid(~group,scales= "free",space= "free")
dev.off() 

箱线图：

#2. 箱线图
# Wilcoxon秩和检验比较两组（PD vs对照样本）免疫细胞浸润的差异（p< 0.05） --------------------------------------------------------------------
tiics_result <- t(cibersort.res) %>% as.matrix()
rownames(group) <- group$id
colnames(group) <- c("group", "sample") 
#新建三个空的矩阵
pvalue = padj = log2FoldChange <- matrix(0, nrow(tiics_result), 1)

group_High <- group$sample[which(group$group == 'pSS')]
# group_High<-as.character(group_High$sample)
group_Low <- group$sample[which(group$group == 'Control')]
# group_Low<-as.character(group_Low$sample)

#wilcox.test计算p值，以及log2FC
for (i in 1:nrow(tiics_result)){
  # print(tiics_result[i][group_High])
  pvalue[i, 1] = p.value = 
    wilcox.test(tiics_result[i, group_High], tiics_result[i, group_Low])$p.value
  
  log2FoldChange[i, 1] = mean(tiics_result[i, group_High]) - 
    mean(tiics_result[i, group_Low])
}

padj <- p.adjust(as.vector(pvalue), "fdr", n = length(pvalue))
rTable <- data.frame(log2FoldChange, 
                     pvalue, 
                     padj,
                     row.names = rownames(tiics_result))

control <- signif(apply(tiics_result[rownames(rTable), group_Low], 
                        1,
                        mean), 4)
case <- signif(apply(tiics_result[rownames(rTable), group_High], 
                     1, 
                     mean), 4)

rTable <- data.frame(control, 
                     case,
                     rTable[, c("padj", "pvalue", "log2FoldChange")])


rTable$immune_cell <- rownames(rTable)
rTable$sig <- ifelse(rTable$pvalue < 0.05,
                     ifelse(rTable$pvalue < 0.01, 
                            ifelse(rTable$pvalue < 0.001,
                                   ifelse(rTable$pvalue < 0.0001,
                                          paste(rTable$immune_cell, "****",  sep = ""),
                                          paste(rTable$immune_cell, "***", sep = "")),
                                   paste(rTable$immune_cell, "**", sep = "")),
                            paste(rTable$immune_cell, "*",  sep = "")), 
                     rTable$immune_cell)  

write.csv(rTable,
          file = "./data/tiics_wilcox_test_2.csv",
          row.names = F)

# diff_Table <- rTable[which(rTable$pvalue<0.05),]
#diff_Table <- rTable
# write.csv(diff_Table,
#           file = "diff_tiics_wilcox_test.csv",
#           row.names = F)


library(tidyr)
library(ggplot2)
library(ggpubr)
#devtools::install_github("kongdd/Ipaper")
library(Ipaper)

xCell2 <- data.frame(Immune_Cell = rownames(tiics_result), 
                     tiics_result, 
                     pvalue=rTable$pvalue)
# xCell2 <- xCell2[which(xCell2$pvalue<0.05),]

diff_tiics <- rownames(xCell2)
violin_dat <- gather(xCell2,
                     key = Group, #新的列名：变量名称
                     value = score, #新的列名：变量值
                     -c("Immune_Cell","pvalue")) #要聚合的列，除了Immune_Cell和pvalue

violin_dat$Group <- ifelse(gsub("\\.","-",violin_dat$Group) %in% group_Low,
                           "Control", "pSS") 
head(violin_dat)
boxplot_diff_TIICs <- ggplot(violin_dat, aes(x=Immune_Cell, 
                                             y=score,
                                             fill=Group)) +
  geom_violin(trim=T,color="black") + #绘制小提琴图, “color=”设置小提琴图的轮廓线的颜色(#不要轮廓可以设为white以下设为背景为白色，其实表示不要轮廓线)#"trim"如果为TRUE(默认值),则将小提琴的尾部修剪到数据范围。如果为FALSE,不修剪尾部。
  stat_boxplot(geom="errorbar",
               width=0.1,
               position = position_dodge(0.9)) +
  geom_boxplot(width=0.7,
               position=position_dodge(0.9),
               outlier.shape = NA)+ #绘制箱线图，此处width=0.1控制小提琴图中箱线图的宽窄
  geom_point(aes(fill = Group),
             size = 0.05,
             position = position_dodge(0.9))+
  scale_fill_manual(values= c("#8DD3C7","#E31A1C"))+ #设置填充的颜色
  labs(title="", x="", y = "Score",size=20) +
  stat_compare_means(data = violin_dat,
                     mapping = aes(group = Group),
                     label ="p.signif",
                     hide.ns = F) +
  theme_bw()+#把背景设置为白底
  theme(plot.title = element_text(hjust =0.5,colour="black",face="bold",size=18), # 将图表标题居中
        axis.text.x=element_text(angle=45,hjust=1,colour="black",face="bold",size=10), #设置x轴刻度标签的字体显示倾斜角度为45度，并向下调整1(hjust = 1)，字体大小为14
        axis.text.y=element_text(hjust=0.5,colour="black",size=12), #设置y轴刻度标签的字体簇，字体大小，字体样式为plain
        axis.title.x=element_text(size=16,face="bold"),#设置x轴标题的字体属性
        axis.title.y=element_text(size=14,face="bold"), #设置y轴标题的字体属性
        legend.text=element_text(face="bold", hjust = 0.5,colour="black", size=11), #设置图例的子标题的字体属性
        legend.title=element_text(face="bold", colour="black", size=11),#设置图例的总标题的字体属性
        #legend.justification=c(-0.1,1.2), #可调整图例的位置。##(1,1)第一个1是调整图例在图的内外(左右移动)，第二个1是在图中上下移动。
        #legend.position=c(0, 1.04), #legend.position=c(0,1)左上角，(1,1)是在右上角。
        panel.grid.major = element_blank(), #不显示网格线
        panel.grid.minor = element_blank()) #不显示网格线
boxplot_diff_TIICs
ggsave('./plots/Fig9_boxplot_exp2.png', boxplot_diff_TIICs, w=10, h=6)

复现结果：

图9 B：

文章原图（上） vs CNSknowall R语言复现（下）

图9 D：

解析：pSS组中初始naïve B cells, memory B cells, naïve CD4 + T cells, activated memory CD4+ T cells, γδT cells, M0macrophages, and M1 macrophages 水平升高，plasma cells, CD8 + T cells, Tregs, resting NK cells, monocytes, activated dendritic cells, 和 activated mast cells 水平降低。

复现结果大部分免疫细胞的变化趋势与原文献一致，memory B cells表现出水平升高趋势但未到达显著，activated dendritic cells和monocytes水平降低但未达到显著，推测可能是上游质控处理有差异，或者样本量太少导致。

图9 G：在pss样本中，CD8 + T cells与Tregs (r = 0.83)和单核monocytes(r = 0.82)呈显著正相关。相比之下，plasma cells 与 γδT cells 呈负相关(r = -0.63)(图9F)。这些发现表明，与对照组相比，患者存在不同的免疫模式，以及各种类型免疫细胞之间的相互作用。

图9 F+H：共享生物标志物与免疫细胞之间的相关性分析也进行了研究。结果显示在PD和pSS样本中，γδT cells, activated memory CD4 + T cells, M0 Macrophages 和 naïve B cells 与 CSR2RB、CXCR4和LYN基因表达呈正相关（图9 F和H）。在验证数据集中也观察到了类似的结果，进一步证实了这些发现（补充图1A和B）。这表明这三种基因可能在PD并发pSS的发病机制中发挥关键作用，影响免疫细胞的浸润。

#3 相关性矩阵
library(corrplot)

#创建空矩阵
cor_mat <- matrix(NA, nrow = 22, ncol = 22)
rownames(cor_mat) <- rownames(tiics_result)
colnames(cor_mat) <- rownames(tiics_result)
#计算相关性
for(i in 1:nrow(tiics_result)){
  for(j in 1:nrow(tiics_result)){
    cor_mat[i,j] <- stats::cor(tiics_result[i,],tiics_result[j,])
  }
}

png("./plots/Fig9_cor_matrix_2.png", width = 1200, height = 1200)

#要先画下面的一半
corrplot(cor_mat,
         mar = c(18, 6, 6, 12),
         method = c('number'), 
         type = c('lower'), 
         col = colorRampPalette(c("#00008B", "#555597", "#FFFFB0", 
                                           "#B1553A", "#8B0000"))(100),
         order = c('AOE'), 
         diag = TRUE, 
         number.cex = 1,
         tl.pos = 'tl', #不显示文本标签
         tl.col = 'black', #横纵坐标文本文字颜色
         tl.cex = 1, #横纵坐标文本文字大小
         tl.srt = 45 #横纵坐标文本旋转角度
)

corrplot(cor_mat, add = TRUE,
         method = c('pie'), 
         type = c('upper'), 
         col = colorRampPalette(c("#00008B", "#555597", "#FFFFB0", 
                                           "#B1553A", "#8B0000"))(100),
         outline = 'grey', 
         order = c('AOE'), 
         diag = FALSE,
         tl.pos = 'n' #不显示文本标签
)

dev.off()

免疫与三个biomarker基因间相关性热图：

#4. 热图
### 关键基因与差异细胞相关性----------------------------------------------------
library(pheatmap)

dat_allcell <- t(cibersort.res) %>% as.matrix()
dim(dat_allcell)

dat_diffcell <- dat_allcell %>% t %>% as.data.frame 

hubgene <- c("CSF2RB", "CXCR4", "LYN")
dat_hubgene <- exp2[hubgene, group$sample] %>% t %>% as.data.frame 


cor_r <- cor(dat_diffcell, dat_hubgene,method = "spearman") 
cor_p <- WGCNA::corPvalueStudent(cor_r,length(rownames(dat_hubgene)))


cor_r2 <- cor_r %>% as.data.frame %>% tibble::rownames_to_column(var = "cell") %>%
  tidyr::gather(., gene,Correlation,-cell) #转换数据长短 cor 用于gglot2可视化
cor_p2 <- cor_p %>% as.data.frame %>% tibble::rownames_to_column(var = "cell") %>% 
  tidyr::gather(., gene, Pvalue, -cell) #转换数据长短 p 用于ggplot2可视化


cor_dat <- cbind(cor_r2, cor_p2)[,c("cell", "gene", "Correlation", "Pvalue")]
colnames(cor_dat) <- c('cell','gene','correlation','p')
cor_dat$cell <- gsub('.',' ',cor_dat$cell, fix=T)
write.csv(cor_dat, "./data/correlation_cor_exp2.csv")

cor_dat$sig <- cut(cor_dat$p,
                   breaks = c(0, 0.001, 0.01, 0.05, 1),
                   labels = c("***","**","*",""),
                   right=FALSE)
#可视化

p <- ggplot(cor_dat,aes(x = gene, y = cell))+
  geom_tile(aes(fill = correlation), size = 0.5)+
  geom_text(aes(label = sig), color = "black", size = 3) +
  scale_x_discrete(expand=c(0,0)) +
  scale_y_discrete(expand=c(0,0),position="left")+
  theme(text=element_text(family="Roboto"),
        axis.text.x=element_text(color="black"),
        axis.text.y=element_text(color="black"),
        axis.ticks.x = element_blank(),
        axis.ticks.y=element_blank(),
        panel.border = element_rect(fill=NA,color="grey80",
        size=1, linetype="solid"))+
  scale_fill_gradient2(low = "#00008B", high = "#8B0000")

ggsave("./plots/Fig9_heatmap_exp2.png", p, height = 8, width = 5)

复现结果：

图9 G：

文章原图（上） vs CNSknowall R语言复现（下）

图9 H：

解析：图9G在pss样本中，CD8 + T cells与Tregs (原文献r = 0.83，复现r = 0.52，红框框出)和单核monocytes(原文献r = 0.82，复现r = 0.59，黄框框出)呈显著正相关。相比之下，plasma cells 与 γδT cells 呈负相关(原文献r = -0.63，复现r = -0.73，蓝框框出)(图9G)。

    在PD和pSS样本中，γδT cells, activated memory CD4 + T cells, M0 Macrophages 和 naïve B cells 与 CSR2RB、CXCR4和LYN基因表达呈正相关。（在图9F和H中均用红框框出），与原文结果几乎完全一致。

 

 长按二维码关注我们，用最短的时间和最高的效率学习更多数据分析方法!

扫描上方二维码或登录平台官网后添加CNSknowall客服微信咨询！官网地址：

https://cnsknowall.com

CNSknowall：新型科研数据（0代码生信）分析平台，同时含有机制图模块（原创数千素材）+SCI文献例句检索模块+OPenAI官方GPT，500款CNS级别图表统统1秒内1键出图，登录即秒变数据分析大神，体验前所未有的高效便捷数据分析之旅，开启科研天骄之路！

可向下滑动批阅！