RNASeq实战练习-DESeq2差异及clusterProfi

2021-08-24 本文已影响0人小小白的jotter

转录组接下来的分析就用 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及富集分析

GO、KEGG富集分析（一）有参情况

clusterProfiler基因功能富集分析 +气泡图

Bioconductor的镜像修改

转换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