今天更新TCGA數(shù)據(jù)庫的利用系列第三篇文章，在對TCGA數(shù)據(jù)進(jìn)行挖掘時(shí)，通常會篩選出來一些表達(dá)量顯著異常的基因，作為后續(xù)研究的對象，這個(gè)篩選過程叫做差異分析；本篇文章將分為三大模塊對差異分析進(jìn)行介紹。
關(guān)于差異分析的官方解釋：
差異分析就是將一組資料的總變動(dòng)量，以可能造成變動(dòng)的因素分解成不同的部分，并且以假設(shè)鑒定的方法來判斷這些因素是否確實(shí)能解釋資料的變動(dòng)。
我自己的一點(diǎn)理解：差異分析就是對總體樣本數(shù)據(jù)中非正常數(shù)據(jù)的篩選。對于轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異分析時(shí)，limma 、edgeR、DESeq2這三種程序包都可以（limma相對較受歡迎），大部分科研性文章基本上是用其中一種方法篩選差異表達(dá)基因，但為了使得結(jié)果更加準(zhǔn)確，在做畢業(yè)課題時(shí)我把三種方法都做了一遍，把它們結(jié)果的交集作為篩選出來的差異表達(dá)基因。
不管用那一種方法做差異分析，基本上要做的步驟就是：一，創(chuàng)建表達(dá)矩陣；二、創(chuàng)建設(shè)計(jì)矩陣（分組矩陣）；三、得到差異表達(dá)分析矩陣。但不同包基于算法、數(shù)據(jù)模型不同，所用的函數(shù)、篩選標(biāo)準(zhǔn)也大不相同，所以用代碼實(shí)現(xiàn)時(shí)結(jié)果有很大的差別。
無論用那種包做差異分析，在做之前必須要保證需要用到的包已經(jīng)安裝成功。在R語言中安裝程序包的代碼（其中的一種方式）如下：
source('https:///biocLite.R')#加載的網(wǎng)址都一樣
biocLite('edgeR')#把雙引號的內(nèi)容換成你所需要程序包的名稱即可
limma做差異分析
傳入原始樣本基因表達(dá)矩陣（表達(dá)矩陣格式如下圖）
接下來就是對基因表達(dá)矩陣進(jìn)行一些處理，讓樣本名變成數(shù)據(jù)矩陣的列名，基因名變成數(shù)據(jù)矩陣的行名，同時(shí)把ensembl_symbl那一列去掉(用express_rec <- express_rec[,(-1)]命令即可)，變成如下這個(gè)格式：
表達(dá)矩陣?yán)锩娴臄?shù)據(jù)太大，但為了使數(shù)據(jù)呈現(xiàn)正態(tài)分布，需要對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，這里我用的是函數(shù)log(express_rec,2)，在標(biāo)準(zhǔn)化之前，需要把表達(dá)矩陣內(nèi)為0的數(shù)據(jù)賦值為1，目的是為了防止取log時(shí)，數(shù)據(jù)變?yōu)樨?fù)無窮大。
（express_rec[express_rec==0<-1）
下面進(jìn)行分組矩陣的組建，首先提前創(chuàng)建好一個(gè)矩陣列表，如下，行名為樣本編碼，列名為樣本類型，如下面這種格式：
而limma包用到的設(shè)計(jì)矩陣是下面這種格式：
實(shí)現(xiàn)代碼如下：
rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]
Group <- factor(group_text$group,levels = c('Tumor','Normal'))
design <- model.matrix(~0+Group)
colnames(design) <- c('Tumor','Normal')
rownames(design) <- rownames(group_text)
接下來的步驟依次進(jìn)行數(shù)據(jù)擬合、經(jīng)驗(yàn)貝葉斯檢驗(yàn)、篩選差異表達(dá)基因。
fit <- lmFit(express_rec,design)
#制作比對標(biāo)準(zhǔn)；
contrast.matrix <- makeContrasts(Tumor - Normal,levels=design)
fit2 <- contrasts.fit(fit,contrast.matrix)
#進(jìn)行經(jīng)驗(yàn)貝葉斯檢驗(yàn)；
fit2 <- eBayes(fit2)
#基于logFC為標(biāo)準(zhǔn)，設(shè)置數(shù)量上限為30000，調(diào)整方法為fdr;
all_diff <- topTable(fit2, adjust.method = 'fdr',coef=1,p.value = 1,lfc <- log(1,2),number = 30000,sort.by = 'logFC')#從高到低排名；
limma包的另一種方法，精確權(quán)重法（voom）,然后把篩選得到的差異表達(dá)基因?qū)懭隿sv文件中。
#limma的另一種方法；
dge <- DGEList(counts = express_rec)
dge <- calcNormFactors(dge)#表達(dá)矩陣進(jìn)行標(biāo)準(zhǔn)化；
v <- voom(dge, design,plot=TRUE)#利用limma_voom方法進(jìn)行差異分析；
fit <- lmFit(v, design)#對數(shù)據(jù)進(jìn)行線性擬合；
fit <- eBayes(fit)#貝葉斯算法組建
all <- topTable(fit, coef=ncol(design),n=Inf)#從高到低排名；
sig.limma <- subset(all_diff,abs(all$logFC)>1.5&all$P.Value<0.05)#進(jìn)行差異基因篩選；
write.csv(sig.limma,'C:/Users/FREEDOM/Desktop/TCGA_data/limm_diff.csv')#寫入csv文件中；
DESeq2做差異分析
第一步跟limma程序包一樣，讀入表達(dá)矩陣，對表達(dá)矩陣進(jìn)行數(shù)據(jù)處理；
express_rec<- read.csv('C:/Users/FREEDOM/Desktop/TCGA_data/after_note2.csv')#讀取數(shù)據(jù)
group_text <- read.csv('C:/Users/FREEDOM/Desktop/TCGA_data/group_text.csv')
library('DESeq2')#加載包；
install.packages('rpart')#有這個(gè)包可忽略，沒有的時(shí)候才安裝；
express_rec <- express_rec[,-1]
rownames(express_rec) <-express_rec[,1]
express_rec <- express_rec[(-1)]#表達(dá)矩陣的處理；
rownames(group_text) <- group_text[,1]
創(chuàng)建分組矩陣；
rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]#分組矩陣的數(shù)據(jù)處理；
all(rownames(group_text)==colnames(express_rec))#確保表達(dá)矩陣的列名與分組矩陣行名相一致；
構(gòu)建 DESeq2 所需的 DESeqDataSet 對象；
dds <- DESeqDataSetFromMatrix(countData=express_rec, colData=group_text, design<- ~ group) #DESeq2的加載
head(dds)
dds <- dds[rowSums(counts(dds)) > 1, ] #過濾一些low count的數(shù)據(jù)；
使用DESeq進(jìn)行差異表達(dá)分析，返回 results可用的DESeqDataSet對象；
> dds <- DESeq(dds)#DESeq進(jìn)行標(biāo)準(zhǔn)化；
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2819 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
可以用summary函數(shù)查看表達(dá)基因上下調(diào)分布基本情況；
> summary(res)#查看經(jīng)過標(biāo)準(zhǔn)化矩陣的基本情況；
out of 24823 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 7590, 31%
LFC < 0 (down) : 5120, 21%
outliers [1] : 0, 0%
low counts [2] : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
最后提取差異表達(dá)基因，存入csv文件中。
edgeR進(jìn)行差異分析
這個(gè)包的操作步驟與limma包類似，首先就是讀入數(shù)據(jù)，創(chuàng)建表達(dá)、分組矩陣。
path = 'C:/Users/FREEDOM/Desktop/TCGA_data/after_note2.csv'
path1 ='C:/Users/FREEDOM/Desktop/TCGA_data/group_text.csv'
express_rec <- read.csv(path,headers <- T)#讀取表達(dá)矩陣;
group_text <- read.csv(path1,headers <- T)#讀取分組矩陣；
library(edgeR)#加載edgeR包
express_rec <- express_rec[,-1]
rownames(express_rec) <- express_rec[,1]
express_rec <- express_rec[(-1)]#創(chuàng)建表達(dá)矩陣；
rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]#加載分組矩陣；
group <-factor(group_text$group)
dge <- DGEList(counts = express_rec,group = group)#構(gòu)建DEList對象；
y <- calcNormFactors(dge)#利用calcNormFactor函數(shù)對DEList對象進(jìn)行標(biāo)準(zhǔn)化(TMM算法)
#創(chuàng)建設(shè)計(jì)矩陣，跟Limma包相似；
Group <- factor(group_text$group,levels = c('Tumor','Normal'))
design <- model.matrix(~0+Group)
colnames(design) <- c('Tumor','Normal')
rownames(design) <- rownames(group_text)#創(chuàng)建分組矩陣；
edgeR包創(chuàng)建的分組矩陣與limma一樣，是以factor因子格式展現(xiàn)出來；
接下來依次進(jìn)行構(gòu)建DGEList對象、利用TMM算法對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化、估計(jì)離散值、數(shù)據(jù)擬合、創(chuàng)建對比矩陣、對數(shù)據(jù)做QL-text檢驗(yàn)、差異表達(dá)基因?qū)懭隿sv文件中。
dge <- DGEList(counts = express_rec,group = group)#構(gòu)建DEList對象；
y <- calcNormFactors(dge)#利用calcNormFactor函數(shù)對DEList對象進(jìn)行標(biāo)準(zhǔn)化(TMM算法)
y <- estimateDisp(y,design)#估計(jì)離散值（Dispersion）
fit <- glmQLFit(y, design, robust=TRUE)#進(jìn)一步通過quasi-likelihood (QL)擬合NB模型
TU.vs.NO <- makeContrasts(Tumor-Normal, levels=design)#這一步主要構(gòu)建比較矩陣；
res <- glmQLFTest(fit, contrast=TU.vs.NO)#用QL F-test進(jìn)行檢驗(yàn)
# ig.edger <- res$table[p.adjust(res$table$PValue, method = 'BH') < 0.01, ]#利用‘BH’方法；
result_diff <- res$table#取出最終的差異基因；
write.csv(edge_diff,'C:/Users/FREEDOM/Desktop/TCGA_data/edgeR_diff2.csv')
