今天更新TCGA數(shù)據(jù)庫的利用系列第三篇文章姆怪，在對TCGA數(shù)據(jù)進行挖掘時，通常會篩選出來一些表達量顯著異常的基因，作為后續(xù)研究的對象稽揭，這個篩選過程叫做差異分析俺附；本篇文章將分為三大模塊對差異分析進行介紹
關(guān)于差異分析的官方解釋：
差異分析就是將一組資料的總變動量，依可能造成變動的因素分解成不同的部份溪掀，并且以假設(shè)檢定的方法來判斷這些因素是否確實能解釋資料的變動事镣。
我自己的一點理解：差異分析就是對總體樣本數(shù)據(jù)中非正常數(shù)據(jù)的篩選。對于轉(zhuǎn)錄組數(shù)據(jù)進行差異分析時揪胃，limma 璃哟、edgeR、DESeq2這三種程序包都可以（limma相對較受歡迎）喊递，大部分科研性文章基本上是用其中一種方法取篩選差異表達基因随闪，但為了使得結(jié)果更加準確，在做畢業(yè)課題時我把三種方法都做了一遍骚勘，把它們結(jié)果的交集作為篩選出來的差異表達基因铐伴。
不管用那一種方法做差異分析，基本上要做的步驟就是：一俏讹，創(chuàng)建表達矩陣当宴；二、創(chuàng)建設(shè)計矩陣（分組矩陣）泽疆；三户矢、得到差異表達分析矩陣。
但不同包基于算法于微、數(shù)據(jù)模型不同逗嫡，所用的函數(shù)、篩選標準也大不相同株依，所以用代碼實現(xiàn)時結(jié)果有很大的差別。
無論用那種包做差異分析延窜，在做之前必須要保證需要用到的包已經(jīng)安裝成功恋腕。在R語言中安裝程序包的代碼（其中的一種方式）如下：
source('https://bioconductor.org/biocLite.R')#加載的網(wǎng)址都一樣
biocLite("edgeR")#把雙引號的內(nèi)容換成你所需要程序包的名稱即可
limma做差異分析
傳入原始樣本基因表達矩陣（表達矩陣格式如下圖）
接下來就是對基因表達矩陣進行一些處理，讓樣本名變成數(shù)據(jù)矩陣的列名逆瑞，基因名變成數(shù)據(jù)矩陣的行名荠藤，同時把ensembl_symbl那一列去掉(用express_rec <- express_rec[,(-1)]命令即可)，變成如下這個格式：
表達矩陣里面的數(shù)據(jù)太大获高，但為了使數(shù)據(jù)呈現(xiàn)正態(tài)分布哈肖，需要對數(shù)據(jù)進行標準化，這里我用的是函數(shù)log(express_rec,2)念秧，在標準化之前淤井，需要把表達矩陣內(nèi)為0的數(shù)據(jù)賦值為1，目的是為了防止取log時，數(shù)據(jù)變?yōu)樨摕o窮大币狠。
（express_rec[express_rec==0<-1）
下面進行分組矩陣的組建游两，首先提前創(chuàng)建好一個矩陣列表，如下漩绵，行名為樣本編碼贱案，列名為樣本類型，如下面這種格式：
而limma包用到的設(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)
接下來的步驟依次進行數(shù)據(jù)擬合止吐、經(jīng)驗貝葉斯檢驗宝踪、篩選差異表達基因
fit?<-?lmFit(express_rec,design)
#制作比對標準；
contrast.matrix?<-?makeContrasts(Tumor?-?Normal,levels=design)
fit2?<-?contrasts.fit(fit,contrast.matrix)
#進行經(jīng)驗貝葉斯檢驗碍扔；
fit2?<-?eBayes(fit2)
#基于logFC為標準瘩燥，設(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）,然后把篩選得到的差異表達基因?qū)懭隿sv文件中。
#limma的另一種方法套鹅；
dge<-?DGEList(counts?=?express_rec)
dge?<-?calcNormFactors(dge)#表達矩陣進行標準化站蝠；
v?<-?voom(dge,?design,plot=TRUE)#利用limma_voom方法進行差異分析；
fit?<-?lmFit(v,?design)#對數(shù)據(jù)進行線性擬合卓鹿；
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)#進行差異基因篩選；
write.csv(sig.limma,'C:/Users/FREEDOM/Desktop/TCGA_data/limm_diff.csv')#寫入csv文件中吟孙；
DESeq2做差異分析
第一步跟limma程序包一樣澜倦，讀入表達矩陣，對表達矩陣進行數(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')#含有這個包可忽略藻治，沒有的時候才安裝；
express_rec?<-?express_rec[,-1]
rownames(express_rec)?<-express_rec[,1]
express_rec?<-?express_rec[(-1)]#表達矩陣的處理巷挥；
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))#確保表達矩陣的列名與分組矩陣行名相一致；
構(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進行差異表達分析雏节，返回 results可用的DESeqDataSet對象
>?dds?<-?DESeq(dds)#DESeq進行標準化；
estimating?size?factors
estimating?dispersions
gene-wise?dispersion?estimates
mean-dispersion?relationship
final?dispersion?estimates
fitting?modelandtesting
--?replacing?outliersandrefittingfor2819genes
--?DESeq?argument'minReplicatesForReplace'=7
--?original?counts?are?preserved?in?counts(dds)
estimating?dispersions
fitting?modelandtesting
可以用summary函數(shù)查看表達基因上下調(diào)分布基本情況
>?summary(res)#查看經(jīng)過標準化矩陣的基本情況高职；
outof24823withnonzero?totalreadcount
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'argumentof?results
[2]?see'independentFiltering'argumentof?results
最后提取差異表達基因钩乍，存入csv文件中。
edgeR進行差異分析
這個包的操作步驟與limma包類似怔锌，首先就是讀入數(shù)據(jù)寥粹，創(chuàng)建表達变过、分組矩陣。
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)#讀取表達矩陣;
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)建表達矩陣牵啦；
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對象進行標準化(TMM算法)?
#創(chuàng)建設(shè)計矩陣哈雏，跟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)出來
接下來依次進行構(gòu)建DGEList對象、利用TMM算法對數(shù)據(jù)進行標準化罪针、估計離散值彭羹、數(shù)據(jù)擬合、創(chuàng)建對比矩陣泪酱、對數(shù)據(jù)做QL-text檢驗派殷、差異表達基因?qū)懭隿sv文件中
dge?<-?DGEList(counts?=?express_rec,group=group)#構(gòu)建DEList對象；
y?<-?calcNormFactors(dge)#利用calcNormFactor函數(shù)對DEList對象進行標準化(TMM算法)?
y?<-?estimateDisp(y,design)#估計離散值（Dispersion）
fit?<-?glmQLFit(y,?design,?robust=TRUE)#進一步通過quasi-likelihood?(QL)擬合NB模型
TU.vs.NO?<-?makeContrasts(Tumor-Normal,?levels=design)#這一步主要構(gòu)建比較矩陣墓阀；
res?<-?glmQLFTest(fit,?contrast=TU.vs.NO)#用QL?F-test進行檢驗
#?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')
以上就是做差異分析三種R包的使用方法斯撮，關(guān)于本篇文章涉及到的完整源碼的獲取方式：關(guān)注公眾號：程序員大飛经伙，后臺回復(fù)關(guān)鍵詞：TCGA差異分析即可。
