本系列为 RNA-Seq 分析具体实现过程，将具体描述流程及可执行代码～
本次复现文章来自 Nature communication，
题目为：Temporal and Spatial Heterogeneity of Host Response to SARS-CoV-2 Pulmonary Infection
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)
拟合图像
显著基因.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))
kegg结果.png
气泡图展示
