前面进行了上游分析，接下来将在R中进行下游分析。Bulk RNA-seq最主要的下游分析就是进行差异基因表达分析，所以这一部分只讲解标准的下游分析流程。

首先读入数据

#设置工作目录
setwd("/home/GSE184771/4.stringtie)
getwd()
#加载包
#DESeq2 的统计模型假设输入的数据必须是原始计数数据（整数值），转换后的数据（如 VST 或 rlog 转换数据、标准化数据都不可以）
library(DESeq2) 
library(pheatmap)
library(ggplot2)
#读取数据
countdata <- read.csv("./gene_count_matrix.csv",row.names="gene_id") #设置行名
View(countdata)
countdata <- as.matrix(countdata) #转化为矩阵
countdata <- countdata[rowMeans(countdata)>1,] #数据中会有大量不表达/表达量为0的基因，可以筛除这些行
condition <- factor(c(rep("control",3),rep("exepriment",3))) #提前构建分组信息，也可以直接输入分组信息表格，然后提取condition列，转化为因子。
coldata<-data.frame(row.names=colnames(countdata), condition) #整理成一个样本信息表格，counts表格的列名与样本信息表格的行名必须相同
coldata$batch <- factor(c(rep("first",6),rep("second",3))) #如果有不同批次的信息也可以加入，方便后续进行主成分分析、绘制热图等
View(coldata)

进行差异表达分析

#dds: DESeq2的对象，由 DESeqDataSetFromMatrix函数创建。可通过模型公式（design formula）将批次信息纳入差异表达分析模型，以矫正批次效应，但效果不一定好
dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata,design = ~ batch + condition)
head(dds)
dds$condition
dds$condition <- relevel(dds$condition, ref ="control")#直接指定对照组作为参考水平(对照组在前)，此种方法更直观一点；
#使用DESeq()函数进行标准化、差异分析
dds <- DESeq(dds,parallel = T)
sizeFactors(dds)##包含了每个样本的大小因子，是一个数值，用于将该样本的计数数据调整到可以与其他样本比较的相同尺度上

需要注意的是如果，coldata$batch中的对照组和实验组是不同批次的，那么运行
dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata,design = ~ batch + condition)，会报错，错误提示设计矩阵不是满秩的，因为batch和condition是成对出现的，只有control first和experiment second，没有control second 和experiment first，所以batch和condition完全共线性

报错

去除批次效应，并进行主成分分析，观察实验组和对照组有没有分开，即组间差异大于组内差异，选择能够分开的组进行差异表达分析

#原始的count和归一化之后的count, 其PCA图是杂乱无序的，而VST或rlog转换之后，生物学重复之间更佳的接近，不同分组也区分的较为明显。
#DEG只是要比较不同样本间同一个基因是否有差异，只把counts在样本内做了标准化，从而使不同样本的同一个基因具有可比性。使用 counts标准化，然后再进行scale，PCA效果接近rlog，但不是完全一样。
#PCA，聚类等分析不仅要比较不同样本间同一个基因的差异，还要计算同一个样本内不同基因的贡献，如果直接用归一化数据做分析，会导致样本内表达量大的基因对结果影响过大。但是如果简单用log+1的方式进行转换，又会导致表达量小的基因影响过大，因此DESeq2提出了rlog和VST的方法
#vst函数需要一个DESeqDataSet对象作为输入，因为vst函数在执行转换时需要利用样本的分组信息和其他元数据来正确地处理和归一化计数数据。而DESeqDataSet对象包含了实验设计的信息（如分组信息、批次效应等），还会对原始计数数据进行一些初步的处理，如数据过滤、标准化等。
vst <- vst(dds) #vst转换后
vst_count <-assay(vst)

####################移除批次效应用于下游分析，如绘制热图、PCA等，使用vst转换后的数据并移除批次效应，但校正后的数据仅用于可视化和探索性分析，而不是差异表达分析。####################
#model.matrix 要基于coldata$condition这个变量来生成设计矩阵。
mod <- model.matrix(~factor(coldata$condition))
dds_batch <- removeBatchEffect(vst_count,batch = coldata$batch,design=mod)
#对计数矩阵进行转置，使得样本作为行，基因作为列。scale. = TRUE: 在进行PCA之前对数据进行标准化处理，使得每个基因的方差为1。
pca_result_batch <- prcomp(t(dds_batch), scale. = TRUE)#在进行PCA之前，检查并移除那些具有恒定值（所有值相同）或全零的列。这些列在统计分析中没有信息量，无法用于计算方差或协方差。
View(pca_result_batch)

#画PCA图
library(ggrepel) 
pca_result_batch <-as.data.frame(pca_result_batch$x)
pca_result_batch$group <- c("C", "C", "C", "E", "E", "E","E","E","E")
#添加样本点，每个点的大小为3。手动设置形状。设置图例位置在顶部，并移除图例标题。使用 ggrepel 包添加样本标签，防止文本重叠

pca_batch=ggplot(data = pca_result_batch, aes(x = PC1, y = PC2,color=group,shape = group)) +
#添加置信椭圆，但是点太少画不成，可能需要自己写椭圆函数
  #stat_ellipse(level = 0.95, show.legend = F) +
  geom_point(size=3) +scale_shape_manual(values = c(1:14))+ 
  theme_classic()+
  theme(legend.position="top", legend.title=element_blank())+
  geom_text_repel(aes(label=rownames(pca_result_batch)), show.legend=F, max.overlaps =30, size=3)+#防止文本重合，如果标签太多且重叠过多（超过 max.overlaps 设置的值），则会隐藏部分标签
  labs(title = "PCA of RNA-seq Data", x = "PC1", y = "PC2")
pca_batch
ggsave("PCA_plot_batch.pdf",plot=pca_batch,width = 8,height = 7)

继续进行差异表达分析

##使用results()函数提取基因差异表达的分析结果；
res <- results(dds)
summary(res)
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
# 依次按照pvalue值log2FoldChange值进行排序。行根据 pvalue 的升序排列（即 p 值小的排在前面），并且在 pvalue 相同的情况下按照 log2FoldChange 的降序排列
res1 <- res1[order(res1$pvalue, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
#由于gene id和gene name以|连接，直接它们按照|分开后，gene name有重复。编辑输出的表格，用Excel的数据分列功能将它们分开，并将重复的gene name写为gene name-gene id，保存为csv格式再上传读入
write.table(res1,file = "res1.xlsx",quote = F,row.names = T)
res1<-read.csv("./res1.csv")
res1 <-as.data.frame(res1)
View(res1)
#去除gene id 中"."及后面的数字
res1$gene_id <-gsub("\\..*","",res1$gene_id)
rownames(res1)<-res1[,2] #将gene_name列作为行名
res1<-res1[,-2] #去除多余gene_name列
View(res1)
# 筛选pvalue和padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
res1_DEG <-  subset(res1, abs(log2FoldChange) > 1 & pvalue < 0.05 & padj < 0.05)
res1_up<- res1_DEG[which(res1_DEG$log2FoldChange >= 1 ),]      # 表达量显著上升的基因
res1_down<- res1_DEG[which(res1_DEG$log2FoldChange <= -1),]    # 表达量显著下降的基因
write.table(res1_up,file = "DEG_up.xlsx",quote = F,row.names = T)
write.table(res1_down,file = "DEG_down.xlsx",quote = F,row.names = T)

res1.csv

res1处理后

绘制火山图和前20个变化明显基因的热图

#画火山图，对象是所有基因res1
genes <- res1
# 根据上调、下调、不变为基因添加颜色信息。ifelse条件语句padj<0.05且log2FC的绝对值大于等于1，在前面条件都满足的条件下再嵌套log2FC>1为red，否则为蓝色，不满足前面条件的为灰色
genes$color <- ifelse(genes$padj<0.05 & abs(genes$log2FoldChange)>= 1,ifelse(genes$log2FoldChange >= 1,'red','blue'),'gray')
color <- c(red = "red",gray = "gray",blue = "blue",green="pink") ###需要额外定义颜色向量，否则可能会导致颜色向量模糊，定义不清
p <- ggplot(
  genes, aes(log2FoldChange, -log10(padj), col = color)) +  
  geom_point(size=1) + 
  theme_bw() + #设置图形主题为白底
  scale_color_manual(values = color) + #手动指定颜色映射，color变量前面已定义
  labs(x="log2 (fold change)",y="-log10 (qvalue)") + #y轴和x轴标签
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) + #添加水平虚线，设置虚线水平位置为q=0.05的位置，线型为虚线，虚线线宽0.6
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) + #添加垂直虚线
  # 图例
  theme(legend.position = "none", #隐藏图例，颜色已通过color映射
        panel.grid=element_blank(), #隐藏网格线
        axis.title = element_text(size = 16), #设置轴标题文本大小
        axis.text = element_text(size = 14) #设置轴刻度标签文本大小
  )
p
ggsave("volcano.pdf",plot=p,width = 12.5,height = 10.5) #保存火山图到当前工作目录

#选择几个感兴趣的基因，将其在图中标注出来
marker_genes <- c("Bcl6","Ccdc8","Asgr2")
genes$gene_name <- rownames(genes)
marker_genes <- genes[marker_genes,]
marker_genes$subset_col <- c("green")
#ggplot是图层绘图，一层一层叠加绘图元素，makergene的火山图是在p的基础上标记出几个基因的位置，颜色用粉色区分出来
p1 <- p +  geom_point(data=marker_genes,aes(x=log2FoldChange, y=-log10(padj),color=subset_col))+
  scale_color_manual(values = color)+
  #添加文本标签，可防止重叠
  geom_text_repel(data = marker_genes, 
    # #文本标签的内容是基因名称，即前面定义的行名
    aes(x=log2FoldChange, y=-log10(padj),label =gene_name),
    #文本标签的大小为4
    size = 4, 
    color="black",
    #控制连线长短
    box.padding = unit(1.5, "lines"), 
    #遮蔽点和标签，越大连线越短
    point.padding = unit(1, "lines"),  
    min.segment.length = 0,
    segment.color = "black",
    max.overlaps = Inf
  )
p1
ggsave("makergenes_volcano.pdf",plot=p1,width = 12.5,height = 10.5)

###绘制前20个变化明显基因的热图

#将上调和下调的前20个基因从差异分析表格中取出子集
sorted_up <- res1_up[order(-res1_up$log2FoldChange,log10(res1_up$pvalue)),]
sorted_down <- res1_down[order(res1_down$log2FoldChange,log10(res1_down$pvalue)),]
top10_up <- sorted_up[1:10,]
top10_down <- sorted_down[1:10,]
#合并行
top_total <- rbind(top10_up,top10_down)
View(top_total)
top_total$gene_name <- rownames(top_total)

#使用count矩阵绘图，将这20个基因在count矩阵中取出子集
countdata <- as.data.frame(countdata)
countdata$gene_name <- rownames(countdata)
#将gene_id去除
countdata$gene_name <- sub(".+\\|","",countdata$gene_name)
common_genes <- intersect(countdata$gene_name, top_total$gene_name)#不为空
df <- countdata[countdata$gene_name %in% common_genes, ]
View(df)
rownames(df) <-df$gene_name
df <- df[,-10]
p_heatmap <-pheatmap(df,
              show_rownames = T,
              show_colnames = T,
              cluster_cols = F,
              cluster_rows=T, #行聚类
              height=100,  
              #对行进行标准化处理
              scale = "row",
              frontsize = 5,
              angle_col=45, 
              color =colorRampPalette(c("#8854d0", "#ffffff","#fa8231"))(100), #生成包含100种颜色的梯度，从紫色过渡到白色再到橙色
              clustering_method = "single",
              border_color = "lightgray",
              border=T
) 
p_heatmap
ggsave("top10_heatmap.pdf",plot=p_heatmap,width = 14,height = 12)

volcano,pdf

makergenes_valcano.pdf

其实可以看出c组和E组的组内重复并不好，所以前面要运行主成分分析，看看最终选哪几组进行下游分析

top10_heatmap.pdf

进行GO注释和KEGG注释

#使用clusterProfiler包
library(clusterProfiler)
library(AnnotationDbi)
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db)

genelist <- res1_down$gene_id
genelist_2 <-res1_up$gene_id
gene_down<-gsub("\\..*","",genelist)
gene_up <-gsub("\\..*","",genelist_2)
#ENSMUSG是ensembl ID格式,需要格式转换，enrichGO函数对输入的基因名格式有要求，一般使用"ENTREZID”
#bitr 函数将基因ID从ENSEMBL类型转换为ENTREZID类型，指定使用的小鼠基因数据库
ids <- bitr(gene_down,fromType = "ENSEMBL",toType = "ENTREZID",OrgDb=org.Mm.eg.db,drop=TRUE)
ids_2 <- bitr(gene_up,fromType = "ENSEMBL",toType = "ENTREZID",OrgDb=org.Mm.eg.db,drop=TRUE)
#enrichGO 函数进行基因本体富集分析。ids$ENTREZID 是基因列表，指定基因ID类型为ENTREZID，ont = "ALL" 对所有GO类别进行富集分析，只有p值和q值小于0.05的结果会被保留，readable = T 将结果中的基因ID转换为可读的基因名称
GO_down <- enrichGO(ids$ENTREZID,OrgDb = org.Mm.eg.db,keyType = "ENTREZID",ont = "ALL",pvalueCutoff = 0.05,qvalueCutoff = 0.05,readable = T)
GO_up <- enrichGO(ids_2$ENTREZID,OrgDb = org.Mm.eg.db,keyType = "ENTREZID",ont = "ALL",pvalueCutoff = 0.05,qvalueCutoff = 0.05,readable = T)
GO_df_down<- as.data.frame(GO_down)
GO_df_up<- as.data.frame(GO_up)
View(GO_df_down)
View(GO_df_up)
GO_df_sorted_down <- GO_df_down[order(GO_df_down$pvalue),]
GO_df_sorted_up <- GO_df_up[order(GO_df_up$pvalue),]
#enrichKEGG函数进行KEGG分析
#富集的柱状图或者气泡图一般会选择显著性高的术语或通路进行作图。y轴一般有几种情况，如count/Gene Ratio/Fold enrichment。越大说明该条目下基因数越多、通路可能越重要。
#####但是需要注意的是，一般ID转换仅能成功转换30%~60%左右。很多基因没有被成功转换id，所以它们不能进行富集分析，这样会损失很多基因。所以可以多多尝试其他的基因富集分析，如使用DAVID网站。

使用DAVID网站进行功能富集分析，可以使用多种gene id进行分析
进入functional annotation tool，提交list

上传gene id

选择背景物种，选择注释类型，GO注释中BP、CC、MF。或者选择pathways中KEGG pathway进行KEGG注释。选择下面的Functional Annotation Chart

GO注释

chart

###DAVID富集分析后将ensembl_gene_id转为对应的gene_name
library(tidyverse)
GO_up <- read.csv("./GO_up.csv")
View(res1_up)
res1_up$gene_name <- rownames(res1_up)
##separate_rows 函数将 Genes 列中的基因ID按照逗号和空格拆分成多行。left_join 函数将拆分后的数据与 res1_up 数据框合并，以获取对应的 gene_name
GO_up_expanded <- GO_up %>% separate_rows(Genes,sep = ", ") %>% left_join(res1_Enew_up,by=c("Genes"="gene_id"))
View(GO_up_expanded)
##group_by 按 Term 列分组，然后使用 summarise 将每个组中的 gene_name 合并为一个字符串，并以斜杠分隔。summarise函数用于对每个分组执行汇总操作。它添加了genes新列
GO_up_combine <-GO_up_expanded %>% group_by(Term) %>% summarise(genes=paste(gene_name,collapse = "/"))
View(GO_up_combine)
#left_join 函数将合并后的数据与原始数据 Enew_GO_down 合并，以保留原始的 Term 列和其他信息。
GO_up_final <- Enew_GO_up_combine %>% left_join(GO_up,by=c("Term"="Term"))
View(GO_up_final)

res1_up

GO_up

GO_up_final

绘制GO注释和KEGG注释的柱状图或气泡图

#将前面转换好的表格保存，在excel中去除多余的列，重新读取
down_pathway_DAVID <- readxl::read_xlsx("./down_kegg.xlsx")
down_GO_DAVID <- readxl::read_xlsx("./down_go.xlsx")
View(down_pathway_DAVID)
View(down_GO_DAVID)
#根据FDR(多重检验矫正过的pvalue)排序。对于1个基因来说错误发生率为0.05，那么100个基因时，就可能有5个基因错误认为成差异基因，我们测序得到大约3万个基因，那就有1500个基因可能被错误认为成差异基因。多重检验矫正的目的就是为了减少这种假阳性
down_pathway_DAVID <- arrange(down_pathway_DAVID,down_pathway_DAVID[,8])
down_GO_DAVID <- arrange(down_GO_DAVID,down_GO_DAVID[,8])

#将description转换为因子类型，否则类型可能不统一，比如有的是字符串
down_pathway_DAVID$Description <- factor(down_pathway_DAVID$Description)
down_GO_DAVID$Description <- factor(down_GO_DAVID$Description)

#数据框已经排序，并且Description已经转换为因子，但为了确保ggplot2使用的是正确的因子顺序，我们仍需显式地设置因子的排列顺序。排序数据框并将列转换为因子并不能保证因子的级别按照排序顺序排列。
down_pathway_DAVID$Description <- factor(down_pathway_DAVID$Description, 
                                         levels = down_pathway_DAVID$Description[order(down_pathway_DAVID$`Fold Enrichment`)])
View(down_pathway_DAVID)
#气泡图可以展示4个变量，横坐标是富集倍数，富集倍数高意味着输入基因集在特定条目中的基因数相对背景基因集的比例更大，表明该条目在输入基因集中具有更高的显著性或重要性。纵坐标是通路名称，颜色根据显著性来定义，气泡大小是下调基因中属于该通路的基因数目
p_down_KEGG <- ggplot(down_pathway_DAVID,aes(x=down_pathway_DAVID$`Fold Enrichment`,y=Description,colour=-1*log10(PValue),size=Count))+
  geom_point()+
  scale_size(range=c(2, 8))+
  #颜色梯度：显著性由低到高为由蓝到红
  scale_colour_gradient(low = "blue",high = "red")+
  # 设置图形主题为白底
  theme_bw()+
  ylab("KEGG Pathway Terms")+
  xlab("Fold Enrichment")+
  #指定了颜色轴的标签。使用expression 函数来创建一个数学表达式，-log[10](PValue)
  labs(color=expression(-log[10](PValue)))+
  #设置图例标题和图例文本的大小为14
  theme(legend.title=element_text(size=14), legend.text = element_text(size=14))+
  #y轴标题的右侧边距为50个像素，x轴标题的顶部边距为20个像素
  theme(axis.title.y = element_text(margin = margin(r = 50)),axis.title.x = element_text(margin = margin(t = 20)))+
 #设置x轴文本的样式。face = "bold" 将文本设置为粗体。color = "black" 将文本颜色设置为黑色。angle = 0 将文本角度设置为0度（水平显示）。vjust = 1 垂直调整设置为1，通常指定文本在其分配空间中的垂直位置。
  theme(axis.text.x = element_text(face ="bold",color="black",angle=0,vjust=1))
p_down_KEGG
ggsave("down_KEGG_bubble.pdf",plot=p_down_KEGG,width = 13,height = 11) 

#GO柱状图
subset_down_GO <-subset(down_GO_DAVID,down_GO_DAVID$FDR < 0.05)
bp=subset_down_GO[which(subset_down_GO$Category %in% "BP"),]
bp=arrange(bp,bp[,11])
bp$Description <- factor(bp$Description, levels = bp$Description)
cc=subset_down_GO[which(subset_down_GO$Category %in% "CC"),]
View(cc)
cc=arrange(cc,cc[,11])
cc$Description <- factor(cc$Description, levels =cc$Description)
mf=subset_down_GO[which(subset_down_GO$Category %in% "MF"),]
mf=arrange(mf,mf[,11])
mf$Description <- factor(mf$Description, levels = mf$Description)
go_enrich_down=rbind(bp,cc,mf)
View(go_enrich_down)

go_enrich_down$Description <- factor(go_enrich_down$Description,levels = unique(go_enrich_down$Description))

#得到的图每一部分都没有按FDR的大小进行排序，问题可能出在 factor 的顺序设置上。需要确保 Description 列在每个部分（BP、CC、MF）中都根据 FDR 排序，然后再将它们合并为一个数据框
p_test <- ggplot(go_enrich_down,aes(y=Count,x=Description, fill=-log10(FDR))) +
  geom_bar(stat="identity",position = "dodge") +
  facet_grid(Category~.,scales = "free",space = "free")+
  coord_flip() + 
  theme_bw() +
  xlab("GO Terms")+
  ylab("Fold Enrichment")+
  labs(fill=expression(-log[10](qvalue)))+
  theme(
    plot.title = element_text(hjust = 0.5),
    strip.text.y = element_text(size = 14),
    legend.position="right",
    legend.title = element_text(size=15),
    legend.text = element_text(size=14),
    axis.text.x = element_text(size=10),
    axis.text.y = element_text(size=10),
    axis.title.x = element_text(size=14),
    axis.title.y = element_text(size=14)
  )+
  scale_fill_gradient(low = "blue", high = "red")
p_test
ggsave("down_GO_bar.pdf",plot=p_test_E45,width =14,height = 12)

down_KEGG_bubble.pdf

down_GO_bar.pdf