复现NC(三)| DEseq2+GO/KEGG富集分析

本系列为 RNA-Seq 分析具体实现过程，将具体描述流程及可执行代码～
本次复现文章来自 Nature communication，
题目为：Temporal and Spatial Heterogeneity of Host Response to SARS-CoV-2 Pulmonary Infection

本篇将介绍使用Deseq2进行差异表达分析和进行G0/KEGG富集分析的过程

Deseq2差异表达分析

得到所有测序结果的reads数后，需要使用Deseq2 进行差异表达分析：
其中涉及多个步骤，首先将对原始计数进行建模，使用归一化因子(大小因子)来考虑库深度的差异。然后，估算基因离散度，并缩小这些估计值，以生成更准确的离散度估计值，从而对计数进行建模。最后，DESeq2将拟合负二项模型，并使用Wald检验或似然比检验进行假设检验
具体在运行时，涉及以下程序：

#构建对象，coldata
dds <- DESeqDataSetFromMatrix(mycounts_data, coldata, design= ~ level)
dds <- DESeq(dds)
#计算，设置padj为0.05
res= results(dds,alpha = 0.01)
#查看离散度估计的曲线图来检查模型与我们的数据的匹配性
plotDispEsts(dds)

拟合图像 在p-value规定为0.05，FoldChange>2的情况下，可以得到5466个上调基因，93个下调基因 显著基因.gif

GO/KEGG富集分析

为了更情况的看清楚上一步筛选出的差异表达基因的功能，需要使用富集分析，功能富集分析是将基因或者蛋白列表分成多个部分，即将一堆基因进行分类，将具有相似功能的基因放到一起，并和生物学表型关联起来：

• GO分析：
  Cellular component，CC 细胞成分，解释基因位置
  Biological process， BP 生物学过程，生物学过程是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制
  Molecular function，MF 分子功能，例如催化某种反应

library(org.Hs.eg.db)
gene.diff_name<-bitr(diff_name$ensembl_gene_id, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = "org.Hs.eg.db")
#以cc为例
ego_cc_name<-enrichGO(gene       = gene.diff_name$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
barplot(ego_cc_name)

柱状图展示.png

• KEGG Pathway富集分析:
  KEGG是一个整合了基因组、化学和系统功能信息的综合数据库，我们使用其pathway部分

library(stringr)
kk<-enrichKEGG(gene = as.numeric(gene.diff_name$ENTREZID),
               organism = 'hsa',
               pvalueCutoff = 0.9,
               qvalueCutoff =0.9)
head(kk)
kk[1:30]
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

下图框中即为pathway的具体描述 kegg结果.png
气泡图展示