以下是B站生信技能樹GEO數(shù)據(jù)庫(kù)挖掘的課程筆記
主要內(nèi)容及學(xué)習(xí)目的：
- 介紹GEO數(shù)據(jù)庫(kù)：了解數(shù)據(jù)存放位置；
- 介紹GSE項(xiàng)目的3種下載方式；
- 介紹ID轉(zhuǎn)換：使用R語(yǔ)言技巧實(shí)現(xiàn)基因ID之間的轉(zhuǎn)換，我們下載的數(shù)據(jù)通常使用的是不同的芯片探針，它們有不同的探針I(yè)D號(hào)我們需要把它轉(zhuǎn)化成ENTREZID或SYMBOLID才能被大眾認(rèn)知；
- 介紹表達(dá)矩陣的相關(guān)可視化及歸一化：從GEO數(shù)據(jù)庫(kù)下載的是作者處理好的矩陣，我們需要會(huì)判別它是否符合要求，并學(xué)會(huì)提取分組信息；
- 比較各組基因的表達(dá)量得到差異表達(dá)基因list或感興趣基因集；
- 得到差異表達(dá)基因list后做富集分析；
- 用GSEA軟件做一些圖；
通過(guò)閱讀文章提煉GEO數(shù)據(jù)庫(kù)挖掘的脈絡(luò)：選擇GSE ---> 得到表達(dá)矩陣 ---> control VS treatment 進(jìn)行差異分析 ---> 得到差異表達(dá)基因list ---> 5大數(shù)據(jù)庫(kù)的注釋 ---> PPI等網(wǎng)絡(luò)
接下來(lái)我們按照上面的分析思路，一步一步進(jìn)行講解
1.了解GEO數(shù)據(jù)庫(kù)，找到文章的GSE編號(hào)
任何一篇GEO數(shù)據(jù)挖掘文章，都可以找到它的GSE編號(hào)，找到后我們把網(wǎng)址最后的GSE編號(hào)修改一下，直接去網(wǎng)頁(yè)粘貼并轉(zhuǎn)到就能看到該編號(hào)在GEO數(shù)據(jù)庫(kù)的詳細(xì)頁(yè)面：
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872
我們看下GEO數(shù)據(jù)庫(kù)的主頁(yè)
我們只需要知道這三個(gè)概念就好：
- GEO Platform (GPL)
- GEO Sample (GSM)
-
GEO Series (GSE)
理解起來(lái)也很容易。一篇文章可以有一個(gè)或者多個(gè)GSE數(shù)據(jù)集，一個(gè)GSE數(shù)據(jù)集里面可以有一個(gè)或者多個(gè)GSM樣本。而每個(gè)數(shù)據(jù)集都有著自己對(duì)應(yīng)的芯片平臺(tái)，就是GPL。
點(diǎn)開后向下滑動(dòng)找到GPL的表格信息
把GPL表格信息向右滑動(dòng)，找到
gene_assignment那列，把//做為切割符，取出第二個(gè)字符就是真正的基因名，這時(shí)探針就和基因完美匹配啦~
知道如何找到任何一篇文章的數(shù)據(jù)存放位置，接下來(lái)就要下載數(shù)據(jù)進(jìn)行分析了。
2.下載數(shù)據(jù)
下載數(shù)據(jù)的3種方式：
一. 直接下載rawdata —— 不推薦使用
二. 從網(wǎng)頁(yè)上直接下載表達(dá)矩陣 ---> 讀入R里
表達(dá)矩陣下載到本地后要讀入到R里：
a = read.table(file="./GSE42872_series_matrix.txt.gz",
header = T,sep = "\t",quote = "",fill = T,
comment.char = "!")
在讀入下載好的表達(dá)矩陣時(shí)，為什么要加那么多參數(shù)才能下載成功？
我們首先需要在電腦上解壓并打開文本文件，根據(jù)文件的樣子選擇參數(shù)
我們看
a的rowname是行號(hào)，沒(méi)有意義的，需要轉(zhuǎn)成探針I(yè)D號(hào)即a的第一列：rownames(a)= a[,1]
> head(a)
X.ID_REF. X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
1 7892501 7.24559 6.80686 7.73301 6.18961 7.05335 7.20371
2 7892502 6.82711 6.70157 7.02471 6.20493 6.76554 6.24252
3 7892503 4.39977 4.50781 4.88250 4.36295 4.18137 4.73492
4 7892504 9.48025 9.67952 9.63074 9.69200 9.91324 9.65897
5 7892505 4.54734 4.45247 5.11753 4.87307 5.15505 3.99340
6 7892506 6.80701 6.90597 6.72472 6.77028 6.77058 6.77685
> rownames(a)
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13"
> rownames(a)= a[,1]
> head(a)
X.ID_REF. X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501 7892501 7.24559 6.80686 7.73301 6.18961 7.05335 7.20371
7892502 7892502 6.82711 6.70157 7.02471 6.20493 6.76554 6.24252
7892503 7892503 4.39977 4.50781 4.88250 4.36295 4.18137 4.73492
7892504 7892504 9.48025 9.67952 9.63074 9.69200 9.91324 9.65897
7892505 7892505 4.54734 4.45247 5.11753 4.87307 5.15505 3.99340
7892506 7892506 6.80701 6.90597 6.72472 6.77028 6.77058 6.77685
此時(shí)，a的列名就是探針I(yè)D號(hào)了。但是現(xiàn)在還是不符合預(yù)期，我們還要把RefSeq ID那一列去掉，也就是去掉此時(shí)的第一列：a = a[,-1]
> a = a[,-1]
> head(a)
X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501 7.24559 6.80686 7.73301 6.18961 7.05335 7.20371
7892502 6.82711 6.70157 7.02471 6.20493 6.76554 6.24252
7892503 4.39977 4.50781 4.88250 4.36295 4.18137 4.73492
7892504 9.48025 9.67952 9.63074 9.69200 9.91324 9.65897
7892505 4.54734 4.45247 5.11753 4.87307 5.15505 3.99340
7892506 6.80701 6.90597 6.72472 6.77028 6.77058 6.77685
這就是由樣本和探針組成的表達(dá)矩陣
三. 在R里使用GSE號(hào)和GEOquery包從GEO數(shù)據(jù)庫(kù)上直接下載——最推薦使用下載方式
library(GEOquery)
eSet <- getGEO("GSE42872",
destdir = '.', #下載在當(dāng)前目錄
getGPL = F) #平臺(tái)信息不要
使用以上代碼就可以將GSE42872數(shù)據(jù)下載到R里當(dāng)前工作目錄并賦值給eSet，下載完成后要注意檢查數(shù)據(jù)文件的完整性——看我們下載的數(shù)據(jù)大小是否大于等于官網(wǎng)上給的大小。如果我們下載的數(shù)據(jù)內(nèi)存大于官網(wǎng)上的那沒(méi)事兒，如果小于官網(wǎng)上的那下載的數(shù)據(jù)就不完整。
2.1得到表達(dá)矩陣
我們用方法2下載的a和用方法3下載的eSet都是GSE42872數(shù)據(jù)，但它們是不一樣的：
> class(a)
[1] "data.frame"
> class(eSet)
[1] "list"
我們可以看到a是一個(gè)數(shù)據(jù)框，eSet是一個(gè)列表這里我們稱它為對(duì)象。
得到eSet對(duì)象里包含著各種各樣的信息：表達(dá)矩陣、芯片如何設(shè)計(jì)的、樣本如何分組 等等~
eSet是一個(gè)大列表，我們需要從中提取出表達(dá)矩陣，才能進(jìn)行后續(xù)的操作。
為什么？因?yàn)橐粋€(gè)GSE號(hào)里面對(duì)應(yīng)多種芯片平臺(tái)數(shù)據(jù)，我們使用GSE號(hào)下載數(shù)據(jù)就會(huì)把所有芯片平臺(tái)的數(shù)據(jù)整合到一個(gè)list里面，每個(gè)list里的元素存放一個(gè)平臺(tái)的表達(dá)矩陣。我們的數(shù)據(jù)只有一個(gè)平臺(tái)所以eSet列表里就只有一個(gè)元素：
使用列表取子集的方法提取
eSet列表里的第一個(gè)元素：eSet[[1]]；并使用exprs函數(shù)把它轉(zhuǎn)化成矩陣：exp <- exprs(eSet[[1]])
> eSet[[1]]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
varLabels: title geo_accession ... cell type:ch1 (34 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 24469106
Annotation: GPL6244
> exp <- exprs(eSet[[1]])
> exp[1:4,1:4]
GSM1052615 GSM1052616 GSM1052617 GSM1052618
7892501 7.24559 6.80686 7.73301 6.18961
7892502 6.82711 6.70157 7.02471 6.20493
7892503 4.39977 4.50781 4.88250 4.36295
7892504 9.48025 9.67952 9.63074 9.69200
這時(shí)我們?cè)倏从煞椒?得到的表達(dá)矩陣a，和由方法3得到的表達(dá)矩陣exp是一模一樣的：
> head(a)
X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501 7.24559 6.80686 7.73301 6.18961 7.05335 7.20371
7892502 6.82711 6.70157 7.02471 6.20493 6.76554 6.24252
7892503 4.39977 4.50781 4.88250 4.36295 4.18137 4.73492
7892504 9.48025 9.67952 9.63074 9.69200 9.91324 9.65897
7892505 4.54734 4.45247 5.11753 4.87307 5.15505 3.99340
7892506 6.80701 6.90597 6.72472 6.77028 6.77058 6.77685
> head(exp)
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
7892501 7.24559 6.80686 7.73301 6.18961 7.05335 7.20371
7892502 6.82711 6.70157 7.02471 6.20493 6.76554 6.24252
7892503 4.39977 4.50781 4.88250 4.36295 4.18137 4.73492
7892504 9.48025 9.67952 9.63074 9.69200 9.91324 9.65897
7892505 4.54734 4.45247 5.11753 4.87307 5.15505 3.99340
7892506 6.80701 6.90597 6.72472 6.77028 6.77058 6.77685
基因ID之間的轉(zhuǎn)換
entrez ID<---probe_id--->symbol ID
分兩步走：
- 過(guò)濾
probe_id，得到每個(gè)基因所對(duì)應(yīng)的唯一的probe_id- 得到
probe_id與symbol ID這件的轉(zhuǎn)換關(guān)系
使用R語(yǔ)言技巧實(shí)現(xiàn)基因ID之間的轉(zhuǎn)換，我們下載的數(shù)據(jù)通常使用的是不同的芯片探針，它們有不同的探針I(yè)D（probe_id）我們需要把它轉(zhuǎn)化成entrez ID或symbol ID才能被大眾認(rèn)知；
所以接下來(lái)，我們學(xué)習(xí)怎樣將探針I(yè)D（probe_id）轉(zhuǎn)換成entrez ID和symbol ID
ID轉(zhuǎn)換的第一步必須要加載特定的R包，下載哪個(gè)包，需要根據(jù)GPL來(lái)定
> eSet[[1]]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
varLabels: title geo_accession ... cell type:ch1 (34 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 24469106
Annotation: GPL6244
我們看到GPL號(hào)是：GPL6244
> eSet[[1]]@annotation
[1] "GPL6244"
用R獲取芯片探針與基因的對(duì)應(yīng)關(guān)系三部曲-bioconductor里搜索GPL6244所對(duì)應(yīng)的R包，發(fā)現(xiàn)是hugene10sttranscriptcluster
所以我們加載此包：
library(hugene10sttranscriptcluster.db) 注意加上后綴 .db
> library(hugene10sttranscriptcluster.db)
> ls("package:hugene10sttranscriptcluster.db")
[1] "hugene10sttranscriptcluster" "hugene10sttranscriptcluster.db"
[3] "hugene10sttranscriptcluster_dbconn" "hugene10sttranscriptcluster_dbfile"
[5] "hugene10sttranscriptcluster_dbInfo" "hugene10sttranscriptcluster_dbschema"
[7] "hugene10sttranscriptclusterACCNUM" "hugene10sttranscriptclusterALIAS2PROBE"
[9] "hugene10sttranscriptclusterCHR" "hugene10sttranscriptclusterCHRLENGTHS"
[11] "hugene10sttranscriptclusterCHRLOC" "hugene10sttranscriptclusterCHRLOCEND"
[13] "hugene10sttranscriptclusterENSEMBL" "hugene10sttranscriptclusterENSEMBL2PROBE"
[15] "hugene10sttranscriptclusterENTREZID" "hugene10sttranscriptclusterENZYME"
[17] "hugene10sttranscriptclusterENZYME2PROBE" "hugene10sttranscriptclusterGENENAME"
[19] "hugene10sttranscriptclusterGO" "hugene10sttranscriptclusterGO2ALLPROBES"
[21] "hugene10sttranscriptclusterGO2PROBE" "hugene10sttranscriptclusterMAP"
[23] "hugene10sttranscriptclusterMAPCOUNTS" "hugene10sttranscriptclusterOMIM"
[25] "hugene10sttranscriptclusterORGANISM" "hugene10sttranscriptclusterORGPKG"
[27] "hugene10sttranscriptclusterPATH" "hugene10sttranscriptclusterPATH2PROBE"
[29] "hugene10sttranscriptclusterPFAM" "hugene10sttranscriptclusterPMID"
[31] "hugene10sttranscriptclusterPMID2PROBE" "hugene10sttranscriptclusterPROSITE"
[33] "hugene10sttranscriptclusterREFSEQ" "hugene10sttranscriptclusterSYMBOL"
[35] "hugene10sttranscriptclusterUNIGENE" "hugene10sttranscriptclusterUNIPROT"
通過(guò)命令：
ls("package:hugene10sttranscriptcluster.db")
我們可以看到這個(gè)包里面有很多數(shù)據(jù)集，想要得到probe_id和symbol的對(duì)應(yīng)關(guān)系要用hugene10sttranscriptclusterSYMBOL數(shù)據(jù)集，用toTable函數(shù)提取數(shù)據(jù)集里面的信息：
> ids=toTable(hugene10sttranscriptclusterSYMBOL)
> head(ids)
probe_id symbol
1 7896759 LINC01128
2 7896761 SAMD11
3 7896779 KLHL17
4 7896798 PLEKHN1
5 7896817 ISG15
6 7896822 AGRN
現(xiàn)在我們查看下一共多少個(gè)基因？一萬(wàn)八千多個(gè)基因
> length(unique(ids$symbol))
[1] 18834
unique函數(shù)是用來(lái)：Extract Unique Elements 去除重復(fù)的symbol只提取不同的元素；length函數(shù)統(tǒng)計(jì)去重之后還有多少個(gè)基因。
再查看每個(gè)基因?qū)?yīng)多少個(gè)探針：
> tail(sort(table(ids$symbol)))
RPL41 UBTFL1 CDK11B UBE2D3 IGKC LRRFIP1
6 6 8 8 10 10
可以看到有的基因設(shè)計(jì)了10個(gè)探針或8個(gè)探針....
table() 函數(shù)可以生成頻數(shù)統(tǒng)計(jì)表，這里就是統(tǒng)計(jì)每個(gè)基因symbol出現(xiàn)的次數(shù)然后將其表格化；sort()函數(shù)將symbol出現(xiàn)的頻率從小到大進(jìn)行排序；tail()取最后6個(gè)即出現(xiàn)頻率最大的6個(gè)。
> table(sort(table(ids$symbol)))
1 2 3 4 5 6 8 10
18072 599 132 16 5 6 2 2
table一下我們可以看到，18072個(gè)基因設(shè)計(jì)了1個(gè)探針；599個(gè)基因設(shè)計(jì)了2個(gè)探針；132個(gè)基因設(shè)計(jì)了3個(gè)探針.....也就是說(shuō)大部分的基因只設(shè)計(jì)了1個(gè)探針。
其實(shí)一般基因都會(huì)設(shè)計(jì)很多探針的，我們下載的表達(dá)矩陣是作者處理之后的，把許多不好的探針都過(guò)濾掉了，我們處理作者的數(shù)據(jù)要默認(rèn)人家做的是對(duì)的，否則就要下載原始數(shù)據(jù)自己處理。
> table(rownames(exp) %in% ids$probe_id)
FALSE TRUE
13470 19827
發(fā)現(xiàn)有13470個(gè)探針沒(méi)有對(duì)應(yīng)的基因名；19827個(gè)探針有對(duì)應(yīng)的基因名。
x %in% y表示 x 的元素在y中嗎？然后返回邏輯值。rownames(exp)即表達(dá)矩陣exp的行名是文章數(shù)據(jù)中用到的所有探針I(yè)D（probe_id）；ids$probe_id是具有對(duì)應(yīng)基因的所有探針。所以返回的TRUE就是文章數(shù)據(jù)中有對(duì)應(yīng)基因的探針數(shù)。
現(xiàn)在我們對(duì)探針進(jìn)行過(guò)濾，把沒(méi)有對(duì)應(yīng)基因名的探針過(guò)濾掉：
> exp = exp[rownames(exp) %in% ids$probe_id,]
過(guò)濾的本質(zhì)就是矩陣取子集，如：matrix[2,]意思就是取矩陣matrix的第2行和所有的列。同理，我們這里exp[rownames(exp) %in% ids$probe_id,]就是取具有對(duì)應(yīng)基因的所有探針（行），和所有的列。
對(duì)比一下過(guò)濾之前和過(guò)濾之后的探針數(shù)量：
> table(rownames(exp) %in% ids$probe_id)
FALSE TRUE
13470 19827
> dim(exp)
[1] 33297 6
> exp = exp[rownames(exp) %in% ids$probe_id,]
> dim(exp)
[1] 19827 6
可以發(fā)現(xiàn)過(guò)濾之前有33297個(gè)探針，過(guò)濾之后就剩下19827個(gè)探針了。
然后，我們使用match函數(shù)把ids里的探針順序改一下，使ids里探針順序和我們表達(dá)矩陣的順序完全一樣：
ids=ids[match(rownames(exp),ids$probe_id),]
match()函數(shù)返回的是一個(gè)位置向量，該向量記錄著第一個(gè)參數(shù)中每個(gè)元素在第二個(gè)參數(shù)中的位置。所以，此時(shí)ids里的探針順序與表達(dá)矩陣exp的探針順序一一對(duì)應(yīng)：
> head(ids)
probe_id symbol
1 7896759 LINC01128
2 7896761 SAMD11
3 7896779 KLHL17
4 7896798 PLEKHN1
5 7896817 ISG15
6 7896822 AGRN
> head(exp)
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
7896759 8.75126 8.61650 8.81149 8.32067 8.41445 8.45208
7896761 8.39069 8.52617 8.43338 9.17284 9.10216 9.14120
7896779 8.20228 8.30886 8.18518 8.13322 8.06453 8.15884
7896798 8.41004 8.37679 8.27521 8.34524 8.35557 8.44409
7896817 7.72204 7.74572 7.78022 7.72308 7.53797 7.73401
7896822 9.19237 9.10929 9.03668 9.94821 9.96994 9.99839
既然已經(jīng)完全對(duì)應(yīng)上，我們就可以通過(guò)probe_id將表達(dá)矩陣exp進(jìn)行分組，將同一個(gè)symbol所對(duì)應(yīng)的多個(gè)探針?lè)殖刹煌慕M，并對(duì)每組探針進(jìn)行統(tǒng)計(jì)：計(jì)算每組中每行探針表達(dá)量的平均值（也就是每個(gè)探針在6個(gè)樣本中表達(dá)量的均值rowMeans(x)），再取平均值最大的那個(gè)探針作為該symbol所對(duì)應(yīng)的唯一探針，該組中的其它探針過(guò)濾掉，這樣每個(gè)symbol就對(duì)應(yīng)一個(gè)探針了，看下代碼是如何操作的：
tmp = by(exp,
ids$symbol,
function(x) rownames(x)[which.max(rowMeans(x))])
probes = as.character(tmp)
dim(exp)
exp = exp[rownames(exp) %in% probes,] # 過(guò)濾有多個(gè)探針的基因
dim(exp)
是不是沒(méi)有理解上面代碼在干些什么？沒(méi)關(guān)系，我們?cè)敿?xì)解釋一下：
by()函數(shù)在這里發(fā)揮的功能就是將表達(dá)矩陣exp中的探針?lè)纸M，同一個(gè)symbol所對(duì)應(yīng)的多個(gè)探針?lè)值揭唤M，并對(duì)每組探針進(jìn)行統(tǒng)計(jì)得到symbol所對(duì)應(yīng)的唯一探針。所以tmp里放著by()函數(shù)的統(tǒng)計(jì)結(jié)果即每個(gè)symbol所對(duì)應(yīng)的唯一探針I(yè)Dprobe_id，用probes = as.character(tmp)將結(jié)果變身為純字符型向量：
> head(tmp)
INDICES
A1CF A2M A2ML1 A3GALT2 A4GALT A4GNT
"7933640" "7960947" "7953775" "7914643" "8076497" "8090955"
> head(probes)
[1] "7933640" "7960947" "7953775" "7914643" "8076497" "8090955"
>
學(xué)習(xí)by()函數(shù)如何完成以上操作的。《R語(yǔ)言實(shí)戰(zhàn)》這本書上是這樣描述的
使用
by()分組計(jì)算描述性統(tǒng)計(jì)量，它可以一次返回若干個(gè)統(tǒng)計(jì)量。格式為：
by(data, INDICES, FUN)
其中data是一個(gè)數(shù)據(jù)框或矩陣；INDICES是一個(gè)因子或因子組成的列表，定義了分組；FUN是任意函數(shù)。
簡(jiǎn)單一句話理解就是：by()函數(shù)就是根據(jù)因子將整個(gè)data分成幾個(gè)小的data.frame，然后進(jìn)行運(yùn)算處理。
同理，我們這里：
by(exp, ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))])
第二個(gè)參數(shù)ids$symbol定義了分組，將第一參數(shù)—exp表達(dá)矩陣分成了若干個(gè)小矩陣，每個(gè)小矩陣?yán)锎娣胖粋€(gè)symbol所對(duì)應(yīng)的所有探針。第三個(gè)參數(shù)是我們自己定義的函數(shù)：計(jì)算每個(gè)小矩陣中每行探針表達(dá)量的平均值（也就是每個(gè)探針在6個(gè)樣本中表達(dá)量的均值rowMeans(x)），再取平均值最大的那個(gè)探針作為該symbol所對(duì)應(yīng)的唯一探針which.max(rowMeans(x))。
by()函數(shù)就可以返回每個(gè)分組里的統(tǒng)計(jì)結(jié)果，即每個(gè)symbol所對(duì)應(yīng)的唯一探針I(yè)Dprobe_id。
這時(shí)，探針I(yè)D和基因symbol就一一對(duì)應(yīng)了，將表達(dá)矩陣探針I(yè)D即exp表達(dá)矩陣的行名（rownames(exp)）換為基因symbol：
rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]
> head(exp)
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
LINC01128 8.75126 8.61650 8.81149 8.32067 8.41445 8.45208
SAMD11 8.39069 8.52617 8.43338 9.17284 9.10216 9.14120
KLHL17 8.20228 8.30886 8.18518 8.13322 8.06453 8.15884
PLEKHN1 8.41004 8.37679 8.27521 8.34524 8.35557 8.44409
ISG15 7.72204 7.74572 7.78022 7.72308 7.53797 7.73401
AGRN 9.19237 9.10929 9.03668 9.94821 9.96994 9.99839
此時(shí)，我們已經(jīng)將探針I(yè)D轉(zhuǎn)化成基因symbol了。
由上面的介紹我們可以看到，在轉(zhuǎn)換ID中最重要的是根據(jù)GPL平臺(tái)號(hào)找到所對(duì)應(yīng)的R注釋包，可是如果找不到GPL平臺(tái)對(duì)應(yīng)的R注釋包怎么辦呢？
答：我們不用GEO號(hào)進(jìn)行下載，而是下載平臺(tái)信息(GPL)，從平臺(tái)信息中選擇我們想要的列：探針名、基因名....
GPL里面的信息量特別大，下載特別考驗(yàn)網(wǎng)速。
gpl <- getGEO('GPL6480', destdir = ".")
colnames(Table(gpl))
head(Table(gpl)[,c(1,6,7)]) #看gpl對(duì)象中哪一列是我們想要的取出來(lái)，發(fā)現(xiàn)1/6/7列是我們想要的
write.csv(Table(gpl)[,c(1,6,7)],"GPL6480.csv") #把我們想要的部分即探針名對(duì)應(yīng)的基因名....存起來(lái)
獲取分組信息—group_list
分組信息就是告訴我們哪些組是control；哪些組是tumor。
使用pData函數(shù)獲取分組信息—group_list：
pd <- pData(eSet[[1]])
pData()函數(shù)可以得到每個(gè)樣本的描述信息，一般來(lái)說(shuō)數(shù)據(jù)框的第一列(title列)描述了哪些是control；哪些是treatment。
根據(jù)第一列所描述的信息我們自己創(chuàng)建分組信息
group_list：方法一：使用
stringr函數(shù)
library(stringr)
group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
group_list
stringr包用于字符串的處理，str_detect是該包里的函數(shù)，用來(lái)確定一個(gè)字符向量能否匹配一種模式。它返回一個(gè)與輸入向量具有同樣長(zhǎng)度的邏輯向量：
> str_detect(pd$title,"Control")
[1] TRUE TRUE TRUE FALSE FALSE FALSE
這里的輸入向量是數(shù)據(jù)框pd的第一列pd$title內(nèi)容，即由6個(gè)元素組成的字符型向量。str_detect()函數(shù)會(huì)自動(dòng)判斷Control，是否存在于pd$title向量的每一個(gè)元素中，存在返回TRUE，否則返回FALSE。
str_detect函數(shù)處理后我們?cè)偈褂?ifelse生成符合要求的分組信息group_list
> group_list
[1] "contorl" "contorl" "contorl" "treat" "treat" "treat"
方法二：自己造一個(gè)
我們已經(jīng)知道了前三個(gè)是control后三個(gè)是treatment，那就自己生成一個(gè)符合要求的分組信息：
group_list=c(rep("control",times=3),rep("treat",times=3))
group_list
3. 檢查表達(dá)矩陣
得到表達(dá)矩陣就是描述了某個(gè)基因在某個(gè)樣本的表達(dá)量。有了這個(gè)表達(dá)矩陣我們可以做后面的分析，第一步就是確定我們得到的表達(dá)矩陣是否正確：
- 查看管家基因的表達(dá)量
- 檢測(cè)分組之間是否有差異：PCA圖、熱圖和hclust圖等等
3.1 檢驗(yàn)常見基因的表達(dá)量
查看典型管家基因（如：GAPDH、ACTB）的表達(dá)量，如果表達(dá)量高于正常值，說(shuō)明我們數(shù)據(jù)沒(méi)問(wèn)題。
此時(shí)表達(dá)矩陣exp的行名已經(jīng)由探針I(yè)D轉(zhuǎn)換成基因名了，所以我們使用exp['GAPDH',]來(lái)提取該基因在所有樣品中的表達(dá)量。
> exp['GAPDH',]
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
14.3187 14.3622 14.3638 14.4085 14.3569 14.3229
> exp['ACTB',]
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
13.8811 13.9002 13.8822 13.7773 13.6732 13.5363
我們可以看到我們數(shù)據(jù)中兩個(gè)管家基因的表達(dá)量都偏高，符合預(yù)期。為什么知道它偏高呢？畫一個(gè)整體樣本所有基因的表達(dá)量的boxplot：boxplot(exp)
發(fā)現(xiàn)大部分基因的表達(dá)量都在8-9，而
GAPDH、ACTB在13-14左右，所以是偏高的。假如，我們發(fā)現(xiàn)管家基因表達(dá)量特別低，那我們就要思考是不是在提取表達(dá)矩陣的時(shí)候哪里出了問(wèn)題：比如ID轉(zhuǎn)換的時(shí)候轉(zhuǎn)換錯(cuò)了等等....
3.2 看表達(dá)矩陣的分布圖—畫圖看各個(gè)樣本的表達(dá)量
使用ggplot2畫各個(gè)樣本表達(dá)量的boxplot圖
# 準(zhǔn)備畫圖所需數(shù)據(jù)exp_L
library(reshape2)
head(exp)
exp_L = melt(exp)
head(exp_L)
colnames(exp_L)=c('symbol','sample','value')
head(exp_L)
# 獲得分組信息
library(stringr)
group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
group_list
exp_L$group = rep(group_list,each=nrow(exp))
head(exp_L)
# ggplot2畫圖
library(ggplot2)
p = ggplot(exp_L,
aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
##boxplot圖精修版
p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)
作圖函數(shù)看起來(lái)復(fù)雜我們拆開：
準(zhǔn)備作圖所需要的數(shù)據(jù)
exp_L---> 獲得分組信息并加到exp_L中 --->ggplot2作圖
我們先理解一下 exp_L 數(shù)據(jù)
> head(exp_L)
symbol sample value group
1 LINC01128 GSM1052615 8.75126 contorl
2 SAMD11 GSM1052615 8.39069 contorl
3 KLHL17 GSM1052615 8.20228 contorl
4 PLEKHN1 GSM1052615 8.41004 contorl
5 ISG15 GSM1052615 7.72204 contorl
6 AGRN GSM1052615 9.19237 contorl
> table(exp_L[,2])
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
18834 18834 18834 18834 18834 18834
> dim(exp_L)
[1] 113004 4
> 18834*6
[1] 113004
由以上代碼我們可以看到exp_L矩陣是這樣分布的：每個(gè)基因（18834個(gè)）在第一個(gè)樣本GSM1052615中的value值，每個(gè)基因（18834個(gè)）在第二個(gè)樣本中的value值....以此類推一共有6個(gè)樣本。
難點(diǎn)攻克：如何得到這樣的
exp_L矩陣呢？？？使用reshape2包
reshape2包是一套重構(gòu)和整合數(shù)據(jù)集的絕妙的萬(wàn)能工具。大致用法就是，需要首先將數(shù)據(jù)融合（melt），以使每一行都是唯一的標(biāo)識(shí)符-變量組合。然后將數(shù)據(jù)重塑（cast）為你想要的任何形狀。在重鑄過(guò)程中，你可以使用任何函數(shù)對(duì)數(shù)據(jù)進(jìn)行整合。
我們這里只用到這個(gè)包里的數(shù)據(jù)融合（melt）功能。
數(shù)據(jù)集的融合（melt）是將它重構(gòu)為這樣一種格式：每個(gè)測(cè)量變量（每個(gè)基因在每個(gè)樣本中的表達(dá)量）獨(dú)占一行，行中帶有要唯一確定這個(gè)測(cè)量所需的標(biāo)識(shí)符變量（基因symbol和樣本sample）。注意，必須指定要唯一確定每個(gè)測(cè)量所需的變量（也就是說(shuō)基因symbol和樣本sample必須對(duì)應(yīng)唯一的表達(dá)量），而表示測(cè)量變量名的變量將由程序?yàn)槟阕詣?dòng)創(chuàng)建（即表達(dá)量獨(dú)占一行后程序會(huì)自動(dòng)創(chuàng)建表達(dá)量所對(duì)應(yīng)的symbol和sample）。
說(shuō)成人話就是，以前exp矩陣是一個(gè)基因在6個(gè)樣本中的表達(dá)量占一行，melt后就會(huì)將表達(dá)量獨(dú)占一行。一個(gè)表達(dá)量的值需要有兩個(gè)定語(yǔ)才能唯一指定，即這個(gè)表達(dá)量是哪個(gè)樣本（sample）中的哪個(gè)基因（symbol）的。
從圖中可以看到兩個(gè)分組
control和treat基本在一條線上，這樣的數(shù)據(jù)說(shuō)明可以進(jìn)行后續(xù)比較，如果不在一條線上說(shuō)明有批次效應(yīng)batch infect，需要用limma包內(nèi)置函數(shù)normalizeBetweenArrays人工校正一下(Normalization)：
library(limma)
exp = normalizeBetweenArrays(exp)
關(guān)于畫樣本表達(dá)量的分布圖，除了上面介紹的boxplot，ggplot2還可以畫別的，看情況使用就好，不同的圖有不同的展現(xiàn)方式但都在展現(xiàn)同一個(gè)問(wèn)題那就是各個(gè)樣本的表達(dá)量，看自己喜歡用就好：
p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exp_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exp_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exp_L,aes(value,col=group))+geom_density()
print(p)
3.3 檢查樣本分組信息
檢查樣本分組信息，一般看PCA圖，hclust圖
hclust
# 更改表達(dá)矩陣列名
head(exp)
colnames(exp) = paste(group_list,1:6,sep='')
head(exp)
# 定義nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
# 聚類
hc=hclust(dist(t(exp)))
par(mar=c(5,5,5,10))
# 繪圖
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
畫出圖后我們發(fā)現(xiàn)，control和treatment很好的分開了，組內(nèi)也很好的聚類到了一起說(shuō)明數(shù)據(jù)過(guò)關(guān)。
PCA
library(ggfortify)
# 互換行和列，再dim一下
df=as.data.frame(t(exp))
# 不要view df，列太多，軟件會(huì)卡??；
dim(df)
dim(exp)
exp[1:6,1:6]
df[1:6,1:6]
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
save(exp,group_list,file = "step2output.Rdata")
同樣發(fā)現(xiàn)該分開的分開了，該聚在一起的聚在一起了，數(shù)據(jù)很好，符合預(yù)期。
經(jīng)過(guò)上面一系列的表達(dá)矩陣可視化，我們檢查了表達(dá)矩陣發(fā)現(xiàn)是正確的，接下來(lái)就要做差異分析啦
4. 差異分析及可視化
芯片數(shù)據(jù)做差異分析最常用的就是limma包
使用這個(gè)包需要三個(gè)數(shù)據(jù)：
- 表達(dá)矩陣(exp)
- 分組矩陣(design)
- 差異比較矩陣(contrast.matrix)
下面我們開始準(zhǔn)備這三個(gè)輸入數(shù)據(jù)：
表達(dá)矩陣(exp)我們?cè)缇偷玫搅?，不用再制作了?br>
我們也得到了存放分組信息的向量group_list，用它來(lái)制作我們的分組矩陣
4.1 limma包做差異分析輸入數(shù)據(jù)的準(zhǔn)備
輸入數(shù)據(jù)—分組矩陣
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "step2output.Rdata")
dim(exp)
library(limma)
# 做分組矩陣
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exp)
design #得到的分組矩陣
輸入數(shù)據(jù)—差異比較矩陣
> contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse = "-"),levels = design)
> contrast.matrix
Contrasts
Levels treat-contorl
contorl -1
treat 1
contrast.matrix 這個(gè)矩陣聲明，我們要把treat組和contor組進(jìn)行差異分析比較：-1和1的意思是contorl是用來(lái)被比的，treat是來(lái)比的即：treat/contorl
到此，做差異分析所需要的三個(gè)矩陣就做好了：表達(dá)矩陣(exp)、分組矩陣(design)、差異比較矩陣(contrast.matrix)
我們已經(jīng)制作好了必要的輸入數(shù)據(jù)，下面開始講如何使用limma包來(lái)進(jìn)行差異分析！
4.2 limma包做差異分析
只有三個(gè)步驟：
- lmFit
- eBayes
- topTable
##step1
fit <- lmFit(exp,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
save(nrDEG,file = "DEGoutput.Rdata")
此時(shí)我們就得到差異分析矩陣（nrDEG），重點(diǎn)看logFC和P值：
差異分析就是對(duì)每個(gè)基因都進(jìn)行檢驗(yàn)，檢驗(yàn)基因的
logFG是多大、平均表達(dá)量是多少、p.value是否顯著等...
4.3 差異表達(dá)基因的可視化
用
limma包得到差異分析表達(dá)矩陣后作圖檢查差異基因是否真的很差異
畫熱圖
選差異最顯著的前25個(gè)基因畫熱圖，查看差異是否真的很顯著
##熱圖
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
load(file = "DEGinput.Rdata")
library(pheatmap)
choose_gene = head(rownames(nrDEG),25)
choose_matrix = exp[choose_gene,]
choose_matrix = t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
火山圖
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
colnames(nrDEG)
plot(nrDEG$logFC,-log10(nrDEG$P.Value))
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
5.富集分析—KEGG、 GO
富集分析就是用常用的數(shù)據(jù)庫(kù)來(lái)注釋基因list。差異分析通過(guò)自定義的閾值挑選了有統(tǒng)計(jì)學(xué)顯著的基因列表后我們其實(shí)是需要對(duì)它們進(jìn)行注釋才能了解其功能，最常見的就是GO/KEGG數(shù)據(jù)庫(kù)注釋，當(dāng)然也可以使用Reactome和Msigdb數(shù)據(jù)庫(kù)來(lái)進(jìn)行注釋。而最常見的注釋方法就是超幾何分布檢驗(yàn)。超幾何分布檢驗(yàn)，運(yùn)用到通路的富集概念就是“總共有多少基因（這個(gè)地方值得注意，主流認(rèn)為只考慮那些在KEGG等數(shù)據(jù)庫(kù)注釋的背景基因），你的通路有多少基因，你的通路被抽中了多少基因（在差異基因里面屬于你的通路的基因）?！?目的就是知道，哪些通路中的哪些基因的表達(dá)因?yàn)樗幬锘蛘吣承┎僮鞯淖饔冒l(fā)生了較大的變化，導(dǎo)致通路有較大改變。
5.1 KEGG pathway analysis
clusterProfiler包作KEGG富集分析
#clusterProfiler作kegg富集分析：
library(clusterProfiler)
gene_up= deg[deg$change == 'up','ENTREZID']
gene_down=deg[deg$change == 'down','ENTREZID']
gene_diff=c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
dim(kk.up)
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
dim(kk.down)
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
class(kk.diff)
#提取出數(shù)據(jù)框
kegg_diff_dt <- kk.diff@result
#根據(jù)pvalue來(lái)選,用于可視化
down_kegg <- kk.down@result %>%
filter(pvalue<0.05) %>%
mutate(group=-1)
up_kegg <- kk.up@result %>%
filter(pvalue<0.05) %>%
mutate(group=1)
#可視化走起
kegg_plot <- function(up_kegg,down_kegg){
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
geom_bar(stat="identity") +
scale_fill_gradient(low="blue",high="red",guide = FALSE) +
scale_x_discrete(name ="Pathway names") +
scale_y_continuous(name ="log10P-value") +
coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
ggtitle("Pathway Enrichment")
}
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
ggsave(g_kegg,filename = 'kegg_up_down.png')
GSEA作KEGG富集分析
GSEA是另一個(gè)常用的富集分析方法，目的是看看基因全局表達(dá)量的變化是否有某些特定的基因集合的傾向性。
data(geneList, package="DOSE")
head(geneList)
length(geneList)
names(geneList)
boxplot(geneList)
boxplot(deg$logFC)
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.9,
verbose = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
gse_kegg=kegg_plot(up_kegg,down_kegg)
print(gse_kegg)
ggsave(gse_kegg,filename ='kegg_up_down_gsea.png')
5.2 GO database analysis
做GO數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解。
GO富集分析生物學(xué)意義
GO富集分析原理：有一個(gè)term注釋了100個(gè)差異表達(dá)基因參與了哪個(gè)過(guò)程，注釋完之后（模式生物都有現(xiàn)成的注釋包，不用我們自己注釋），計(jì)算相對(duì)于背景它是否顯著集中在某條通路、某一個(gè)細(xì)胞學(xué)定位、某一種生物學(xué)功能。
對(duì)GO database analysis一般從三個(gè)層面進(jìn)行：
- Cellular component，CC 細(xì)胞成分
- Biological process， BP 生物學(xué)過(guò)程
- Molecular function，MF 分子功能
這三個(gè)層面具體是指：
- Cellular component解釋的是基因存在在哪里，在細(xì)胞質(zhì)還是在細(xì)胞核？如果存在細(xì)胞質(zhì)那在哪個(gè)細(xì)胞器上？如果是在線粒體中那是存在線粒體膜上還是在線粒體的基質(zhì)當(dāng)中？這些信息都叫Cellular component。
- Biological process是在說(shuō)明該基因參與了哪些生物學(xué)過(guò)程，比如，它參與了rRNA的加工或參與了DNA的復(fù)制，這些信息都叫Biological process
-
Molecular function在講該基因在分子層面的功能是什么？它是催化什么反應(yīng)的？
立足于這三個(gè)方面，我們將得到基因的注釋信息。
GO富集分析的R代碼
#go富集分析--耗費(fèi)時(shí)間灰常長(zhǎng)，很正常
library(clusterProfiler)
#輸入數(shù)據(jù)
gene_up= deg[deg$change == 'up','ENTREZID']
gene_down=deg[deg$change == 'down','ENTREZID']
gene_diff=c(gene_up,gene_down)
head(deg)
#**GO分析三大塊**
#細(xì)胞組分
ego_CC <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#生物過(guò)程
ego_BP <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#分子功能：
ego_MF <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
save(ego_CC,ego_BP,ego_MF,file = "ego_GPL6244.Rdata")
rm(list = ls()) ## 魔幻操作，一鍵清空~
load(file = "ego_GPL6244.Rdata")
#作圖
#第一種，條帶圖，按p從小到大排的
barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")
barplot(ego_BP, showCategory=20,title="EnrichmentGO_CC")
#如果運(yùn)行了沒(méi)出圖，就dev.new()
#第二種，點(diǎn)圖，按富集數(shù)從大到小的
dotplot(ego_CC,title="EnrichmentGO_BP_dot")
#保存
pdf(file = "dotplot_GPL6244.pdf")
dotplot(ego_CC,title="EnrichmentGO_BP_dot")
dev.off()
純代碼版：
#GEO B站視頻 純代碼篇
#下載加載包
cran_packages <- c('tidyr',
'tibble',
'dplyr',
'stringr',
'ggplot2',
'ggpubr',
'factoextra',
'FactoMineR',
'WGCNA')
Biocductor_packages <- c('GEOquery',
'hgu133plus2.db',
"KEGG.db",
"limma",
"impute",
"GSEABase",
"GSVA",
"clusterProfiler",
"genefu",
"org.Hs.eg.db",
"preprocessCore",
"hugene10sttranscriptcluster.db")
for (pkg in c(Biocductor_packages,cran_packages)){
require(pkg,character.only=T)
}
# 下載數(shù)據(jù)
rm(list = ls())
library(GEOquery)
eSet <- getGEO("GSE42872",
destdir = '.',
getGPL = F)
# 從eSet中提取表達(dá)矩陣exp
exp <- exprs(eSet[[1]])
head(exp)
# ID轉(zhuǎn)換
##探針I(yè)D（probe_id）轉(zhuǎn)換成symbol ID
eSet[[1]]@annotation
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids=toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
table(rownames(exp) %in% ids$probe_id)
dim(exp)
exp = exp[rownames(exp) %in% ids$probe_id,]
dim(exp)
ids=ids[match(rownames(exp),ids$probe_id),]
head(ids)
head(exp)
tmp = by(exp,
ids$symbol,
function(x) rownames(x)[which.max(rowMeans(x))])
probes = as.character(tmp)
head(tmp)
head(probes)
dim(exp)
exp = exp[rownames(exp) %in% probes,]
dim(exp)
rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]
head(exp)
pd <- pData(eSet[[1]]) # pData函數(shù)得到每個(gè)樣本的描述信息
head(pd)
save(pd,exp,file = "step1output.Rdata")
save(exp,file = "DEGinput.Rdata")
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "step1output.Rdata")
#####
#檢查表達(dá)矩陣
##畫典型基因表達(dá)量的boxplot
exp['GAPDH',]
exp['ACTB',]
boxplot(exp)
boxplot(exp['GAPDH',])
boxplot(exp['ACTB',])
#各個(gè)樣本表達(dá)量的boxplot
##準(zhǔn)備畫圖所需數(shù)據(jù)exp_L
library(reshape2)
head(exp)
exp_L = melt(exp)
head(exp_L)
colnames(exp_L)=c('symbol','sample','value')
head(exp_L)
# 獲得分組信息
library(stringr)
group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
group_list
exp_L$group = rep(group_list,each=nrow(exp))
head(exp_L)
table(exp_L[,2])
dim(exp_L)
### ggplot2畫圖
library(ggplot2)
p = ggplot(exp_L,
aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
# 對(duì)表達(dá)矩陣進(jìn)行聚類繪圖，并添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)
## hclust
# 更改表達(dá)矩陣列名
head(exp)
colnames(exp) = paste(group_list,1:6,sep='')
head(exp)
# 定義nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
# 聚類
hc=hclust(dist(t(exp)))
par(mar=c(5,5,5,10))
# 繪圖
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
## PCA
library(ggfortify)
# 互換行和列，dim一下
head(exp)
df=as.data.frame(t(exp))
# 不要view df，列太多，軟件會(huì)崩掉；
dim(df)
dim(exp)
exp[1:6,1:6]
df[1:6,1:6]
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
save(exp,group_list,file = "step2output.Rdata")
###################################################
############用limma對(duì)芯片數(shù)據(jù)做差異分析############
###################################################
#差異分析——limma
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "step2output.Rdata")
dim(exp)
library(limma)
# 做分組矩陣
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exp)
design #分組矩陣
# 做比較矩陣
# contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
# contrast.matrix ##這個(gè)矩陣聲明，我們要把treat組和contorl組進(jìn)行差異分析比較
# -1和1的意思是contorl是用來(lái)被比的，treat是來(lái)比的
contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse = "-"),levels = design)
contrast.matrix
#到此，做差異分析所需要的三個(gè)矩陣就做好了：表達(dá)矩陣、分組矩陣、差異比較矩陣
#我們已經(jīng)制作好了必要的輸入數(shù)據(jù)，下面開始講如何使用limma這個(gè)包來(lái)進(jìn)行差異分析了！
##step1
fit <- lmFit(exp,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
save(exp,group_list,nrDEG,file = "DEGoutput.Rdata")
#用limma包得到差異分析表達(dá)矩陣后作圖檢查差異基因是否真的很差異
##熱圖
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
load(file = "DEGinput.Rdata")
library(pheatmap)
choose_gene = head(rownames(nrDEG),25)
choose_matrix = exp[choose_gene,]
choose_matrix = t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
##火山圖
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
colnames(nrDEG)
plot(nrDEG$logFC,-log10(nrDEG$P.Value))
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
#富集分析
#富集分析準(zhǔn)備工作：
##首先對(duì)差異表達(dá)矩陣nrDEG，進(jìn)行加工
###1.把行名變成SYMBOL列
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
library(dplyr)
deg = nrDEG
deg <- mutate(deg,symbol = rownames(deg))
head(deg)
###2.加change列：上調(diào)或下調(diào)，火山圖要用
logFC_t = 1 #不同的閾值，篩選到的差異基因數(shù)量就不一樣，后面的超幾何分布檢驗(yàn)結(jié)果就大相徑庭。
change=ifelse(deg$P.Value>0.01,'stable',
ifelse( deg$logFC >logFC_t,'up',
ifelse( deg$logFC < -logFC_t,'down','stable') )
)
deg <- mutate(deg,change)
head(deg)
table(deg$change)
###3.加ENTREZID列，后面富集分析要用
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(unique(deg$symbol), fromType = "SYMBOL", #ID轉(zhuǎn)換核心函數(shù)bitr
toType = c( "ENTREZID"),
OrgDb = org.Hs.eg.db)
head(s2e)
head(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
head(deg)
save(exp,group_list,deg,file = "enrich_input.Rdata")
#####################
######富集分析#######
#####################
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'enrich_input.Rdata')
## 1.KEGG pathway analysis
#上調(diào)、下調(diào)、差異、所有基因
#clusterProfiler作kegg富集分析：
library(clusterProfiler)
gene_up= deg[deg$change == 'up','ENTREZID']
gene_down=deg[deg$change == 'down','ENTREZID']
gene_diff=c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
dim(kk.up)
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
dim(kk.down)
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
class(kk.diff)
#提取出數(shù)據(jù)框
kegg_diff_dt <- kk.diff@result
#根據(jù)pvalue來(lái)選,用于可視化
down_kegg <- kk.down@result %>%
filter(pvalue<0.05) %>%
mutate(group=-1)
up_kegg <- kk.up@result %>%
filter(pvalue<0.05) %>%
mutate(group=1)
#可視化
kegg_plot <- function(up_kegg,down_kegg){
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
geom_bar(stat="identity") +
scale_fill_gradient(low="blue",high="red",guide = FALSE) +
scale_x_discrete(name ="Pathway names") +
scale_y_continuous(name ="log10P-value") +
coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
ggtitle("Pathway Enrichment")
}
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
ggsave(g_kegg,filename = 'kegg_up_down.png')
#gsea作kegg富集分析：
data(geneList, package="DOSE")
head(geneList)
length(geneList)
names(geneList)
boxplot(geneList)
boxplot(deg$logFC)
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.9,
verbose = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
gse_kegg=kegg_plot(up_kegg,down_kegg)
print(gse_kegg)
ggsave(gse_kegg,filename ='kegg_up_down_gsea.png')
### 2.GO database analysis
#go富集分析
library(clusterProfiler)
#輸入數(shù)據(jù)
gene_up= deg[deg$change == 'up','ENTREZID']
gene_down=deg[deg$change == 'down','ENTREZID']
gene_diff=c(gene_up,gene_down)
head(deg)
#**GO分析三大塊**
#細(xì)胞組分
ego_CC <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#生物過(guò)程
ego_BP <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#分子功能：
ego_MF <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
save(ego_CC,ego_BP,ego_MF,file = "ego_GPL6244.Rdata")
rm(list = ls())
load(file = "ego_GPL6244.Rdata")
#第一種，條帶圖，按p從小到大排的
barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")
barplot(ego_BP, showCategory=20,title="EnrichmentGO_CC")
#如果運(yùn)行了沒(méi)出圖，就dev.new()
#第二種，點(diǎn)圖，按富集數(shù)從大到小的
dotplot(ego_CC,title="EnrichmentGO_BP_dot")
#保存
pdf(file = "dotplot_GPL6244.pdf")
dotplot(ego_CC,title="EnrichmentGO_BP_dot")
dev.off()
特別感謝小潔老師激發(fā)了我學(xué)習(xí)GEO數(shù)據(jù)庫(kù)挖掘的興趣；有些圖片還有富集分析的代碼就來(lái)自小潔老師的課件哦
