一文解决RNA测序资料的差异分析（limma,deseq,edg
本文目标：
（1）使用edger包做TCGA数据库RNA-seq数据差异分析
（2）使用deseq包做TCGA数据库RNA-seq数据差异分析
（3）使用limma包做TCGA数据库RNA-seq数据差异分析
（4）如何在没有生物学重复的情况下（比如说只有两个样本，来求取差异基因）
DESeq2和EdgeR都可用于做基因差异表达分析，主要也是用于RNA-Seq数据，同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。
这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。
edgeR 使用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因。 特别地，经验贝叶斯用于通过在基因之间来调节跨基因的过度离散程度。 使用类似于Fisher精确检验但适应过度分散数据的精确检验用于评估每个基因的差异表达。edgeR 在默认情况下，执行TMM标准化程序以考虑样本之间的不同测序深度，通过Benjamini-Hochberg用于控制FDR 。
Limma包基于线性模型建模。 它最初设计用于分析微阵列数据，但最近已扩展到RNA-seq数据。 根据limma用户指南的当前建议是使用edgeR包的TMM标准化和“voom”转换，其本质上将标准化数据取对数（基数2）并估计它们的均值 - 方差关系以确定在线性建模之前每次观察的权重。 默认情况下，Benjamini-Hochberg程序用于估计FDR 。
DESeq使用类似于edgeR的负二项式模型，与edgeR类似，执行缩放因子归一化以考虑不同样本的变化的测序深度，并且Benjamini-Hochberg用于控制FDR。 DESeq能够分析具有少量重复的实验。DESeq技术上可以在没有任何生物学重复的情况下进行实验。DESeq2是在DESeq基础上更新的软件。
（1）edgeR包的差异分析代码。
rm(list = ls())
library("DESeq")
library("limma")
library("edgeR")
foldChange=1
padj=0.05
setwd("D:\\train\\diff")
#读取数据，表达量数据每一列是一个样本，第一列为基因名
rt = read.csv("mRNA_exprSet.csv",sep=",",header=T)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
#删除掉表达量较低的基因
data=data[rowMeans(data)>1,]
#构建分组，前60列为正常样本，后63列为肿瘤样本
group =c(rep("normal",60),rep("tumor",63))
design = factor(group)
design <- model.matrix(~group)
y <- DGEList(counts=data,group=group)
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)
et <- exactTest(y,pair = c("normal","tumor"))
topTags(et)
ordered_tags <- topTags(et, n=100000)
allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
diff=allDiff
newData=y$pseudo.counts
# write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
write.csv(diffSig, file="edger_diffSig.csv")
#
# diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]
# write.table(diffUp, file="edger_up.xls",sep="\t",quote=F)
# diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
# write.table(diffDown, file="edger_down.xls",sep="\t",quote=F)
#
