在使用Bioconductor 分析芯片數(shù)據(jù)（一）中我們使用的是芯片的原始數(shù)據(jù)（即*.cel文件）進(jìn)行分析废境，但在大部分情況下我們是可以利用GEO數(shù)據(jù)庫(kù)中已有的表達(dá)矩陣來(lái)進(jìn)行數(shù)據(jù)分析的忍疾。
本次我們選用GSE17975數(shù)據(jù)穷劈，這是一個(gè)雙通道芯片數(shù)據(jù)鳖藕，研究Kawasaki disease的致病源頭大磺。但是奇怪的是作者沒(méi)有上傳芯片源文件（*.cel文件）憔四，因此我們直接使用 表達(dá)矩陣數(shù)據(jù)進(jìn)行分析赋咽。除了芯片外篙螟，我們還選擇了一個(gè)Kawasaki disease研究相關(guān)轉(zhuǎn)錄組數(shù)據(jù)，利用這兩種方法得到的差異基因的交集進(jìn)行后續(xù)分析芯肤。
1巷折、下載數(shù)據(jù)
這里依舊使用GEOquery包來(lái)下載數(shù)據(jù)：
> setwd('D:/TCGA/microarray_analysis/GSE17975/data/')
# 直接下載表達(dá)矩陣
> gmat <- getGEO('GSE17975',destdir = '.')
# 查看下里面是什么
> gmat
$GSE17975_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 22839 features, 3 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM450153 GSM450154 GSM450155
varLabels: title geo_accession ... disease state:ch2 (48 total)
varMetadata: labelDescription
featureData
featureNames: AGhsA010102 AGhsA010103 ... AGhsS010824 (22839 total)
fvarLabels: ID Description ... SPOT_ID (13 total)
fvarMetadata: Column Description labelDescription
experimentData: use 'experimentData(object)'
Annotation: GPL1291
> class(gmat)
[1] "list"
> class(gmat[[1]])
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
# 得到表達(dá)矩陣
> exprset <- exprs(gmat[[1]]); dim(exprset)
[1] 22839 3 # 可以看到表達(dá)矩陣中共有22839個(gè)基因（features），3個(gè)樣本
# 使用pData還可以查看更詳細(xì)的實(shí)驗(yàn)相關(guān)信息
> pdata <- pData(exprset)纷妆； dim(pdata)
[1] 3 48 # 有3行48列實(shí)驗(yàn)數(shù)據(jù)信息
### 最后查看一下表達(dá)矩陣信息
> head(exprset)
GSM450153 GSM450154 GSM450155
AGhsA010102 0.8335 0.0257 0.0841
AGhsA010103 0.0649 -0.3715 -0.2740
AGhsA010105 -0.3237 -0.7539 -0.2793
AGhsA010106 0.4844 -0.6897 -0.4521
AGhsA010108 -1.7418 -0.8783 -3.5804
AGhsA010109 -0.3401 -1.2311 -0.5332
> class(exprset)
[1] "matrix"
# 轉(zhuǎn)換成data.frame數(shù)據(jù)結(jié)構(gòu)方便后面使用
> exprset <- as.data.frame(exprset)
從上面可以發(fā)現(xiàn)盔几，第一晴弃，雖然我們是直接下載的表達(dá)矩陣掩幢，但是其中也含有許多的實(shí)驗(yàn)數(shù)據(jù)信息以及芯片平臺(tái)信息；第二通過(guò)直接查看表達(dá)矩陣他已經(jīng)經(jīng)過(guò)了標(biāo)準(zhǔn)化上鞠。
GEO Platform (GPL) 芯片平臺(tái)
GEO Sample (GSM) 樣本ID號(hào)
GEO Series (GSE) study的ID號(hào)
GEO Dataset (GDS) 數(shù)據(jù)集的ID號(hào)
2际邻、篩選差異基因
limma進(jìn)行差異表達(dá)分析需要準(zhǔn)備好3個(gè)數(shù)據(jù)：表達(dá)矩陣，分組矩陣芍阎，差異比較矩陣世曾。表達(dá)矩陣我們上面已經(jīng)從GEO下載得到，只需要將芯片ID注釋為基因ID即可（也可在差異分析后進(jìn)行注釋ID）谴咸，但這里的芯片平臺(tái)比較罕見(jiàn)轮听，找了半個(gè)小時(shí)沒(méi)找到對(duì)應(yīng)的注釋包，因此直接下載其對(duì)應(yīng)的注釋包寫代碼進(jìn)行注釋岭佳。
> gpl1291 <- getGEO(filename='GPL1291.soft')
> colnames(Table(gpl1291))
[1] "ID" "Description" "CDS_ID"
[4] "GB_ACC" "Entrez Gene ID" "Symbol"
[7] "GeneName" "IPI ID" "GO Molecular Function"
[10] "GO Biological Process" "GO Cellular Component" "ORF"
[13] "SPOT_ID"
# 從上面可以看到GPL1291.soft文件中第一列為芯片ID血巍，第5列和第6列為我們常用的ID，這里使用第6列Gene Symbol進(jìn)行注釋
> symbol <- Table(gpl1291)[c('ID','Symbol')]; head(symbol)
ID Symbol
1 AGhsA010101 IFNA8
2 AGhsA010102 OGG1
3 AGhsA010103 KIR2DL2
4 AGhsA010104 FGFR2
5 AGhsA010105 GAGE1
6 AGhsA010106 VCX
### 最后我們只需要將兩個(gè)數(shù)據(jù)框merge一下就可以了
> exprset1 <- merge(x=exprset,y=symbol,by='ID',all.x=T);head(exprset1)
ID GSM450153 GSM450154 GSM450155 Symbol
1 AGhsA010102 0.8335 0.0257 0.0841 OGG1
2 AGhsA010103 0.0649 -0.3715 -0.2740 KIR2DL2
3 AGhsA010105 -0.3237 -0.7539 -0.2793 GAGE1
4 AGhsA010106 0.4844 -0.6897 -0.4521 VCX
5 AGhsA010108 -1.7418 -0.8783 -3.5804 EEF1A1
6 AGhsA010109 -0.3401 -1.2311 -0.5332 DAZ1
### 查看一下注釋有沒(méi)有一對(duì)多或者未注釋上的情況
> dim(exprset);dim(exprset1)
[1] 22839 4
[1] 22839 5
> exprset1.na <- exprset1[is.na(exprset1$Symbol)]; dim(exprset1.na)
[1] 22839 0
### 最后刪除ID列珊随，并把symbol列調(diào)到第一列
> exprset2 <- exprset1[,-1]
> exprset3 <- exprset2[c(4,1,2,3)]; head(exprset3)
Symbol GSM450153 GSM450154 GSM450155
1 OGG1 0.8335 0.0257 0.0841
2 KIR2DL2 0.0649 -0.3715 -0.2740
3 GAGE1 -0.3237 -0.7539 -0.2793
4 VCX 0.4844 -0.6897 -0.4521
5 EEF1A1 -1.7418 -0.8783 -3.5804
6 DAZ1 -0.3401 -1.2311 -0.5332
從上面可以看到芯片ID注釋的結(jié)果很好述寡，全部注釋上了而且沒(méi)有一對(duì)多的情況柿隙。
但是真的沒(méi)有問(wèn)題嗎？會(huì)不會(huì)不是NA而是別的錯(cuò)誤沒(méi)有被我們發(fā)現(xiàn)鲫凶？讓我們?cè)賮?lái)做一次檢查禀崖。
### 查看一下有沒(méi)有空白字符
> failed <- exprset3[exprset3$Symbol=='',]
> head(failed)
Symbol GSM450153 GSM450154 GSM450155
16 -0.6712 -0.9269 -1.8417
50 -1.2345 0.1686 -0.5564
56 -0.5735 -1.0439 -1.5821
100 2.1925 -0.2059 0.4510
102 -2.1779 -0.6038 -0.3111
120 -0.4600 0.3663 -0.5886
> dim(failed)
[1] 7025 4 # (?Д?≡?Д?) 驚呆了有沒(méi)有，竟然有這么多沒(méi)有注釋上
### 難道是merge出錯(cuò)了螟炫？不應(yīng)該呀波附，查看一下官方給的注釋文件
> official.failed <- symbol[symbol$Symbol=='',]
> dim(official.failed)
[1] 10985 2 # 再次驚呆，官方給的注釋文件竟然有這么多空白昼钻。叶雹。。
### 最后换吧，去除這些未注釋的ID
> exprset4 <- exprset3[exprset3$Symbol!='',]; dim(exprset4)
[1] 15814 4
最后我們得到了15814個(gè)注釋到的基因折晦，相比之前少了7025個(gè)空的。所以沾瓦，看起來(lái)簡(jiǎn)簡(jiǎn)單單的一步注釋满着，也到處都是坑啊。
由于數(shù)據(jù)比較老贯莺，這里各列的數(shù)據(jù)其實(shí)都是表示logFC风喇，原文作者也是直接根據(jù)logFC選擇的差異基因進(jìn)行后續(xù)分析的，因此這里我們也是直接使用3組平行試驗(yàn)均值的|logFC|>1篩選差異基因進(jìn)行后續(xù)分析缕探。
> exprset4$mean <- (exprset4$GSM450153+exprset4$GSM450154+exprset4$GSM450155)/3
# 設(shè)定閾值|mean| >1 來(lái)篩選差異基因
> microexprgenes.differ <- exprset4[abs(exprset4$mean)>1,]
> microexprgenes.differ <- microexprgenes.differ[c(1,2,3,4)]
> dim(microexprgenes.differ)
[1] 938 4
這里我們篩選得了938個(gè)差異基因魂莫。然后去一下重復(fù)的gene symbol：
> microexprgenes.differ.dup <- microexprgenes.differ[duplicated(microexprgenes.differ$Symbol),]
> dim(microexprgenes.differ.dup)
[1] 99 4
> head(microexprgenes.differ.dup)
Symbol GSM450153 GSM450154 GSM450155
28 EEF1A1 -1.5649 -1.5187 -4.8573
677 FCGR1A 2.7853 2.5568 2.1375
1304 MATR3 -0.6394 -1.8890 -1.5522
2662 CCNB1IP1 -1.3921 -1.7959 -1.0923
3349 DHX9 -0.6461 -1.2042 -3.1203
4523 PSMB3 1.4895 1.1023 1.1603
> microexprgenes.differ <- microexprgenes.differ[!duplicated(microexprgenes.differ$Symbol),]
> dim(microexprgenes.differ)
[1] 839 4
最后我們得到了839個(gè)差異基因。
3爹耗、轉(zhuǎn)錄組差異基因篩選
這里直接使用的數(shù)據(jù)同樣也是關(guān)于Kawasaki disease研究的耙考，GEO數(shù)據(jù)集為GSE64486，我們這里直接采用已經(jīng)處理好的表達(dá)矩陣潭兽，但是我們這里使用的表達(dá)矩陣數(shù)據(jù)為GSE64486_untreated_vs_control_foldgenes_broad.txt.gz倦始，即未治療的與對(duì)照組。
由于這里的數(shù)據(jù)是處理好的山卦，我們可以直接設(shè)置logFC和P-value閾值來(lái)篩選差異基因鞋邑。這里設(shè)置的閾值為P-value<0.05, |logFC|>1。
# 先讀取數(shù)據(jù)
> exprgenes <- read.table('GSE64486_untreated_vs_control_foldgenes_broad.txt',header = T, stringsAsFactors = F,sep='\t');head(exprgenes)
ID P.Value Q.Value adj.P.Val Gene.Symbol RefSeq.ID Entrez.ID Fold.Change
1 ENSG00000211896 7.601798e-25 1.645219e-20 1.645219e-20 IGHG1 NA 3500 101.77408
2 ENSG00000138755 3.384964e-09 1.495083e-06 1.495083e-06 CXCL9 NA 4283 45.22891
3 ENSG00000211897 3.600406e-16 8.657976e-13 8.657976e-13 IGHG3 NA 3502 45.06493
4 ENSG00000211899 9.557181e-17 2.711637e-13 2.711637e-13 IGHM NA 3507 41.89091
5 ENSG00000132465 1.444665e-14 2.233298e-11 2.233298e-11 IGJ NA 3512 36.47332
6 ENSG00000196735 3.890287e-21 2.806518e-17 2.806518e-17 HLA-DQA1 NA 3117 35.58685
# 篩選差異基因
> exprgenes.differ <- exprgenes[exprgenes$P.Value<0.05 & abs(exprgenes$Fold.Change)>1,]
> dim(exprgenes);dim(exprgenes.differ)
[1] 43285 8
[1] 7099 8
可以看到账蓉，這里一共篩選得到了7099個(gè)枚碗，得到的基因比較多，可能閾值太寬了铸本，暫且先這樣接著做吧肮雨。
最后，還是檢查一下作者提供的數(shù)據(jù)注釋的gene symbol有沒(méi)有空的或者重復(fù)的：
> failed.symbol <- exprgenes.differ[duplicated(exprgenes.differ$Gene.Symbol),]
> head(failed.symbol)
ID P.Value Q.Value adj.P.Val Gene.Symbol RefSeq.ID Entrez.ID Fold.Change
120 ENSG00000252493 1.048966e-09 5.470419e-07 5.470419e-07 NA NA 7.509095
124 ENSG00000235300 3.582412e-09 1.566310e-06 1.566310e-06 NA NA 7.389821
132 ENSG00000251320 4.767764e-13 5.291607e-10 5.291607e-10 NA NA 7.135815
160 ENSG00000252862 7.716183e-16 1.757868e-12 1.757868e-12 NA NA 6.585053
167 ENSG00000251301 1.863953e-05 1.750134e-03 1.750134e-03 NA NA 6.419032
202 ENSG00000262151 4.473742e-05 3.500927e-03 3.500927e-03
> dim(failed.symbol)
[1] 2281 8
### 去除那些沒(méi)有g(shù)ene symbol 信息的
> exprgenes.differ <- exprgenes.differ[!duplicated(exprgenes.differ$Gene.Symbol),]
> dim(exprgenes.differ)
[1] 4818 8
可以看到归敬，一共2281個(gè)沒(méi)有注釋到gene symbol酷含，去除后一共得到了4818個(gè)有g(shù)ene symbol信息的差異基因鄙早。
我們可以發(fā)現(xiàn)，無(wú)論是芯片數(shù)據(jù)還是轉(zhuǎn)錄組數(shù)據(jù)椅亚，ID注釋的空白或者重復(fù)都是一個(gè)容易忽略的問(wèn)題限番。有機(jī)會(huì)再進(jìn)行較為深入的研究一下。
4呀舔、共有差異基因選取
前面我們分別得到了關(guān)于川崎病研究的芯片和轉(zhuǎn)錄組的差異基因弥虐，現(xiàn)在我們將這兩種數(shù)據(jù)合并取交集，以期得到與川崎病相關(guān)更為密切的基因媚赖。這里操作也比較簡(jiǎn)單霜瘪，直接使用R的merge函數(shù)即可。
> intersectgenes <- merge(x=microexprgenes.differ, y=exprgenes.differ, all=F,by.x='Symbol',by.y = 'Gene.Symbol')
> dim(intersectgenes)
[1] 128 11
> head(intersectgenes)
Symbol GSM450153 GSM450154 GSM450155 ID P.Value Q.Value adj.P.Val RefSeq.ID Entrez.ID Fold.Change
1 AIM1 -0.7155 -0.8391 -2.3808 ENSG00000112297 4.518321e-02 0.287315306 0.287315306 NA 202 2.555354
2 AK5 -1.9269 -0.6967 -0.5628 ENSG00000154027 7.907520e-05 0.005415775 0.005415775 NA 26289 1.350051
3 ANXA3 4.2110 0.8687 -0.1016 ENSG00000138772 3.122002e-02 0.229588586 0.229588586 NA 306 -2.328736
4 ARHGAP15 -0.0725 -1.5821 -2.0469 ENSG00000075884 6.321948e-05 0.004553170 0.004553170 NA 55843 5.064183
5 ASGR2 1.5563 1.4054 1.2066 ENSG00000161944 1.474010e-02 0.146976591 0.146976591 NA 433 2.061535
6 ATM -1.4344 -0.4961 -1.9324 ENSG00000149311 1.403235e-02 0.142613364 0.142613364 NA 472 2.313447
可以看到惧磺，這里我們一共得到了128個(gè)共有基因颖对，我們將使用這128個(gè)共有基因進(jìn)行后續(xù)分析。講這些基因?qū)懭胛募４娣奖愫竺媸褂谩?/p>
> write.table(intersectgenes,'../intersectgenes_logFC_broad.txt',row.names = F,sep='\t')
> write.table(intersectgenes$Symbol,'../intersectgenes.txt',row.names = F,sep='\t')
做個(gè)韋恩圖可視化看看：
vn <- venn.diagram(list(tissu=exprgenes.differ$Gene.Symbol, peripheral.blood=microexprgenes.differ$Symbol),filename='../intersectgenes.png',fill=c('red','green'),alpha=c(0.5,0.5),lty=2,cat.fontface=2,col='black',cat.col=c('red','green'))
總結(jié)
在這一部分个唧，我們利用了兩份GEO數(shù)據(jù)庫(kù)中的內(nèi)容，其中一份是芯片數(shù)據(jù)设预，另一個(gè)是轉(zhuǎn)錄組數(shù)據(jù)徙歼，對(duì)這兩個(gè)實(shí)驗(yàn)數(shù)據(jù)線選定閾值篩選出差異基因，然后綜合兩個(gè)數(shù)據(jù)尋找共同的差異基因鳖枕。
這部分操作比較簡(jiǎn)單魄梯，都是直接用的原作者給出的表達(dá)矩陣，然后進(jìn)行merge取交集即可耕魄。但是在操作過(guò)程中仍然有兩點(diǎn)是值得注意的：第一是各種基因ID之間的相互轉(zhuǎn)換關(guān)系画恰，這里可以利用各個(gè)數(shù)據(jù)庫(kù)提供的ID對(duì)應(yīng)文件進(jìn)行轉(zhuǎn)換即可；第二點(diǎn)吸奴，在本次實(shí)踐中我們發(fā)現(xiàn)了大量的未成功注釋的、空白的或重復(fù)的gene symbol缠局，這可以利用R中的duplicated函數(shù)在每次進(jìn)行數(shù)據(jù)分析時(shí)進(jìn)行一次檢查则奥。
