以Macaque样本为例进行Bulk/Smart-seq2数据一般分析流程的演示（二）——使用DESeq2获得组间差异基因以及KEGG通路富集分析

设置工作路径，清空工作环境

setwd("/data1/guolhh/proj/monkeysNPC/")
rm(list = ls())

读取上游处理好的counts数据，并调整表达矩阵列名等项目

counts <- fread("./counts/featureCount/raw_counts_Macaque.csv",header = T)#[,-c(2:6)]
write.xlsx(counts, "test.xlsx")
counts <- read.xlsx("./test.xlsx",sheet = 1)#我的Rstudio不能正确读取csv，正常情况下这两步可省略
library(stringr)#调用需要的包
#整理列名得到表达矩阵
colnames(counts)[2:19] <- str_split_fixed(colnames(counts)[2:19], "[_]", 2)[,1]
colnames(counts)[2:19] <- str_split_fixed(colnames(counts)[2:19], "[/]", 2)[,2]
colnames(counts)[8] <- "M12YB"
colnames(counts)[4] <- "E152B"
colnames(counts)[9:10] <- c("M17YP1B","M17YP1L")
colnames(counts)[14:15] <- c("M11YP1B","M11YP1L")

根据不同需求，对表达矩阵进行分组，进行DESeq2并储存相应结果

age_group <- counts[,c(1:3,5:7,11:13,16,19)]

#调整表达矩阵列的顺序，使其按照年龄顺序排列
col <- colnames(age_group)
new_col <- c(col[1],col[5:6],col[3:4],col[2],col[11],col[8:10],col[7])
new_col
age_group <- age_group[,new_col]
age_group <- age_group[,-c(2:4)]
age_vs_newborn <- age_group[,c(1:6)]

#DEseq2需要的表达矩阵需要是纯数字，所以将基因信息转到表达矩阵行名中
df <- age_vs_newborn[,-1]
rownames(df) <- as.character(age_vs_newborn$Geneid)
#对counts进行筛选，全样本counts和大于20并且都不为0被留下
count <- df[rowSums(df)>20&apply(df,1,function(x){all(x>0)}),]

#调用DESeq2
library(DESeq2)
#根据样本重复情况设置metadata
metadata = data.frame(
  dex = factor(rep(c("new_born","adult"), c(2,3))),
  id = factor(colnames(age_vs_newborn[2:6]))
)

#DESeq2主流流程
dds <- DESeqDataSetFromMatrix(countData=count, 
                              colData=metadata, 
                              design=~dex, 
)
dds$dex <- relevel(dds$dex, ref = "new_born")#规定对照组
dds$dex
colData(dds)
dds <- DESeq(dds)
#查看DESeq结果
resultsNames(dds) 
#[1] "Intercept"        "dex_adult_vs_new_born"
#提取DESeq运算结果
res <- results(dds,tidy = T)
#挑选表达矩阵中的基因进行test，以确保ctrl和treat组没有设置反
list <- c("ENSMMUG00000000634","ENSMMUG00000000632")
#查看随机挑选的基因在counts中和在res中的趋势是否一致
df[rownames(df) %in% list,]
res[res$row %in% list,]

#获取logFC > 1的基因为接下来的富集分析做准备
LFC1 <- res[abs(res$log2FoldChange) > 1&res$pvalue < 0.05,]
#dir.create("./result/Mmul10/aged_vs_young/DESeq2")
write.xlsx(res,"./result/Mmul10/aged_vs_young/DESeq2/DESeq(adult_vs_new_born)-result.xlsx")
write.xlsx(LFC1,"./result/Mmul10/aged_vs_young/DESeq2/LFC1_pV0.05(adult_vs_new_born)-result.xlsx")

接下来是根据差异基因进行KEGG通路富集

#读取保存好的差异基因
LFC1 <- read.xlsx("./result/Mmul10/aged_vs_young/DESeq2/LFC1_pV0.05(adult_vs_new_born)-result.xlsx",sheet = 1)
###KEGG
group <- "adult_vs_new_born"
#创建保存结果的文件夹
dir.create(paste("./result/Mmul10/aged_vs_young/DESeq2/",group,sep = ""))
df <- LFC1
num <- nrow(LFC1)
#调用需要的包
library("ggplot2")
library( "clusterProfiler")
library("org.Mmu.eg.db")
#load富集库，也可以使用线上的库
load("/data1/guolhh/GO_ref/enrichment_annotation/crab_eating_enrichment.RData")
columns(org.Mmu.eg.db)
#对基因名进行转换，使其和库中的基因名一致
geneinfo_K = AnnotationDbi::select(org.Mmu.eg.db, keys=keys(org.Mmu.eg.db), 
                                   columns = c('ENSEMBL',"SYMBOL"))
geneinfo_K <- unique(geneinfo_K)
df <- unique(df)
colnames(df)[1] <- "Gene.Name"
df <- merge(df, geneinfo_K, by.x = "Gene.Name", by.y = "ENSEMBL")
up_gene <- df[df$log2FoldChange > 0,]
#up
genename <- up_gene$SYMBOL
term2gene <- kegg[, c(1, 5)]
term2name <- kegg[, c(1, 2)]
ekegg_up <- enricher(gene = genename, pvalueCutoff = 0.05,
                     pAdjustMethod = "BH", TERM2GENE = term2gene,
                     TERM2NAME = term2name)
kegg_up <- as.data.frame(ekegg_up@result)
kegg_top20_up <- kegg_up[kegg_up$pvalue < 0.05,][1:30,]
kegg_top20_up$change <- 'up'

#down
down_gene <- df[df$log2FoldChange < 0,]
genename <- down_gene$SYMBOL
ekegg_down <- enricher(gene = genename, pvalueCutoff = 0.05,
                       pAdjustMethod = "BH", TERM2GENE = term2gene,
                       TERM2NAME = term2name)
kegg_down <- as.data.frame(ekegg_down@result)
kegg_top20_down <- kegg_down[order(kegg_down$pvalue),][1:10,]
kegg_top20_down$change <- 'down'
kegg_result <- rbind(kegg_top20_up, kegg_top20_down)
kegg_result$change <- factor(kegg_result$change, levels = c('up', 'down'))
kegg_result$Description <- factor(kegg_result$Description, levels = rev(unique(kegg_result$Description)))
kegg_result$change <- factor(kegg_result$change, levels = c('up', 'down'))

#富集结果图像绘制
p <- ggplot(kegg_result,aes(x = -log10(pvalue), y = Description,
                            color = -log10(pvalue), size = Count)) + geom_point() + 
  theme_bw() + scale_color_gradient(low = "blue",high = "red") +
  scale_size_continuous(range = c(1,5))+
  labs(y='',x='-log10(pvalue)',title='',color = "-log.p")+
  theme() +
  geom_vline(xintercept = -log10(0.05),linetype =2)+
  theme(plot.title=element_text(size=22,color="black",hjust = 0.5),
        axis.text=element_text(size=12, color="black"),
        axis.text.x=element_text(angle = 0,hjust = 0.5),
        axis.text.y = element_text(size = 18),
        axis.title=element_text(size=20),
        axis.line=element_line(size=1,colour='black'),
        axis.ticks=element_line(size=1,colour='black'),
        strip.text = element_text(size = 18),
        legend.text=element_text(size=16),legend.title = element_text(size=20))+
  facet_grid(~change ,scales = 'free')
#保存需要的结果（图片，基因列表等）
write.xlsx(kegg_result,file = paste("./result/Mmul10/aged_vs_young/DESeq2/",group,"/KEGG_",group,"_LFC1_",num," genes.xlsx",sep = ""))
ggsave(paste("./result/Mmul10/aged_vs_young/DESeq2/",group,"/KEGG_",group,"_LFC1_",num," genes.png",sep = ""), p, width = 16, height = 16)