Hello，大家好！此前，小編給大家介紹了：如何從GEO數(shù)據(jù)庫(kù)下載芯片數(shù)據(jù)及相關(guān)樣本處理信息等等（不知道的可以戳這里哦）。這些芯片數(shù)據(jù)下載下來(lái)干嘛呢？下載必然是為了挖數(shù)據(jù)啦！指不定什么有意思的東西我就發(fā)現(xiàn)了呢？指不定老天爺眼睛一閉讓我發(fā)了文章了呢？
言歸正傳，拿到數(shù)據(jù)，分析的第一步往往是進(jìn)行基因差異表達(dá)分析，所以，針對(duì)芯片數(shù)據(jù)，我們就來(lái)給大家介紹一款基因差異表達(dá)分析的常用方法——R包limma。
數(shù)據(jù)簡(jiǎn)介與設(shè)置
為了方便演示，這里選擇了人的早幼粒細(xì)胞白血病細(xì)胞系NB4細(xì)胞的六個(gè)樣本數(shù)據(jù)（GSE2600），分析的輸入文件是下載的表達(dá)矩陣文件，而分析之前需要確保正確安裝和加載limma，同時(shí)需要對(duì)工作路徑進(jìn)行設(shè)置。
library('limma')
workdir="F:/GEO/20180520"
setwd(workdir)
數(shù)據(jù)處理
1、表達(dá)矩陣
數(shù)據(jù)為六個(gè)樣本，讀取數(shù)據(jù)之后，大家可以利用head()簡(jiǎn)單查看數(shù)據(jù)的情況等。
>?expreSet=read.csv2("GSE2600expressionMatrix.csv",?header?=T,?row.names?=1,check.names?=F)
>head(exprSet,3)
GSM49939GSM49940GSM49941GSM49942GSM49943GSM49944
1007_s_at23.013.826.575.994.984.6
1053_at1449.91826.72242.81508.81523.02355.5
117_at109.271.5106.7128.884.179.6
針對(duì)表達(dá)矩陣，需要查看其整體分布情況，可以利用boxplot()繪制box分布圖，GEO下載的表達(dá)矩陣數(shù)據(jù)基本上都是標(biāo)準(zhǔn)化的數(shù)據(jù)，可以由箱線圖的分布特點(diǎn)看出這些樣本的數(shù)據(jù)基本分布一致（中位數(shù)、上四分位數(shù)、下四分位數(shù)等等），如下圖結(jié)果：
#獲取樣品數(shù)量，并設(shè)置圖片顏色
n.sample?=?ncol(exprSet)
cols?=?rainbow(n.sample)
#利用boxplot()繪圖
pdf(file=paste(workdir,"/","Probe_expressionDistribution.pdf",sep=""),?width=24,?height=18)
par(cex?=0.7)
if(n.sample>40)?par(cex?=0.5)
boxplot(exprSet,col?=?cols,?main?="expression",?las?=2)
dev.off()
2、分組矩陣
確認(rèn)表達(dá)矩陣之后，可以由下載保存的樣本處理信息進(jìn)行分組，例如此處的樣本處理分組：CONTROL/INFECTED，經(jīng)過(guò)整理，分組信息大致如下,并基于分組信息構(gòu)建分組矩陣（design):
>group
Treatment
GSM49939???CONTROL
GSM49940???CONTROL
GSM49941???CONTROL
GSM49942??INFECTED
GSM49943??INFECTED
GSM49944??INFECTED
>?design?=?model.matrix(~?Treatment?+0,group)
>?colnames(design)?=?levels(as.factor(c("CONTROL","INFECTED")))
>?design
CONTROL?INFECTED
GSM4993910
GSM4994010
GSM4994110
GSM4994201
GSM4994301
GSM4994401
attr(,"assign")
[1]11
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1]"contr.treatment"
3、差異比較矩陣
基于分組矩陣的信息構(gòu)建差異比較矩陣（cont.matrix)，由差異比較矩陣顯示結(jié)果可知，是進(jìn)行INFECTED 與CONTROL之間的差異分析。
>cont.matrix = makeContrasts(INFECTED-CONTROL, levels=design)
>cont.matrix
Contrasts
LevelsINFECTED-CONTROL
CONTROL-1
INFECTED1
差異表達(dá)分析
差異表達(dá)分析主要是基于lmFit()、eBayes()、topTable()完成分析過(guò)程，并提取了主要的結(jié)果(tT)。
>?fit?=?lmFit(exprSet,?design)
>?fit2?=?contrasts.fit(fit,?cont.matrix)
>fit2?=?eBayes(fit2,0.01)
>tT?=?topTable(fit2,?adjust="fdr",?sort.by="logFC",?resort.by="P",n=Inf)
>
tT?=?subset(tT,?select=c("adj.P.Val","P.Value","logFC"))
>head(tT,15)
adj.P.ValP.ValuelogFC
223020_at0.999642.196175e-05746.10000
1555758_a_at0.999646.467722e-05-540.53333
218676_s_at0.999641.352768e-04-280.86667
237249_at0.999642.669173e-04-93.53333
225100_at0.999642.836527e-04-124.96667
217825_s_at0.999642.903446e-04-143.73333
222099_s_at0.999643.425427e-04493.13333
212634_at0.999644.221452e-04-166.06667
211499_s_at0.999644.391776e-04-129.56667
221098_x_at0.999644.805746e-0495.16667
208974_x_at0.999645.060448e-04947.76667
209670_at0.999645.113338e-04374.20000
202088_at0.999645.262646e-04-594.40000
219394_at0.999645.307063e-04-117.56667
212221_x_at0.999645.393084e-04347.43333
以上，就完成了分析過(guò)程，之后可以根據(jù)結(jié)果設(shè)定合適的閾值篩選出差異表達(dá)基因，并基于結(jié)果進(jìn)行圖形化展示，或者進(jìn)行富集分析、蛋白質(zhì)互作網(wǎng)絡(luò)分析等等，如此可以豐富分析結(jié)果，感興趣的同學(xué)可以去探索一下哦！
更多生物信息課程：
1. 文章越來(lái)越難發(fā)？是你沒(méi)發(fā)現(xiàn)新思路，基因家族分析發(fā)2-4分文章簡(jiǎn)單快速，學(xué)習(xí)鏈接：基因家族分析實(shí)操課程、基因家族文獻(xiàn)思路解讀
2. 轉(zhuǎn)錄組數(shù)據(jù)理解不深入？圖表看不懂？點(diǎn)擊鏈接學(xué)習(xí)深入解讀數(shù)據(jù)結(jié)果文件，學(xué)習(xí)鏈接：轉(zhuǎn)錄組（有參）結(jié)果解讀；轉(zhuǎn)錄組（無(wú)參）結(jié)果解讀
3. 轉(zhuǎn)錄組數(shù)據(jù)深入挖掘技能-WGCNA，提升你的文章檔次，學(xué)習(xí)鏈接：WGCNA-加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析
4. 轉(zhuǎn)錄組數(shù)據(jù)怎么挖掘？學(xué)習(xí)鏈接：轉(zhuǎn)錄組標(biāo)準(zhǔn)分析后的數(shù)據(jù)挖掘、轉(zhuǎn)錄組文獻(xiàn)解讀
5.?微生物16S/ITS/18S分析原理及結(jié)果解讀、OTU網(wǎng)絡(luò)圖繪制、cytoscape與網(wǎng)絡(luò)圖繪制課程
6. 生物信息入門到精通必修基礎(chǔ)課，學(xué)習(xí)鏈接：linux系統(tǒng)使用、perl入門到精通、perl語(yǔ)言高級(jí)、R語(yǔ)言畫圖
7. 醫(yī)學(xué)相關(guān)數(shù)據(jù)挖掘課程，不用做實(shí)驗(yàn)也能發(fā)文章，學(xué)習(xí)鏈接：TCGA-差異基因分析、GEO芯片數(shù)據(jù)挖掘、GSEA富集分析課程、TCGA臨床數(shù)據(jù)生存分析、TCGA-轉(zhuǎn)錄因子分析、TCGA-ceRNA調(diào)控網(wǎng)絡(luò)分析
8.其他課程鏈接：二代測(cè)序轉(zhuǎn)錄組數(shù)據(jù)自主分析、NCBI數(shù)據(jù)上傳、二代測(cè)序數(shù)據(jù)解讀。
