目前差异分析主要利用的包有：edgeR、Deseq2、limma; 本质就是对表达量矩阵做一个归一化，让不同组样本的表达量具有可比性，然后利用理想的统计分布检验函数来计算差异的显著性。
limma
# 准备三个数据：表达矩阵、分组矩阵、差异比较矩阵
# 运行三个步骤：lmFit、eBayes、topTable
# 探索sCLLex数据
suppressPackageStartupMessages(library(CLL))
data(sCLLex) # 加载数据
exprSet=exprs(sCLLex) # 提取表达矩阵
samples=sampleNames(sCLLex) # 提取样本名
pdata=pData(sCLLex) # 分组信息
group_list=as.character(pdata[,2]) # 提取分组信息
dim(exprSet)
# 简单QC检测
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main="expression value",las=2)
# 制作分组矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# limma进行差异分析
## step1
fit <- lmFit(exprSet,design)
## step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
## step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
# write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
limma voom
limma包之前只能用于芯片数据，但是已经开发了voom方法来分析RNA-seq数据。
# 利用pasillaGenes数据
## 加载数据，获得表达矩阵
data(pasillaGenes)
exprSet=counts(pasillaGenes)
head(exprSet)
group_list=pasillaGenes$condition
# 构建分组矩阵
suppressMessages(library(limma))
design <- model.matrix(~factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# 根据分组信息和表达矩阵进行normalization
v <- voom(exprSet,design,normalize="quantile")
# 差异分析
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
tempOutput = topTable(fit2, coef=2, n=Inf)
DEG_voom = na.omit(tempOutput)
head(DEG_voom)
Deseq2
主要分三步：构建dds对象、进行normalization、使用reaults函数提取差异结果 简单版：
data(pasillaGenes)
exprSet=counts(pasillaGenes)
head(exprSet)
group_list=pasillaGenes$condition
# 构建dds对象
colData <- data.frame(row.names=colnames(exprSet),
group_list=group_list
)
## 这是一个复杂的方法构造这个对象！
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
## design 其实也是一个对象，还可以通过design(dds)来继续修改分组信息，但是一般没有必要。
dds
# 运行normalization
suppressMessages(dds2 <- DESeq(dds))
# 提取差异分析结果
resultsNames(dds2)
res <- results(dds2, contrast=c("group_list","treated","untreated")
resOrdered <- res[order(res$padj),]
resOrdered=as.data.frame(resOrdered)
head(resOrdered)
复杂版：
library(DESeq)
#加载数据
K1=read.table("742KO1.count",row.names=1)
K2=read.table("743KO2.count",row.names=1)
W1=read.table("740WT1.count",row.names=1)
W2=read.table("741WT2.count",row.names=1)
#列名
data=cbind(K1,K2,W1,W2)
#如果是htseq的结果，则删除data最后四行
n=nrow(data)
data=data
[c language="(-n+4:-n),"][/c]
#如果是bedtools的结果，取出统计个数列和行名
kk1=cbind(K1$V5)
rownames(kk1)=rownames(K1)
K1=kk1
#差异分析
colnames(data)=c("K1","K2","W1","W2")
type=rep(c("K","W"),c(2,2))
de=newCountDataSet(data,type)
de=estimateSizeFactors(de)
de=estimateDispersions(de)
res=nbinomTest(de,"K","W")
#res就是我们的表达量检验结果
到这里，理论上差异基因的分析已经结束啦！后面只是关于R的bioconductor包的一些简单结合使用而已
library(org.Mm.eg.db)
tmp=select(org.Mm.eg.db, keys=res$id, columns=c("ENTREZID","SYMBOL"), keytype="ENSEMBL")
#合并res和tmp
res=merge(tmp,res,by.x="ENSEMBL",by.y="id",all=TRUE)
#go
tmp=select(org.Mm.eg.db, keys=res$ENSEMBL, columns="GO", keytype="ENSEMBL")
ensembl_go=unlist(tapply(tmp[,2],as.factor(tmp[,1]),function(x) paste(x,collapse ="|"),simplify =F))
#为res加入go注释，
res$go=ensembl_go[res$ENSEMBL]#为res加入一列go
#写入all——data
all_res=res
write.csv(res,file="all_data.csv",row.names =F)
uniq=na.omit(res)#删除无效基因
sort_uniq=uniq[order(uniq$padj),]#按照矫正p值排序
#写入排序后的all_data
write.csv(res,file="all_data.csv",row.names =F)
#标记上下调基因
sort_uniq$up_down=ifelse(sort_uniq$baseMeanA>sort_uniq$baseMeanB,"up","down")
#交换上下调基因列位置
final_res=sort_uniq[,c(12,1:11)]
#写出最后数据
write.csv(final_res,file="final_annotation_gene_bedtools_sort_uniq.csv",row.names =F)
#然后挑选出padj值小于0.05的数据来做富集
tmp=select(org.Mm.eg.db, keys=sort_uniq[sort_uniq$padj<0.05,1], columns="ENTREZID", keytype="ENSEMBL")
diff_ENTREZID=tmp$ENTREZID
require(DOSE)
require(clusterProfiler)
diff_ENTREZID=na.omit(diff_ENTREZID)
ego <- enrichGO(gene=diff_ENTREZID,organism="mouse",ont="CC",pvalueCutoff=0.05,readable=TRUE)
ekk <- enrichKEGG(gene=diff_ENTREZID,organism="mouse",pvalueCutoff=0.05,readable=TRUE)
write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(summary(ego),"GO-enrich.csv",row.names =F)
