前言
转录组入门(7)： 差异基因分析这个步骤推荐在R里面做，载入表达矩阵，然后设置好分组信息，统一用DEseq2进行差异分析，当然也可以走走edgeR或者limma的voom流程。基本任务是得到差异分析结果，进阶任务是比较多个差异分析结果的异同点。生信技能树
实验操作
rm(list=ls())
library(DESeq2)
setwd("E:/0ngs/sam_out/")
expr <- read.table("Akap95_expr.txt",header = T,stringsAsFactors = F)
rownames(expr) <- expr$Gene
expr <- expr[,-1]
表达矩阵expr如下：
1504003725778.png
# library(edgeR)
# keepGene <- rowSums(cpm(expr)>0) >=2
# table(keepGene); dim(expr)
group=c('con','con','treat','treat')
colData <- data.frame(row.names=colnames(expr), group=group)
dds <- DESeqDataSetFromMatrix(countData = expr,
colData = colData,
design = ~ group)
dds <- DESeq(dds)
png("qc_dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds, main="Dispersion plot") # 散度估计
dev.off()
# 提取差异基因
res <- results(dds, contrast = c('group','treat','con'))
resOrdered <- res[order(res$padj),]
DEG <- as.data.frame(resOrdered)
write.csv(DEG,"Akap95-DEG.csv")
差异表达显著的基因：
1504004764836.png
# normalizedCounts1与normalizedCounts2相同
normalizedCounts1 <- t(t(counts(dds)) / sizeFactors(dds))
exprMatrix_rpm=as.data.frame(normalizedCounts1)
normalizedCounts2 <- counts(dds, normalized=T)
# exprMatrix_rpm[grep("Akap8",rownames(exprMatrix_rpm)),]
# 单独用DESeq2的normlization
rld <- rlogTransformation(dds)
exprMatrix_rlog=assay(rld)
# 原始表达值与DESeq2的normlization后的表达值比较
png("DEseq_RAWvsNORM.png",height = 800,width = 800)
par(cex = 0.7)
n.sample=ncol(expr)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(expr, col = cols,main="expression value",las=2)
boxplot(exprMatrix_rlog, col = cols,main="expression value",las=2)
hist(as.matrix(expr))
hist(exprMatrix_rlog)
dev.off()
PS: 本篇都是用的建明群主的代码里DESeq2的部分完成的。
