need:表達(dá)矩陣 分組矩陣，值要是整數(shù)
DESeq2和EdgeR都可用于做基因差異表達(dá)分析，主要也是用于RNA-Seq數(shù)據(jù)推盛，同樣也可以處理類似的ChIP-Seq,shRNA以及質(zhì)譜數(shù)據(jù)乃正。
這兩個(gè)都屬于R包，其相同點(diǎn)在于都是對(duì)count data數(shù)據(jù)進(jìn)行處理牢屋，都是基于負(fù)二項(xiàng)分布模型掰邢。因此會(huì)發(fā)現(xiàn)，用兩者處理同一組數(shù)據(jù)伟阔，最后在相同閾值下篩選出的大部分基因都是一樣的辣之，但是有一部分不同應(yīng)該是由于其估計(jì)離散度的不同方法所導(dǎo)致的。
library(DESeq2)
library(limma)
library(pasilla)#示例數(shù)據(jù)包
data(pasillaGenes)
exprSet=counts(pasillaGenes)? ##做好表達(dá)矩陣
group_list=pasillaGenes$condition##做好分組因子即可
(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list))
dds <- DESeqDataSetFromMatrix(countData = exprSet,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? colData = colData,#Rows of colData correspond to columns of countData
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? design = ~ group_list)
##上面是第一步第一步皱炉，構(gòu)建dds這個(gè)對(duì)象怀估，需要一個(gè)表達(dá)矩陣和分組矩陣！：辖痢多搀！
dds2 <- DESeq(dds)? ##第二步，直接用DESeq函數(shù)即可
resultsNames(dds2)#results extracts a result table from a DESeq analysis giving base means across samples,? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? log2 fold changes, standard errors, test statistics, p-values and adjusted p-values;? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? resultsNames returns the names of the estimated effects (coefficents) of the model;? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? removeResults returns a DESeqDataSet object with results columns removed.
res <-? results(dds2, contrast=c("group_list","treated","untreated"))#this argument specifies what comparison to extract from the object to build a results table
## 提取你想要的差異分析結(jié)果灾部，我們這里是treated組對(duì)untreated組進(jìn)行比較
resOrdered <- res[order(res$padj),]#padj排序
resOrdered=as.data.frame(resOrdered)
可以看到程序非常好用康铭！
它只對(duì)RNA-seq的基因的reads的counts數(shù)進(jìn)行分析，請(qǐng)不要用RPKM等經(jīng)過(guò)了normlization的表達(dá)矩陣來(lái)分析赌髓。
值得一提的是DESeq2軟件獨(dú)有的normlization方法从藤！
rld <- rlogTransformation(dds2)? ## 得到經(jīng)過(guò)DESeq2軟件normlization的表達(dá)矩陣！
exprSet_new=assay(rld)
par(cex = 0.7)#指定符號(hào)的大小
n.sample=ncol(exprSet)#樣本數(shù)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
hist(exprSet)
hist(exprSet_new)
第一張圖也是箱型圖
看這個(gè)圖就知道了锁蠕，它把本來(lái)應(yīng)該是數(shù)據(jù)離散程度非常大的RNA-seq的基因的reads的counts矩陣經(jīng)過(guò)normlization后變成了類似于芯片表達(dá)數(shù)據(jù)的表達(dá)矩陣夷野，然后其實(shí)可以直接用T檢驗(yàn)來(lái)找差異基因了！
但是荣倾，如果你的分組不只是兩個(gè)悯搔，就復(fù)雜了，你需要再仔細(xì)研讀說(shuō)明書(shū)舌仍，甚至你可能需要咨詢實(shí)驗(yàn)設(shè)計(jì)人員或者統(tǒng)計(jì)人員妒貌！
