【原】不同的GSE數(shù)據(jù)集有不同的臨床信息，不同的分組技巧

健明 2021-07-14

展開全文

最近，我發(fā)現(xiàn)學(xué)徒在學(xué)習(xí)GEO數(shù)據(jù)挖掘的過程中，遇到了第一個(gè)也是至關(guān)重要的一個(gè)難題就是對下載后的數(shù)據(jù)集進(jìn)行合適的分組，因?yàn)橹挥袑颖具M(jìn)行合適的分組，才有可能得到我們想要的信息。但是不同的GSE數(shù)據(jù)集有不同的臨床信息，那么我們應(yīng)該挑選合適的臨床信息來進(jìn)行分組呢？
這里面涉及到兩個(gè)問題，首先是能否看懂?dāng)?shù)據(jù)集配套的文章，從而達(dá)到正確的生物學(xué)意義的分組，其次能否通過R代碼實(shí)現(xiàn)這個(gè)分組。同樣的我也是安排學(xué)徒完成了部分任務(wù)并且總結(jié)出來了！
下面看學(xué)徒的表演（PS: 圖片較多的推文，排版真的是嚇?biāo)廊耍。?/p>

Jimmy大神怎么說過，只有多做、多錯(cuò)，才能真正的掌握。所以下面通過幾個(gè)實(shí)戰(zhàn)來說明。

首先是通過對一篇文獻(xiàn)Identification of potential core genes in triple negative breast cancer using bioinformatics analysis所用到的三個(gè)TNBC（Triple-Negative Breast Cancer）三陰性乳腺癌的三個(gè)數(shù)據(jù)集：GSE38959、GSE45827以及GSE62194進(jìn)行分組，首先對GSE38959進(jìn)行分析。

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動(dòng)化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
gset <- getGEO('GSE38959', destdir=".",
              AnnotGPL = F,     ## 注釋文件
              getGPL = F)       ## 平臺文件

class(gset)  #查看數(shù)據(jù)類型
length(gset)  #
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]] #
dat=exprs(a) #a現(xiàn)在是一個(gè)對象，取a這個(gè)對象通過看說明書知道要用exprs這個(gè)函數(shù)
dim(dat)#看一下dat這個(gè)矩陣的維度
dat[1:4,1:4] #查看dat這個(gè)矩陣的1至4行和1至4列，逗號前為行，逗號后為列
pd=pData(a) #通過查看說明書知道取對象a里的臨床信息用pData

pd就是這個(gè)數(shù)據(jù)集的臨床信息，查看后如下

會(huì)發(fā)現(xiàn)有些信息是冗余的，有些是有效信息可以用來分組，但是表型記錄太多，看起來會(huì)混淆，所以需要去除那些冗余信息，就是在所有樣本里面表型記錄都一致的列。如何去冗余,見原文對表型數(shù)據(jù)框進(jìn)行去冗余。

閱讀文章后發(fā)現(xiàn)原文把樣本分為2組：腫瘤與正常，而且總共只有43個(gè)樣本，而臨床信息有47個(gè)樣本，說明有效信息列含有3個(gè)或3個(gè)以上元素，可以再縮小范圍。（注意！如果樣本數(shù)剛好去冗余就行！）

pd1=pd[,apply(pd, 2, function(x){
length(unique(x))>2})]  #縮小范圍
dim(pd1)
#[1] 47  9

發(fā)現(xiàn)范圍已經(jīng)縮小很多。對數(shù)據(jù)框再用apply循環(huán)去查找文章作者是用哪一列來分組的

apply(pd1,2,table)

通過循環(huán)，就可以清楚的知道該用哪一列來進(jìn)行分組啦

然后是搜索關(guān)鍵字進(jìn)行分組

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)] #對表達(dá)矩陣取子集
group_list=c(rep('TNBC',length(TNBC)) ,
         rep('NOR',length(NOR))) #分組信息
table(group_list)
#group_list
#NOR TNBC 
#13   30

第二個(gè)數(shù)據(jù)集GSE45827同樣的方法，重復(fù)的地方不贅述，從有差異的地方開始。

文中為了2組，用了52個(gè)樣本。

然后是搜索關(guān)鍵字進(jìn)行分組

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

第三個(gè)數(shù)據(jù)集同理

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三個(gè)數(shù)據(jù)集

這里主要說一下GSE31056這一個(gè)數(shù)據(jù)集，需要一定的背景知識與細(xì)心才能正常分組，原文里

如果用我們之前的方法找是找不到的，因?yàn)榧?xì)心點(diǎn)你會(huì)發(fā)現(xiàn)GSE給的位置不止tongue，還有mouth等，而文章只需要tongue。所以我們需要對數(shù)據(jù)集取子集。

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)

根據(jù)生物學(xué)相關(guān)知識，margin其實(shí)就屬于normal，所以就找出了分組。

所以可以看到生物學(xué)知識多么重要：沒有生物學(xué)背景的數(shù)據(jù)分析很危險(xiǎn)

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.分析關(guān)于ccRCC (clear cell renal cell carcinoma)腎透明細(xì)胞癌的一個(gè)GSE子集GSE53757

下載數(shù)據(jù)、提取表達(dá)矩陣與臨床信息方法與前面一直，這里就不贅述，也是從有差異的地方開始。這里文章作者只需要第三期的正常與腫瘤組的樣本

table(pd$`tissue:ch1`)
table(pd$`tumor stage:ch1`)
table(pd$source_name_ch1)

通過table函數(shù)，我們看到總共144個(gè)樣本，其中有72個(gè)正常與72個(gè)腫瘤樣本；第三期腫瘤和正常樣本總各有14個(gè)，下面我們就需要提取我們需要的數(shù)據(jù)

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]) #表達(dá)矩陣的列名就是臨床信息的第一列
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