前言

轉(zhuǎn)錄組入門(7)： 差異基因分析這個步驟推薦在R里面做，載入表達矩陣，然后設(shè)置好分組信息，統(tǒng)一用DEseq2進行差異分析，當然也可以走走edgeR或者limma的voom流程。基本任務(wù)是得到差異分析結(jié)果，進階任務(wù)是比較多個差異分析結(jié)果的異同點。生信技能樹

實驗操作

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的部分完成的。