RNA-seq(7) : DEseq2筛选差异表达基因并注释(bioMart)
转录组接下来的分析就用 R 进行
DESeq2差异分析
数据准备
比对得到的 countdata.csv 文件
image-20210806103139849
样本信息表，保存为 coldata.csv 文件
image-20210809233843177
载入数据及分析
# 安装 DESeq2 包
install.packages('BiocManager') #已安装就不需要再安装
BiocManager::install('DESeq2')
# 读取数据
countdata <- read.csv('countdata.csv', row.names = 1, header = T)
coldata <- read.csv('coldata.csv')
# 载入 DESeq2 包
library(DESeq2)
# 构建 DESeqDataSet 对象
dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design= ~condition)
# 计算差异倍数并获得 p 值
dds1 <- DESeq(dds, fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)
# 将 NSR 在前，WT 在后，意为 NSR 相较于 WT 中哪些基因上调/下调
res <- results(dds1, contrast = c('condition', 'NSR', 'WT'))
#输出表格至本地
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
write.table(res1, 'DESeq2.txt', col.names = NA, sep = '\t', quote = FALSE)
image-20210820102408384
筛选差异表达基因
# 按 padj 值升序排序，相同 padj 值下继续按 log2FC 降序排序
res1 <- res1[order(res1$padj, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
# 将 up，down，none 的基因筛选并标记出来
res1[which(res1$log2FoldChange >= 1 & res1$padj < 0.01),'sig'] <- 'up'
res1[which(res1$log2FoldChange <= -1 & res1$padj < 0.01),'sig'] <- 'down'
res1[which(abs(res1$log2FoldChange) <= 1 | res1$padj >= 0.01),'sig'] <- 'none'
#输出选择的差异基因总表
res1_diff <- subset(res1, sig %in% c('up', 'down'))
write.table(res1_diff, file = 'DESeq2.diff.txt', sep = '\t', col.names = NA, quote = FALSE)
#根据 up 和 down 分开输出
res1_up <- subset(res1, sig == 'up')
res1_down <- subset(res1, sig == 'down')
write.table(res1_up, file = 'DESeq2.up.txt', sep = '\t', col.names = NA, quote = FALSE)
write.table(res1_down, file = 'DESeq2.down.txt', sep = '\t', col.names = NA, quote = FALSE)
clusterProfiler转换ID及富集分析
转换ID
# 修改镜像
chooseBioCmirror()
image-20210820101849879
BiocManager::install("clusterProfiler")
library(clusterProfiler)
BiocManager::install("org.At.tair.db") # 镜像改一下安装很快
library(org.At.tair.db)
gene <- row.names(res1_diff)
columns(org.At.tair.db)
image-20210818235111858
# 转换 ID，bitr 会把有缺失的行删掉
tansid <- select(org.At.tair.db,keys = gene,columns = c("GENENAME","SYMBOL","ENTREZID"),keytype = "TAIR")
tansid1 <- bitr(gene,fromType = "TAIR",toType = c("GENENAME","SYMBOL","ENTREZID"),OrgDb = "org.At.tair.db")
write.table(tansid, file = 'IDS.txt', sep = '\t', col.names = NA, quote = FALSE)
write.table(tansid1, file = 'ID.bitr.txt', sep = '\t', col.names = NA, quote = FALSE)
GO分析
# 为了有结果，参数设置的有点高
go.all <- enrichGO(gene = tansid$ENTREZID,OrgDb = org.At.tair.db,keyType = 'ENTREZID',ont = 'ALL',pAdjustMethod = "BH",pvalueCutoff = 0.3,qvalueCutoff = 0.3)
#随后对富集结果进行总览，查看BP，CC，MF的个数
dim(go.all[go.all$ONTOLOGY=='BP',]);dim(go.all[go.all$ONTOLOGY=='CC',]);dim(go.all[go.all$ONTOLOGY=='MF',])
#保存结果
write.csv(go.all@result,'DEG_go.all.result.csv',row.names=F)
image-20210823131753615
KEGG分析
查找KEGG物种简写：https://www.genome.jp/kegg/catalog/org_list.html
拟南芥的物种简写是 ath
# 注意一下基因使用的 ID，我之前一直以为是 ENTREZID，结果一直报错
enrich.kegg <- enrichKEGG(gene =tansid$TAIR,
organism ="ath",
keyType = "kegg",
pvalueCutoff = 1,
pAdjustMethod = "BH",
qvalueCutoff = 1,
use_internal_data =FALSE)
dim(enrich.kegg)
write.csv(enrich.kegg@result,'DEG_KEGG.result.csv',row.names=F)
image-20210823132032952
