4.1 Bulk RNA-seq 下游分析全流程实战

# ================================================================================
# step01 数据标准化与低表达基因过滤
# 输入: GSE52778 原始 counts + 基因长度注释 + 样本分组
# 输出: data/step01-airway.RData (counts, fpkm, tpm, group)
# ================================================================================

# 加载包
rm(list = ls())
library(stringr) # 字符处理包
library(data.table)
library(org.Hs.eg.db)
library(GEOquery)
library(tidyverse)
library(rtracklayer)
library(txdbmaker)

## 1. 读取counts数据
# 来源：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778
counts <- fread("data/GSE52778_raw_counts_GRCh38.p13_NCBI.tsv.gz", data.table = F)
# 列名转行名
counts <- column_to_rownames(counts, "GeneID")
head(counts)
# 过滤在所有样本中表达量均为0的基因
counts <- counts[rowSums(counts)>0,]

## 2. counts转化为标准化值
### 2.1 获取基因长度
# 可以通过上游分析featureCounts的结果来直接获取基因长度，这里我们直接读取
anno <- fread("data/Human.GRCh38.p13.annot.tsv.gz",data.table = F)
colnames(anno)
# 提取需要的列
anno <- anno[,c("GeneID","Symbol","EnsemblGeneID","Length","GeneType")]
head(anno)

### 2.2 counts2fpkm
counts[1:4,1:4] #counts
head(anno) # 基因长度
rownames(anno) <- anno$GeneID
effLen <- anno[row.names(counts), "Length"]
head(effLen)
range(effLen)
# 转化函数
# FPKM/RPKM (Fragments/Reads Per Kilobase Million )  每千个碱基的转录每百万映射读取的Fragments/reads
# RPKM与FPKM分别针对单端与双端测序而言，计算公式是一样的
counts2FPKM <- function(count=count, efflength=efflen){
  PMSC_counts <- sum(count)/1e6   #counts的每百万缩放因子 (“per million” scaling factor) 深度标准化
  FPM <- count/PMSC_counts        #每百万reads/Fragments (Reads/Fragments Per Million) 长度标准化
  FPM/(efflength/1000)
}
fpkm <- apply(counts, 2, counts2FPKM, efflength = effLen)
fpkm[1:5,1:5]

### 2.3 counts2tpm
# TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的Transcripts
counts2TPM <- function(count=count, efflength=efflen){
  RPK <- count/(efflength/1000) #每千碱基reads (reads per kilobase) 长度标准化
  PMSC_rpk <- sum(RPK)/1e6      #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化
  RPK/PMSC_rpk
}
tpm <- apply(counts, 2, counts2TPM, efflength = effLen)
tpm[1:5,1:5]
colSums(tpm) # 每列加和为10^6

## 3. 读取样本分组
gset = getGEO(filename = "data/GSE52778_series_matrix.txt.gz", getGPL = F)
gset
gset[[1]]
pd = pData(gset)
colnames(pd)
group <- pd[,c("geo_accession","title","treatment:ch1")]
group <- group[colnames(tpm), ]
colnames(group)[3] <- "group"
group
table(group$group)
## 这里只用两个分组的样本：Untreated 和 Dexamethasone
group <- group[group$group %in% c("Untreated","Dexamethasone"), ]
group <- group[order(group$group), ]
group

## 4. 低表达过滤
counts <- counts[,rownames(group)]
#查看过滤之前基因数量：
nrow(counts)
# 保留在75%以上样本中表达量>0的基因
counts <- counts[apply(counts,1,function(x) sum(x > 0) > 0.75*ncol(counts)),]
# 查看过滤后的基因数量
nrow(counts)
tpm <- tpm[rownames(counts),colnames(counts)]
fpkm <- fpkm[rownames(counts),colnames(counts)]
# 保存一下
save(counts, fpkm, tpm,group, file = "data/step01-airway.RData")

# ================================================================================
# step02 异常样本检查 (QC): 箱线图 / PCA / 相关性热图
# 输入: data/step01-airway.RData (counts, fpkm, tpm, group)
# 输出: result/qc.boxplot.pdf, result/qc.pca.pdf, result/qc.Cor.pdf
# ================================================================================

rm(list = ls())
library(ggplot2)
library(ggsci) # 配色包
library(cols4all)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(pheatmap)
library(tidyverse)
library(patchwork)
library(ggrepel)

## 加载表达的数据
load(file = "data/step01-airway.RData")
ls()

## 1. 箱线图
## 构造绘图数据,数据变换：宽变长
range(tpm)
exp <- log10(tpm+1)
range(exp)
head(exp)
data <- exp %>%
  as.data.frame() %>%
  pivot_longer(cols = everything(), names_to = "sample",values_to = "expression") %>%
  as.data.frame()
data <- merge(group,data,by.x="geo_accession",by.y="sample")
data <- data[order(data$group), ]
data$geo_accession <- factor(data$geo_accession, levels = unique(data$geo_accession))
head(data)

c4a_gui() #查看色板
mycol <- c4a('area9',4) #配色挑选
mycol

p <- ggplot(data = data, aes(x = geo_accession, y = expression,fill = group, color = group)) +
  geom_boxplot(alpha = 0.6) + #添加箱线图
  scale_fill_manual(values = mycol) + #填充颜色自定义
  scale_color_manual(values = mycol) + #描边颜色自定义
  xlab(label = "") +
  theme_classic() +
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.text.x = element_text(angle = 90))
p
ggsave(filename = "result/qc.boxplot.pdf",height = 6,width = 7)

## 2. PCA分析
data <- log10(tpm+1)
data <- t(data)
pca <- prcomp(data, scale. = F)  # 进行标准化 PCA 计算
# 计算 PCA 解释的方差比例
var_explained <- pca$sdev^2 / sum(pca$sdev^2)
var_explained
sum(var_explained) # 累积和为1
# 提取 PCA 结果
pca_df <- as.data.frame(pca$x[, 1:2])
pca_df$group <- group[rownames(data), "group"]
pca_df$sample <- rownames(pca_df)

# 绘图
p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = group)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_point(data = pca$x, aes(x = PC1,y = PC2), shape = 21, size = 3,
             stroke = 0.5, alpha = 0.8, color = "black" ) + #添加边圈
  geom_text_repel(aes(label = sample), size = 3) +  # 使用 ggrepel 避免重叠
  stat_ellipse(aes(fill = group), alpha = 0.2, geom = "polygon") +
  labs(x = paste0("PC1 (", round(var_explained[1]*100, 2), "%)"),
       y = paste0("PC2 (", round(var_explained[2]*100, 2), "%)")) +
  scale_fill_manual(values = mycol) + #填充颜色自定义
  scale_color_manual(values = mycol) + #描边颜色自定义
  theme_bw() +
  theme(panel.background = element_rect(fill = 'white', colour = 'black'),  # 设置背景为白色，边框为黑色
        panel.border = element_rect( colour = "black", fill = NA, linewidth = 0.8),         # 透明填充 # 边框颜色
        axis.title = element_text(colour = "black", size = 12, margin = margin(t = 12)),
        axis.text = element_text(color = "black"),  # 设置轴刻度文字颜色
        plot.title = element_blank(),  # 不显示标题
        legend.title = element_blank(),  # 不显示图例标题
        legend.key = element_blank(),  # 使图例背景透明
        legend.text = element_text(color = "black", size = 9),  # 设置图例文本颜色和大小
        legend.spacing.x = unit(0.06, 'cm'),  # 设置图例文本之间的水平间距
        legend.key.width = unit(0.01, 'cm'),  # 设置图例的水平大小
        legend.key.height = unit(0.01, 'cm'),  # 设置图例的垂直大小
        legend.background = element_blank(),  # 使图例背景透明
        legend.position = c(0.2, 0.1) # 图例位置
        )
p
ggsave(filename = "result/qc.pca.pdf",height = 4.7,width = 5,plot = p)

## 3. 相关性热图
# 选择差异变化大的基因算样本相关性
# mad: 绝对中位差
# sd:标准差
exp <- tpm
exp <- exp[names(sort(apply(exp, 1, mad),decreasing = T)[1:1000]),]
dim(exp)
# 计算相关性
# 组内的相关性一定大于组间的相关性
# pearson：这种方法对异常值敏感，分析之前往往需要做log10转换
# spearman：这种方法比较的是值与值之间的相对大小关系，log与否对结果没有影响
M <- cor(log10(exp+1),method = "pearson")
M
# 构造注释条
anno <- group[,c("group"),drop=F]
anno

# 保存结果
p <- pheatmap(M,display_numbers = T,
              annotation_col = anno,
              fontsize = 10,
              cellheight = 30,
              cellwidth = 30,
              cluster_rows = T,
              cluster_cols = T)
p <- ggplotify::as.ggplot(p)
p
ggsave(filename = "result/qc.Cor.pdf",height = 6,width = 7.3,plot = p)

# ================================================================================
# 03-2 差异分析: edgeR (TMM 归一化 + GLM 拟合 + 准似然F检验)
# 输入: data/step01-airway.RData
# 输出: result/diff.edgeR.csv, result/diff.edgeR.Rdata
# ================================================================================

rm(list = ls())
library(edgeR)
library(ggplot2)

# 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
load(file = "data/step01-airway.RData")
ls()

# 使用数据为count值
counts[1:4,1:4]

# 分组信息
group
group_list <- group[colnames(counts),"group"]
group_list
# treat vs control
group_list <- factor(group_list,levels = c("Untreated","Dexamethasone")) # 对照在前
group_list
table(group_list)

# 构建线性模型。0代表x线性模型的截距为0
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(counts)
colnames(design) <- levels(factor(group_list))
design

# 构建edgeR的DGEList对象
DEG <- DGEList(counts=counts, group=factor(group_list))
# 归一化基因表达分布
# 对count值标准化的算法  method = "TMM"
DEG <- calcNormFactors(DEG, method = "TMM")
# 计算线性模型的参数
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)
# 拟合线性模型
fit <- glmFit(DEG, design)
# 进行差异分析，c("Untreated","Dexamethasone")) # 对照在前 与 contrast=c(-1,1) 对应
lrt <- glmLRT(fit, contrast=c(-1,1))

# 提取差异分析结果
# 对pvalue进行 fdr 校正 ，方法 adjust.method = "BH"
DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG),adjust.method = "BH"))
head(DEG_edgeR)

# 筛选显著上下调基因，设定阈值
fc_cutoff <- 2 # 1 1.2  1.5  2
#pvalue <- 0.05  # 0.05, 0.01, 0.001
fdr <- 0.05  # 0.05, 0.01, 0.001, 0.2

DEG_edgeR$regulated <- "normal"
DEG_edgeR$regulated[ DEG_edgeR$logFC > log2(fc_cutoff) & DEG_edgeR$FDR < fdr ] <- "up"
DEG_edgeR$regulated[ DEG_edgeR$logFC < -log2(fc_cutoff) & DEG_edgeR$FDR < fdr ] <- "down"

table(DEG_edgeR$regulated)

## 添加一列gene symbol
library(clusterProfiler)
library(org.Hs.eg.db) #不同物种 homo sapiens ，包不一样
keytypes(org.Hs.eg.db)

id2symbol <- bitr(rownames(DEG_edgeR),
                  fromType = "ENTREZID",
                  toType = "SYMBOL",
                  OrgDb = org.Hs.eg.db)
head(id2symbol)

DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
DEG_edgeR_symbol <- merge(id2symbol,DEG_edgeR,
                          by.x="ENTREZID",by.y="GeneID",all.y=T)
head(DEG_edgeR_symbol)

# 选择显著差异表达的结果
library(tidyverse)
DEG_edgeR_symbol_Sig <- filter(DEG_edgeR_symbol,regulated!="normal")

# 保存
write.csv(DEG_edgeR_symbol,"result/diff.edgeR.csv", row.names = F)
save(DEG_edgeR_symbol,file = "result/diff.edgeR.Rdata")

# ================================================================================
# 03-3 差异分析: DESeq2 (RLE 标准化 + 负二项分布 + Wald检验)
# 输入: data/step01-airway.RData
# 输出: result/diff.deseq2.csv, result/diff.deseq2.Rdata
# ================================================================================

rm(list = ls())
library(ggplot2)
library(ggstatsplot)
library(ggsci)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(limma)
library(edgeR)
library(DESeq2)

# 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
load(file = "data/step01-airway.RData")
ls()

# 使用数据为count值
counts[1:4,1:4]

# 分组信息
group
group_list <- group[colnames(counts),"group"]
group_list
# treat vs control
group_list <- factor(group_list,levels = c("Untreated","Dexamethasone")) # 对照在前
group_list
table(group_list)

# 构造deseq样本表型矩阵
colData <- data.frame(row.names=colnames(counts), group_list=group_list)
colData

# 构造deseq2对象
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds)
res <- results(dds, contrast=c("group_list", levels(group_list)[2], levels(group_list)[1]))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG <- as.data.frame(resOrdered)
DEG_deseq2 <- na.omit(DEG)

# 筛选显著上下调基因，设定阈值
fc_cutoff <- 2 # 1 1.2  1.5  2
#pvalue <- 0.05  # 0.05, 0.01, 0.001
fdr <- 0.05  # 0.05, 0.01, 0.001, 0.2

DEG_deseq2$regulated <- "normal"
DEG_deseq2$regulated[ DEG_deseq2$log2FoldChange > log2(fc_cutoff) & DEG_deseq2$padj < fdr ] <- "up"
DEG_deseq2$regulated[ DEG_deseq2$log2FoldChange < -log2(fc_cutoff) & DEG_deseq2$padj < fdr ] <- "down"

table(DEG_deseq2$regulated)

## 添加一列gene symbol
library(clusterProfiler)
library(org.Hs.eg.db) #不同物种 homo sapiens ，包不一样
keytypes(org.Hs.eg.db)

id2symbol <- bitr(rownames(DEG_deseq2),
                  fromType = "ENTREZID",
                  toType = "SYMBOL",
                  OrgDb = org.Hs.eg.db)
head(id2symbol)

DEG_deseq2 <- cbind(GeneID=rownames(DEG_deseq2),DEG_deseq2)
DEG_deseq2_symbol <- merge(id2symbol,DEG_deseq2,
                          by.x="ENTREZID",by.y="GeneID",all.y=T)
head(DEG_deseq2_symbol)

# 选择显著差异表达的结果
library(tidyverse)
DEG_deseq2_symbol_Sig <- filter(DEG_deseq2_symbol,regulated!="normal")

# 保存
write.csv(DEG_deseq2_symbol,"result/diff.deseq2.csv", row.names = F)
save(DEG_deseq2_symbol,file = "result/diff.deseq2.Rdata")

# ================================================================================
# 03-4 差异分析: limma-voom (TMM 标准化 + voom 转换 + 经验贝叶斯)
# 输入: data/step01-airway.RData
# 输出: result/diff.limma.csv, result/diff.limma.Rdata
# ================================================================================

rm(list = ls())
library(ggplot2)
library(ggstatsplot)
library(ggsci)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(limma)
library(edgeR)
library(DESeq2)

# 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
load(file = "data/step01-airway.RData")
ls()
# 使用数据为count值
counts[1:4,1:4]
# 分组信息
group
group_list <- group[colnames(counts),"group"]
group_list
group_list <- factor(group_list,levels = c("Untreated","Dexamethasone")) # 对照在前
group_list
table(group_list)

## 构造对比矩阵
design <- model.matrix(~0+group_list)
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(counts)
design

g1 <- levels(group_list)[1]
g2 <- levels(group_list)[2]
con <- paste0(g2,'-',g1)
cat(con)
cont.matrix <- makeContrasts(contrasts=c(con),levels = design)

## counts标准化以及差异分析
dge <- DGEList(counts=counts)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
v <- voom(dge,design,plot=TRUE, normalize="none")
# 因为 TMM 一般已经足够处理文库大小和组成偏倚额外使用
# quantile normalization 可能过度强制各样本具有相同分布
# 尤其当两组之间确实存在广泛生物学表达差异时，可能削弱真实信号
fit <- lmFit(v, design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)

## 提取差异结果
DEG_limma_voom <- topTable(fit2, coef=con, n=Inf,adjust.method = "BH")
DEG_limma_voom <- na.omit(DEG_limma_voom)
head(DEG_limma_voom)

## 筛选显著上下调基因，设定阈值
fc_cutoff <- 2 # 1 1.2  1.5  2
#pvalue <- 0.05  # 0.05, 0.01, 0.001
fdr <- 0.05  # 0.05, 0.01, 0.001, 0.2

DEG_limma_voom$regulated <- "normal"
DEG_limma_voom$regulated[ DEG_limma_voom$logFC > log2(fc_cutoff) & DEG_limma_voom$adj.P.Val < fdr ] <- "up"
DEG_limma_voom$regulated[ DEG_limma_voom$logFC < -log2(fc_cutoff) & DEG_limma_voom$adj.P.Val < fdr ] <- "down"

table(DEG_limma_voom$regulated)

## 添加一列gene symbol
library(clusterProfiler)
library(org.Hs.eg.db) #不同物种 homo sapiens ，包不一样
keytypes(org.Hs.eg.db)

id2symbol <- bitr(rownames(DEG_limma_voom),
                  fromType = "ENTREZID",
                  toType = "SYMBOL",
                  OrgDb = org.Hs.eg.db)
head(id2symbol)

DEG_limma_voom <- cbind(GeneID=rownames(DEG_limma_voom),DEG_limma_voom)
DEG_limma_voom_symbol <- merge(id2symbol,DEG_limma_voom,
                           by.x="ENTREZID",by.y="GeneID",all.y=T)
head(DEG_limma_voom_symbol)

# 选择显著差异表达的结果
library(tidyverse)
DEG_limma_voom_symbol_Sig <- filter(DEG_limma_voom_symbol,regulated!="normal")

# 保存
write.csv(DEG_limma_voom_symbol,"result/diff.limma.csv", row.names = F)
save(DEG_limma_voom_symbol,file = "result/diff.limma.Rdata")

# ================================================================================
# 03-5 差异结果可视化: 热图 / 火山图 / 箱线图
# 输入: data/step01-airway.RData + result/diff.edgeR.Rdata
# 输出: result/DEG_heatmap.pdf, result/DEG_heatmap-1.pdf,
#       result/DEG_Volcano.pdf, result/DEG_Volcano2.pdf, result/CD163.boxplot.pdf
# ================================================================================

# ---------------------------------------------------------------
# 5.1 热图
# 为什么不可以用原始counts值，为什么可以用FPKM/TPM值？
# 热图用于基因之间比较 + 样本之间比较，需双重校正：
#   基因间比较 → 对基因长度矫正
#   样本间比较 → 对文库大小矫正
# FPKM/TPM 对二者都做了矫正，适用于绘制热图
# ---------------------------------------------------------------
rm(list = ls())
library(pheatmap) # complexheatmap
library(tidyverse)
library(patchwork)

# 加载原始表达矩阵
load(file = "data/step01-airway.RData")
# 读取差异分析结果
load(file = "result/diff.edgeR.Rdata")
ls()

# 提取所有差异表达的基因名
diff <- DEG_edgeR_symbol
head(diff)
diff <- diff[diff$regulated!="normal","ENTREZID"]
head(diff)
# 提取对应的表达矩阵
data <- tpm[diff,]
head(data)
data[1:4,1:4]

# 分组注释信息
group
anno <- group[,"group",drop=F]
anno

# 定义注释条颜色
library(cols4all)
c4a_gui() #查看色板
mycol <- c4a('light',2) #配色挑选
names(mycol) <- c('Dexamethasone','Untreated')
mycol
#"#77AADD"     "#EE8866"
ann_colors = list( group=mycol )
ann_colors

# 绘制热图
p <- pheatmap(data,scale = "row",show_colnames =F,
         show_rownames = F, cluster_cols = T,
         annotation_col=anno,
         annotation_colors = ann_colors,
         col = colorRampPalette(c("#ceb0ff","white","#ff4628"))(100),
         treeheight_row = 0,treeheight_col = 0,
         main = "edgeR's DEG")
p <- ggplotify::as.ggplot(p)
p
ggsave(filename = "result/DEG_heatmap.pdf",height = 6,width = 5)

## 热图款式2
data[1:4,1:4]
anno
ann_color <- list(group =  c(Dexamethasone="#ed1a22",Untreated="#00a651"))
library(RColorBrewer)
color <- colorRampPalette(c('#1f67ad','white','#bb2f34'))(100)
# 看gap区域
table(group$group)
p <- pheatmap(data, scale = "row",show_rownames = F,show_colnames = F,
         color = color,annotation_col = anno, annotation_colors = ann_color,
         cluster_cols = F,
         border_color = "black",gaps_col=c(4))
p <- ggplotify::as.ggplot(p)
p
ggsave(filename = "result/DEG_heatmap-1.pdf",height = 6,width = 5)

# ---------------------------------------------------------------
# 5.2 火山图
# ---------------------------------------------------------------
rm(list = ls())
library(ggplot2)
library(tidyverse)
library(ggrepel)
library(ggpubr)
library(ggthemes)

# 读差异分析结果
load(file = "result/diff.edgeR.Rdata")

# 根据需要修改DEG的值
data <- DEG_edgeR_symbol
head(data)

# 绘制火山图
p <- ggplot(data=data, aes(x=logFC, y=-log10(FDR),color=regulated)) +
  geom_point(alpha=0.6, size=1.5) +
  theme_bw(base_size=20) +
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank()) +
  xlab("log2FC") +
  ylab("-log10(FDR)") +
  scale_colour_manual(values = c(down="#34bfb5",normal="#828586",up="#ff6633")) +
  geom_vline(xintercept=c(-(log2(2)),log2(2)),lty=4,col="grey",lwd=0.6) +
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)
p

# 添加top基因
#label <- data[order(abs(data$logFC),decreasing = T)[1:10],] # 通过FC选取TOP10
label <- data[order(abs(data$PValue),decreasing = F)[1:10],] # 通过pvalue选取TOP10
label <- na.omit(label)
label

p1 <- p +
  geom_point(size = 3, shape = 1, data = label,color="red") +
  geom_text_repel( aes(label = SYMBOL), data = label, color="black" )
p1
ggsave(filename = "result/DEG_Volcano.pdf", width = 6.4, height = 5, plot = p1)

# 还可以是指定的感兴趣的关键基因
cg <- c("TP53","EGFR")
label <- data[data$SYMBOL %in% cg, ]
label
p2 <- p +
  geom_point(size = 3, shape = 1, data = label,color="red") +
  geom_text_repel( aes(label = SYMBOL), data = label, color="black" )
p2

## 火山图款式2
edgeR.res <- DEG_edgeR_symbol %>% mutate(regulated=factor(regulated,levels = c("up","down","normal")),negLog10FDR = -log10(pmax(FDR, .Machine$double.xmin)))
head(edgeR.res)

p <- ggscatter(edgeR.res, x = "logFC", y = "negLog10FDR", color = "regulated", size = 2,
               main = paste0("Differential expressed genes\nbetween Dexamethasone and Untreated"), # ***
               xlab = "",
               ylab = "-log10(adjusted P.value)",
               palette = c("#D01910","#00599F","#CCCCCC"),
               ylim = c(-1, 70), xlim=c(-8,8)) +
  geom_hline(yintercept = -log10(0.05), linetype="dashed", color = "#222222") +
  geom_vline(xintercept = 0.58 , linetype="dashed", color = "#222222") +
  geom_vline(xintercept = -0.58, linetype="dashed", color = "#222222") +
  xlab(label = expression(log[2](fold~change))) +
  labs(color = "Significance") +  # 修改颜色图例的标题
  annotate("text", x = -5, y = 66, label = "496 genes", size = 5, color = "#00599F", hjust = 0.5) + # 设置文本样式
  annotate("text", x = 5, y = 66, label = "539 genes", size = 5, color = "#D01910", hjust = 0.5) + # 设置文本样式
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),  # 去掉主要格子线
    panel.grid.minor = element_blank(),  # 去掉次要格子线
    plot.background = element_blank(),
    legend.position = "bottom"
  )
p
ggsave("result/DEG_Volcano2.pdf", p, width = 5.8, height = 6)

# ---------------------------------------------------------------
# 5.3 箱线图
# ---------------------------------------------------------------
rm(list = ls())
library(ggplot2)
library(ggsci) # 配色包
library(cols4all)

## 加载表达的数据
load(file = "data/step01-airway.RData")
ls()

## 构造绘图数据,数据变换：宽变长
range(tpm)
exp <- log10(tpm+1)
range(exp)
head(exp)
data <- exp %>%
  as.data.frame() %>%
  filter(rownames(.) == "9332") %>% # 以CD163为例
  pivot_longer(cols = everything(), names_to = "sample",values_to = "expression") %>%
  as.data.frame()
data <- merge(group,data,by.x="geo_accession",by.y="sample")
data <- data[order(data$group), ]
data$group <- factor(data$group,levels = unique(data$group))

# 进行t检验
t_test_result <- t.test(expression ~ group, data = data)
p_value <- t_test_result$p.value
p_value # 0.08008757
p_sci <- sprintf("%.2e", p_value)  # "8.01e-02"
p_sci
p_expression <- bquote(italic(p) == .(p_sci))
p_expression

# 手动添加显著性标记
y_max <- max(data$expression) * 1.1
y_line <- y_max * 0.95
y_line
y_text <- y_max * 1.02
y_text

p <- ggplot(data, aes(x = group, y = expression, fill = group)) +
  # 箱线图
  geom_boxplot(alpha = 0.6, outlier.shape = NA, width = 0.4,lwd=1.2) +
  # 抖动散点
  geom_jitter(aes(color = group),  width = 0.15, alpha = 0.8, size = 3) +
  # 手动添加短横线
  geom_segment(aes(x = 1, xend = 2, y = y_line, yend = y_line), color = "grey20", linewidth = 0.8) +
  # 手动添加p值文本（斜体p）
  annotate("text", x = 1.5, y = y_text, label =  deparse(p_expression), parse = TRUE, size = 7, color = "grey20",fontface = "bold") +
  # 美化设置
  scale_fill_manual(values = c("Dexamethasone" = "#fad5a0", "Untreated" = "#b5bcd6")) +
  scale_color_manual(values = c("Dexamethasone" = "#f4b75e",  "Untreated" = "#8492bb")) +
  scale_y_continuous(limits = c(0, NA)) +  # y轴从0开始
  labs(title = "GSE52778", x = "Group", y = "CD163 Expression") +
  theme_classic() +
  theme(
    legend.position = "none",
    # 标题灰黑色加粗加大
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, color = "grey20"),
    # 坐标轴标题灰黑色加粗加大
    axis.title.x = element_text(size = 16, face = "bold", margin = margin(t = 10), color = "grey20"),
    axis.title.y = element_text(size = 16, face = "bold", margin = margin(r = 10), color = "grey20"),
    # 坐标轴刻度标签灰黑色加粗加大
    axis.text.x = element_text(size = 14, face = "bold", color = "grey20"),
    axis.text.y = element_text(size = 12, face = "bold", color = "grey20"),
    # 坐标轴线灰黑色加粗
    axis.line = element_line(linewidth = 1.2, color = "grey20"),
    # 坐标轴刻度线灰黑色加粗
    axis.ticks = element_line(linewidth = 1.2, color = "grey20"),
    axis.ticks.length = unit(0.2, "cm")
  )
p
ggsave(filename = "result/CD163.boxplot.pdf",height = 6,width = 7)

# ================================================================================
# 04-2 功能富集分析: ORA (过代表分析) — GO / KEGG / Hallmark
# 输入: result/diff.edgeR.Rdata + data/MsigDB/*.gmt
# 输出: result/Enrich.ORA_*.csv + Enrich.ORA_KEGG_GO.pdf + Enrich.ORA_HALLMARK.pdf
# ================================================================================
# 过代表分析 (Over-Representation Analysis, ORA)
# 在一组显著表达的基因中，推断受干扰的生物途径/过程
# 统计学原理：累计超几何分布

rm(list = ls())
library(clusterProfiler) #功能富集分析包，统计学原理累计超几何分布
packageVersion("clusterProfiler")
# devtools::install_local("packages/GOSemSim-devel.zip")
# devtools::install_local("packages/clusterProfiler-devel.zip")
library(org.Hs.eg.db) # 数据包，做基因ID转换，是物种特异性包Hs表示人
library(GSEABase) # 读取gmt格式数据
library(ggplot2)
library(tidyverse)

## 1.读取差异分析结果----
load(file = "result/diff.edgeR.Rdata")
ls()
# 提取所有差异表达的基因名
DEG <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal",2]
DEG <- na.omit(DEG)
head(DEG)

## 2.基于GO数据库的功能富集分析----
# 输出所有结果,包括显著和不显著，后续可根据pvalue挑选结果
# org.Hs.eg.db 为物种特异性包,保存的是通路与基因的对应关系，是一个数据包
# ont="ALL"，一次性分析上面三个本体论
ego_CC <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL',
                   ont="CC", pvalueCutoff= 1,qvalueCutoff= 1,
                   minGSSize = 1, maxGSSize = 10000)
ego_MF <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL',
                   ont="MF", pvalueCutoff= 1,qvalueCutoff= 1,
                   minGSSize = 1, maxGSSize = 10000)
ego_BP <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL',
                   ont="BP", pvalueCutoff= 1,qvalueCutoff= 1,
                   minGSSize = 1, maxGSSize = 10000)
ego_BP

colnames(ego_BP@result)
## 绘制条形图
# showCategory = 10：绘制top10通路
# color = "p.adjust", 按照p.adjust从小到大排序的top10
# label_format：通路名字达到多少个字符之后会换行，默认是30
# dotplot(ego_BP,showCategory = 10,label_format = 100)
p_BP <- barplot(ego_BP,showCategory = 10,label_format = 100) + ggtitle("Biological process")
p_CC <- barplot(ego_CC,showCategory = 10,label_format = 100) + ggtitle("Cellular component")
p_MF <- barplot(ego_MF,showCategory = 10,label_format = 100) + ggtitle("Molecular function")
p_BP
p_CC
p_MF

# 设置color参数按照pvalue展示top15条通路
colnames(ego_MF@result)
barplot(ego_MF,showCategory = 15,label_format = 100,color = "pvalue")

## 3.基于KEGG数据库的功能富集----
#genelist <- bitr(gene=DEG, fromType="SYMBOL", toType="ENTREZID", OrgDb='org.Hs.eg.db')
genelist <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal",1]
head(genelist)
options(timeout = 9999)
# https://www.kegg.jp/kegg/tables/br08606.html organism：不同物种参数
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa',
                    pvalueCutoff = 1, qvalueCutoff = 1)
ekegg <- DOSE::setReadable(ekegg, OrgDb='org.Hs.eg.db', keyType='ENTREZID')#ENTREZID替换symbol
#saveRDS(ekegg,file = "result/ekegg.Rds")
ekegg <- readRDS("result/ekegg.Rds")

## 解决联网报错1
if(F) {
  # enrichKEGG 函数需要连网去KEGG数据库下载基因和通路的对应关系
  # 三种方法解决KEGG功能富集总是出现的断网报错：cannot read from connection
  # https://mp.weixin.qq.com/s/W4qIj7Q2kvDndJpxUmsrrQ
  # 需要去msigdb数据库下载
  geneset <- read.gmt("data/MsigDB/c2.cp.kegg_legacy.v2026.1.Hs.symbols.gmt")
  table(geneset$term)
  ## TERM2GENE： 可以自定义基因集,基因与通路的对应关系,这里DEG与geneset里面的id类型相同
  ekegg <- enricher(gene=DEG, pvalueCutoff = 1, qvalueCutoff = 1, TERM2GENE=geneset,
                      minGSSize = 1,maxGSSize = 10000)
}

## 绘制条形图
p_kegg <- barplot(ekegg, showCategory=10,label_format = 100) + ggtitle("KEGG Pathway")
p_kegg

## 合并保存
library(patchwork)
## 看看有没有不显著的通路，有的话，后面可以做个性化绘图
p <- (p_BP/p_CC) |(p_MF/p_kegg)
p
ggsave('result/Enrich.ORA_KEGG_GO.pdf', plot = p, width = 16, height = 9)

### 保存全部的富集结果
ego_BP <- data.frame(ego_BP@result)
ego_CC <- data.frame(ego_CC@result)
ego_MF <- data.frame(ego_MF@result)
ekegg  <- data.frame(ekegg@result)
write.csv(ego_BP,'result/Enrich.ORA_GO_BP.csv')
write.csv(ego_CC,'result/Enrich.ORA_GO_CC.csv')
write.csv(ego_MF,'result/Enrich.ORA_GO_MF.csv')
write.csv(ekegg ,'result/Enrich.ORA_KEGG.csv')

# 挑选有统计学意义的通路,可以方便的更改挑选阈值
ego_BP_sig <- ego_BP[ego_BP$p.adjust < 0.05, ]
ego_CC_sig <- ego_CC[ego_CC$p.adjust < 0.05, ]
ego_MF_sig <- ego_MF[ego_MF$p.adjust < 0.05, ]
ekegg_sig <- ekegg[ekegg$p.adjust<0.05, ]
write.csv(ego_BP_sig,'result/Enrich.ORA_GO_BP.sig.csv')
write.csv(ego_CC_sig,'result/Enrich.ORA_GO_CC.sig.csv')
write.csv(ego_MF_sig,'result/Enrich.ORA_GO_MF.sig.csv')
write.csv(ekegg_sig ,'result/Enrich.ORA_KEGG.sig.csv')

## 4.其他数据库功能富集----
# hs 为人，其他物种需要去msigdb数据库下载，随便选择一个通路
geneset <- read.gmt("data/MsigDB/h.all.v2026.1.Hs.symbols.gmt")
table(geneset$term)
## TERM2GENE： 可以自定义基因集,基因与通路的对应关系
my_path <- enricher(gene=DEG, pvalueCutoff = 1, qvalueCutoff = 1, TERM2GENE=geneset,
                    minGSSize = 1,maxGSSize = 10000)
p1 <- barplot(my_path, showCategory=15,color = "pvalue",label_format = 100) +
  ggtitle("hallmark gene sets")
p1
ggsave("result/Enrich.ORA_HALLMARK.pdf", plot = p1, width = 7.5, height = 4.8)
# 保存
my_path <- data.frame(my_path)
write.csv(my_path,"result/Enrich.ORA_HALLMARK.csv")

# ================================================================================
# 04-3 功能富集分析: GSEA (基因集富集分析) + Science 风格美化点图
# 输入: result/diff.edgeR.Rdata + data/MsigDB/h.all.v2026.1.Hs.symbols.gmt
# 输出: result/Enrich_GSEA.pdf + result/Enrich_GSEA_Dotplot.pdf
#       result/Enrich_GSEA_fc.csv + result/Enrich_GSEA_fc_sig.csv
# ================================================================================
# 原理：利用表达数据计算基因在AB组的差异，按差异排序（log2FC），
# 检验基因集是随机分布还是更集中在某一端（KS 检验，加权排序）

# 清空当前环境变量
rm(list = ls())
library(GSEABase)
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(stats)
library(ggrepel)

## 1.加载数据----
load("result/diff.edgeR.Rdata")
ls()
DEG <- DEG_edgeR_symbol

## 2.构造GSEA分析数据----
# 去掉没有配对上symbol的行
DEG <- DEG[!is.na(DEG$SYMBOL),]
# 去掉重复行
DEG <- DEG[!duplicated(DEG$SYMBOL),]
# 获取排序的基因列表，排序指标log2FC
geneList <- DEG$logFC
names(geneList) <- DEG$SYMBOL
head(geneList)
geneList <- sort(geneList,decreasing = T) # 降序排列
head(geneList)
tail(geneList)

## 3.选择gmt文件----
geneset <- read.gmt("data/MsigDB/h.all.v2026.1.Hs.symbols.gmt")
stat <- as.data.frame(table(geneset$term))
stat

## 4.运行,输出全部结果----
#？GSEA 查看默认值
egmt <- GSEA(geneList, TERM2GENE=geneset, pvalueCutoff = 1, minGSSize = 1,maxGSSize = 10000)
colnames(egmt@result)

#出点图
dotplot(egmt, label_format = 100)
# 按p值出点图,设置color参数
dotplot(egmt, color="pvalue")
# 单个通路图,按照通路名
egmt@result$Description[1]
gseaplot2(egmt, egmt@result$Description[1],  title = egmt@result$Description[1])

# 按照行数中的第五行
egmt@result$Description[5]
egmt@result[5,]
gseaplot2(egmt, 5, color="black", pvalue_table = T)

# 按egmt@result中的第1:6行出图，不显示p值
p <- gseaplot2(egmt, 1:6, color =  RColorBrewer::brewer.pal(6, "Set1") )
p
ggsave('result/Enrich_GSEA.pdf', plot = p, width = 8, height = 6)

## 6.结果美化----
# Science杂志高颜值GSEA打分排序图
# https://mp.weixin.qq.com/s/udzFxcwSpqJACkPM26oHlw
# 取出绘图需要的列，并进行相关设置：
data <- egmt[,c("ID", "NES","setSize","pvalue")]
data$setSize_1 <- data$setSize/10
head(data)
# 给Y轴的通路名设置为因子，排序
data <- data[order(data$NES, decreasing = T),]
data$ID <- factor(data$ID, levels = data$ID)
data$xlab <- 1:nrow(data)
head(data)
summary(data$NES)
summary(data$setSize_1)

# 图中标出的通路名字
label <- c( "HALLMARK_ADIPOGENESIS", "HALLMARK_OXIDATIVE_PHOSPHORYLATION",
            "HALLMARK_IL2_STAT5_SIGNALING", "HALLMARK_KRAS_SIGNALING_UP",
            "HALLMARK_NOTCH_SIGNALING" )
data_label <- data[data$ID %in% label,]
data_label
# 图中通路的颜色
data_label$col <- c("black", "grey80","grey80","#ffb882","#ff5eff")

p <- ggplot(data = data, aes(x = xlab, y = NES)) +
  geom_point(aes(size = setSize_1, alpha = -log10(pvalue)), shape = 21, stroke = 0.7, fill = "#0000ff", colour = "black") + # stroke：设置点的边框宽度
  scale_size_continuous(range = c(0.2, 8)) +
  xlab(label = "Hallmark gene sets") +
  ylab(label = "Normalized enrichment score (NES)") +
  theme_classic(base_size = 15) +
  scale_x_continuous(breaks = seq(0, 50, by = 10), labels = seq(0, 50, by = 10) ) + # 设置x轴的刻度线和刻度标签
  scale_y_continuous( breaks = seq(-4, 2.3, by = 1), labels = seq(-4, 2.3, by = 1) ) +
  guides(size = guide_legend(title = "setSize"),
         alpha = guide_legend(title = "Significance\n(-log10 p-val.)") ) +  # 修改图例标题
  theme(
    axis.line = element_line(color = "black", size = 0.6), # 加粗x轴和y轴的线条
    axis.text = element_text(face = "bold"), # 加粗x轴和y轴的标签
    axis.title = element_text( size = 13)    # 加粗x轴和y轴的标题
  )
p

# 添加label：vjust（垂直调整）或hjust（水平调整）
p2 <- p +
  geom_text_repel(data= data_label,
                  aes(x = xlab, y = NES,label = ID),
                  size = 3, color = data_label$col,
                  force=20,                # 将重叠的文本标签相互推开的强度。force 参数的值越大，标签之间的排斥力度也越大
                  point.padding = 0.5,     # 设置文本标签与对应点之间的最小距离
                  min.segment.length = 0,  # 长度大于0就可以添加引线
                  hjust = 1.2,             # 文本标签的右侧与指定位置对齐
                  segment.color="grey20",
                  segment.size=0.3,        # 设置引导线的粗细
                  segment.alpha=0.8,       # 文本标签中连接线段的透明度
                  nudge_y=-0.1             # 在y轴方向上微调标签位置
  )
p2
ggsave(filename = "result/Enrich_GSEA_Dotplot.pdf", plot = p2, width = 6.2, height = 5)

## 7.保存结果
gsea <- as.data.frame(egmt@result)
# 富集到的通路
gsea_sig <- gsea[gsea$pvalue<0.05,]
write.csv(gsea,"result/Enrich_GSEA_fc.csv",row.names = F)
write.csv(gsea_sig,"result/Enrich_GSEA_fc_sig.csv",row.names = F)