不同的GSE數據集有不同的臨床資訊，不同的分組技巧
- 2019 年 11 月 11 日
- 筆記
最近，我發現學徒在學習GEO數據挖掘的過程中，遇到了第一個也是至關重要的一個難題就是對下載後的數據集進行合適的分組，因為只有對樣本進行合適的分組，才有可能得到我們想要的資訊。但是不同的GSE數據集有不同的臨床資訊，那麼我們應該挑選合適的臨床資訊來進行分組呢？ 這裡面涉及到兩個問題，首先是能否看懂數據集配套的文章，從而達到正確的生物學意義的分組，其次能否通過R程式碼實現這個分組。同樣的我也是安排學徒完成了部分任務並且總結出來了！ 下面看學徒的表演（PS: 圖片較多的推文，排版真的是嚇死人！）
Jimmy大神怎麼說過，只有多做、多錯，才能真正的掌握。所以下面通過幾個實戰來說明。
- 首先是通過對一篇文獻Identification of potential core genes in triple negative breast cancer using bioinformatics analysis所用到的三個TNBC（Triple-Negative Breast Cancer）三陰性乳腺癌的三個數據集：GSE38959、GSE45827以及GSE62194進行分組，首先對GSE38959進行分析。
library(GEOquery) # 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。 #Setting options('download.file.method.GEOquery'='auto') #Setting options('GEOquery.inmemory.gpl'=FALSE) gset <- getGEO('GSE38959', destdir=".", AnnotGPL = F, ## 注釋文件 getGPL = F) ## 平台文件 class(gset) #查看數據類型 length(gset) # class(gset[[1]]) # 因為這個GEO數據集只有一個GPL平台，所以下載到的是一個含有一個元素的list a=gset[[1]] # dat=exprs(a) #a現在是一個對象，取a這個對象通過看說明書知道要用exprs這個函數 dim(dat)#看一下dat這個矩陣的維度 dat[1:4,1:4] #查看dat這個矩陣的1至4行和1至4列，逗號前為行，逗號後為列 pd=pData(a) #通過查看說明書知道取對象a里的臨床資訊用pData
pd就是這個數據集的臨床資訊，查看後如下
會發現有些資訊是冗餘的，有些是有效資訊可以用來分組，但是表型記錄太多，看起來會混淆，所以需要去除那些冗餘資訊，就是在所有樣本裡面表型記錄都一致的列。如何去冗餘,見原文對錶型數據框進行去冗餘。
閱讀文章後發現原文把樣本分為2組：腫瘤與正常，而且總共只有43個樣本，而臨床資訊有47個樣本，說明有效資訊列含有3個或3個以上元素，可以再縮小範圍。（注意！如果樣本數剛好去冗餘就行！）
pd1=pd[,apply(pd, 2, function(x){ length(unique(x))>2})] #縮小範圍 dim(pd1) #[1] 47 9
發現範圍已經縮小很多。對數據框再用apply循環去查找文章作者是用哪一列來分組的
apply(pd1,2,table)
通過循環，就可以清楚的知道該用哪一列來進行分組啦
然後是搜索關鍵字進行分組
TNBC=rownames(pd1[grepl('triple negative breast cancer cells',as.character(pd$`cell type:ch1`)),]) #腫瘤組 NOR=rownames(pd1[grepl('mammary gland ductal cells',as.character(pd$`cell type:ch1`)),]) #正常組 dat=dat[,c(TNBC,NOR)] #對表達矩陣取子集 group_list=c(rep('TNBC',length(TNBC)) , rep('NOR',length(NOR))) #分組資訊 table(group_list) #group_list #NOR TNBC #13 30
第二個數據集GSE45827同樣的方法，重複的地方不贅述，從有差異的地方開始。
文中為了2組，用了52個樣本。
然後是搜索關鍵字進行分組
TNBC=rownames(pd[grepl('Human Basal Tumor Sample',as.character(pd$source_name_ch1)),]) #腫瘤組 NOR=rownames(pd[grepl('Human Normal',as.character(pd$source_name_ch1)),]) #正常組 dat=dat[,c(TNBC,NOR)] group_list=c(rep('TNBC',length(TNBC)), rep('NOR',length(NOR))) #分組資訊 table(group_list) #group_list #NOR TNBC # 11 41
第三個數據集同理
2.選取OTSCC(Oral Tongue Squamous Cell Carcinoma)口腔舌鱗癌Integrative analysis of gene expression profiles reveals distinct molecular characteristics in oral tongue squamous cell carcinoma的GSE13601, GSE31056 and GSE78060三個數據集
這裡主要說一下GSE31056這一個數據集，需要一定的背景知識與細心才能正常分組，原文里
如果用我們之前的方法找是找不到的，因為細心點你會發現GSE給的位置不止tongue，還有mouth等，而文章只需要tongue。所以我們需要對數據集取子集。
pd1=pd[,apply(pd, 2, function(x){ length(unique(x))>1})] pd1=pd1[grepl('tongue',as.character(pd$characteristics_ch1.2)),] apply(pd1,2,table)
根據生物學相關知識，margin其實就屬於normal，所以就找出了分組。
所以可以看到生物學知識多麼重要：沒有生物學背景的數據分析很危險
TU=rownames(pd1[grepl('tumor',as.character(pd1$`site:ch1`)),]) #腫瘤 NOR=rownames(pd1[grepl('margin',as.character(pd1$`site:ch1`)),])#正常 dat=dat[,c(TU,NOR)]#取子集 group_list=c(rep('TU',length(TU)), #分組資訊 rep('NOR',length(NOR))) table(group_list) #group_list #NOR TU # 39 12
3.分析關於ccRCC (clear cell renal cell carcinoma)腎透明細胞癌的一個GSE子集GSE53757
下載數據、提取表達矩陣與臨床資訊方法與前面一直，這裡就不贅述，也是從有差異的地方開始。這裡文章作者只需要第三期的正常與腫瘤組的樣本
table(pd$`tissue:ch1`) table(pd$`tumor stage:ch1`) table(pd$source_name_ch1)
通過table函數，我們看到總共144個樣本，其中有72個正常與72個腫瘤樣本；第三期腫瘤和正常樣本總各有14個，下面我們就需要提取我們需要的數據
patient_t = pd[pd$`tissue:ch1`=='clear cell renal cell carcinoma' & pd$`tumor stage:ch1`=='stage 3', ] #通過向量取交集的方式來取 patient_n=pd[pd$source_name_ch1=='normal match to Stage 3 ccRCC',] colnames(dat)=as.character(pd[,1]) #表達矩陣的列名就是臨床資訊的第一列 dat=dat[,c(patient_t[,1],patient_n[,1])]#取子集 group_list=rep(c('ccRCC','normal'),each=14) #分組資訊 table(group_list) #group_list #ccRCC normal #14 14
總結一下，我們可以根據自己的需求選取合適的程式碼去進行有效的分組，在不同的情況下選取最合適當下的方法，方便自己去做後續的數據分析。
