來自數(shù)據(jù)集GSE33532,對于第一步獲得表達(dá)矩陣早抠，做了還蠻詳盡的注解過程，希望能對其他人有幫助撬讽。雖然看上去全是黑色代碼快蕊连，但注解理解示例全在這里面了哦！

rm(list = ls())  ## 魔幻操作游昼，一鍵清空~
options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(stringsAsFactors = F)#在調(diào)用as.data.frame的時(shí)甘苍，將stringsAsFactors設(shè)置為FALSE可以避免character類型自動轉(zhuǎn)化為factor類型


##################     1. 下載表達(dá)矩陣和臨床信息         ###################
##################  得到表達(dá)矩陣dat和分組信息group_lis   ##################
f='GSE33532_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33532
library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的烘豌。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE33532', destdir=".",
                 AnnotGPL = F,     ## 注釋文件
                 getGPL = F)       ## 平臺文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE33532_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)  #查看數(shù)據(jù)類型
length(gset)  #
class(gset[[1]])
gset
# 因?yàn)檫@個GEO數(shù)據(jù)集只有一個GPL平臺载庭，所以下載到的是一個含有一個元素的list
a=gset[[1]] #
dat=exprs(a) #a現(xiàn)在是一個對象，取a這個對象通過看說明書知道要用exprs這個函數(shù)
dim(dat) 

dat[1:4,1:4] #查看dat這個矩陣的1至4行和1至4列廊佩，逗號前為行囚聚，逗號后為列
pd=pData(a)



tumor=rownames(pd[grepl('tumor',as.character(pd$title)),])
normal=rownames(pd[grepl('normal',as.character(pd$title)),])

dat=dat[,c(tumor,normal)]
group_list=c(rep('tumor',length(tumor)),
             rep('normal',length(normal)))
table(group_list)
dim(dat)
dat[1:4,1:4] 
dat <- as.data.frame(dat)
save(dat,group_list,file = 'step1.1_dat_group.Rdata')




##################   2.獲得注釋信息   ##################
##################       三種方法     ##################
#1.用getGEO下載平臺GPL的注釋信息
#2.在GEO網(wǎng)站上下載好平臺注釋信息后倒入R studio
#3.用bioconductor上的包


#方法一 用getGEO下載平臺GPL的注釋信息
if(F){
  library(GEOquery)
  gpl<- getGEO('GPL570', destdir=".")
  dim(gpl)
  colnames(Table(gpl)) #查一下列明
  head(Table(gpl)[,c(1,11)]) 
  ids=Table(gpl)[,c(1,11)]
  ##踩坑放松1:會不會有NA值，檢查一下,沒問題就不管那些個警告了
  tmpids <- na.omit(ids)
  dim(tmpids)
  head(ids)
  save(ids,file='step1.2_ids.Rdata') #三種方法均獲得的ids均設(shè)置的是'step1.2_ids.Rdata'
}
rm(tmpids)


#方法二 在GEO網(wǎng)站上下載好平臺注釋信息后導(dǎo)入Rstudio
if(F){
  gpl=read.table('GPL570-55999.txt', #提前下載好了平臺GPL570的注釋信息
               sep = '\t',quote = '',
               fill=T,header = T)
  dim(gpl)
  colnames(gpl)
  ids=gpl[,c(1,11)]
  ids=ids[ids$Gene.Symbol != '',] #去除空的
  dim(ids)
  head(ids) 
  save(ids,file='step1.2_ids.Rdata')  #三種方法均獲得的ids均設(shè)置的是'step1.2_ids.Rdata'
}


#方法三 用bioconductor網(wǎng)站上對應(yīng)平臺的基因注釋包
#https://mp.weixin.qq.com/s/5_ceRNaRFP3h08Pjd015Wg
if(T){
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db")
  ids<- toTable(hgu133plus2SYMBOL)
  dim(ids)
  head(ids)
  save(ids,file='step1.2_ids.Rdata')  #三種方法均獲得的ids均設(shè)置的是'step1.2_ids.Rdata'
}

################三種獲得注釋信息的方法标锄，所得到的對應(yīng)的行數(shù)如下：
#> dim(gpl)
#[1] 54675    16
#> dim(ids)
#[1] 45782     2
#> dim(ids)
#[1] 41932     2



##################        3. ID轉(zhuǎn)換         ##################
##################   intersect/merge/%in%   ##################
##################         match函數(shù)        ##################
##################  （1）用intersect函數(shù)    ##################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata') #此處按照第三種從bioconductor獲得的注釋包來進(jìn)行id轉(zhuǎn)換
dim(dat)
dim(ids)
table(rownames(dat) %in% ids$probe_id)
intersect(rownames(dat),ids$probe_id) 
length(intersect(rownames(dat),ids$probe_id))
#xt意味著xt(相同)的
xt <- intersect(rownames(dat),ids$probe_id)
dim(xt) 
##踩坑放松1:為什么用dim是得到NULL?

length(xt)######xt現(xiàn)在是一個character顽铸，所以用length來查看向量的長度
head(xt)
ids_tmp <- ids[ids$probe_id==xt,]#注：想將ids按照取交集后的向量（xt）從ids中取交集
##踩坑放松2:報(bào)錯-長的對象長度不是短的對象長度的整倍數(shù)

length(ids$probe_id)
length(xt)

###tips1：上面會提示報(bào)錯，因?yàn)楫?dāng)向量用'=='進(jìn)行判斷的時(shí)候料皇，就是將兩個向量中的元素進(jìn)行一對一地判斷是否相同谓松，之所以報(bào)‘長的對象長度不是短的對象長度的整倍數(shù)’的報(bào)錯，是因?yàn)閕ds$probe_id有ids的行數(shù)那么多（41932個）践剂，而xt則是取交集后的23738個鬼譬。
###tips2：接下來用match函數(shù)取子集，就是按照xt的順序從大子集ids中取舷手，那么取出來的探針的順序肯定是和xt的順序是一樣的拧簸，然后再從dat中同樣按照xt的順序來取子集，因?yàn)閕ds和dat都是按照xt中的探針的順序取子集男窟，那么二者的dat中的探針名肯定和ids中的探針名是一一對應(yīng)的了盆赤。


####################  重磅嘉賓-match函數(shù)第一次出場  ####################

ids <- ids[match(xt,ids$probe_id),]##ids同樣按照probe_id這一列提取子集贾富，是按照坐標(biāo)位置提取
#用identical函數(shù)和setequal函數(shù)驗(yàn)證一下(這兩個函數(shù)都可以)，是不是完全一樣牺六，也就是說表達(dá)矩陣dat的探針名和交集（xt）是否是一一對應(yīng)的
identical(ids$probe_id,xt)
setequal(ids$probe_id,xt)
save(ids,file='step1.3_ids.Rdata')
  

####################  （2）用merge函數(shù)    ##################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata')
ids<- toTable(hgu133plus2SYMBOL)
dat_tmp <- dat #重新賦值一個dat_tmp
ids_tmp <- ids #重新賦值一個ids_tmp
dim(dat_tmp)
dim(ids_tmp)
#用新賦值后的做探索
dat_tmp <- dat_tmp[ids_tmp$probe_id,] #注：想從dat_tmp表達(dá)矩陣中按照ids_tmp$probe_id取出的探針名提取子集颤枪，dat_tmp的rownames就是探針名，表面看是沒問題呀

##踩坑放松3: dat_tmp <- dat_tmp[ids_tmp$probe_id,]無報(bào)錯淑际，表面看似風(fēng)平浪靜畏纲，背后其實(shí)翻江倒海，用dim查看大小
dim(dat_tmp)  #注：剛才的dat_tmp只有25906個探針春缕，現(xiàn)在竟然和ids_tmp的探針名一樣多的個數(shù)盗胀！表達(dá)矩陣竟然變大了！
is.na(dat_tmp) #注：我就想看個行名是不是NA锄贼，結(jié)果在表達(dá)矩陣中返回了TRUE和FALSE票灰，所以又學(xué)一了is.na函數(shù)
table(is.na(rownames(dat_tmp))) #注：還是不對，改成行名也不對宅荤。
dat_tmp[1:4,1:4]
is.na(dat_tmp)#注：觀察在dat_tmp中為TRUE也就是為NA值的屑迂，行名也都是NA,因此用grepl函數(shù)抓取下
table(grepl('NA',rownames(dat_tmp))) #注意區(qū)分grep和grepl函數(shù)的差別。注意此時(shí)的結(jié)果
#FALSE  TRUE 
#23738 18194
#注意觀察：上面的FALSE正好為23738個冯键，同前面用intersect函數(shù)獲得的ids是不是行數(shù)相同惹盼？
dat_tmp <- na.omit(dat_tmp) #注：用na.omit去掉NA值
dim(dat_tmp) 

##為用merge函數(shù)作準(zhǔn)備-dat新建一列，列名是probe
dat$probe_id <- rownames(dat) #注：dat新建一列惫确，列名是probe手报，此步是為了后面用merge函數(shù)
dim(dat)
####小提問：如何知道dat中第101列的列名
colnames(dat[,101]) #注：dat新建probe這一列后，dat現(xiàn)在有101列雕薪，我想看看這第101例的列名




##踩坑放松4:colnames(dat[,101])為什么返回NULL,如何獲得第101列的列名昧诱？
#分解代碼：先看dat[,101]
dat[,101]

##報(bào)錯原因：把上面的代碼分解晓淀，dat[,101]所袁，根據(jù)“行在左，列在右”凶掰，此時(shí)101在“燥爷，”右邊，因此其實(shí)取出的是第101這列的所有的探針名懦窘，是一個向量前翎，一個長度同行數(shù)相等的向量。向量取子集怎么能用colnanmes來提取呢畅涂？
#正確操縱如下：
colnames(dat)
colnames(dat)[101]#注港华，從向量中取子集，可以按照邏輯值和坐標(biāo)位置午衰，現(xiàn)在我們想取坐標(biāo)位置為第101位的內(nèi)容


#關(guān)于merge函數(shù)立宜，依然用前兩步得到的dat和ids
ids_dat_merge <- merge(ids,dat,by.x ='probe_id',by.y='probe_id')

##踩坑放松5:
#(1)merge函數(shù)：dat新建的一列冒萄，如果列名是probe，也就是說列名不一樣橙数，可不可以尊流？
#(2)merge函數(shù)：by.x=probe,by.y=probe_id,怎么改？出來的矩陣是什么樣的灯帮？
#(3)merge函數(shù)：如果想要merge的那一列的向量順序不一樣崖技，merge后的矩陣是怎樣的,有什么關(guān)系？


#下面是merge函數(shù)小理解钟哥，新建模擬ids_tmp和dat_tmp
if(F){
  #第1迎献、2個merge問題
  ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
                        symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
  dat_tmp <- data.frame(row.names = c('1053_at','117_at','121_at','1316_at','1320_at'),
                        GSM835268=c(1,1,1,1,1),
                        GSM835269=c(2,2,2,2,2),
                        probe = c('1053_at','117_at','121_at','1316_at','1320_at'))
  ids_tmp
  dat_tmp
  merge1 <- merge(ids_tmp,dat_tmp,by.x='probe_id',by.y='probe')
  merge1
  merge2 <- merge(dat_tmp,ids_tmp,by.x='probe',by.y='probe_id')
  merge2
  
  #第3個merge問題，如果要merge的探針名順序不一樣
  ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
                        symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
  dat_tmp <- data.frame(row.names = c('121_at','117_at','1053_at','1316_at','1320_at'),
                        GSM835268=c(1,1,1,1,1),
                        GSM835269=c(2,2,2,2,2),
                        probe = c('121_at','117_at','1053_at','1316_at','1320_at'))
  ids_tmp
  dat_tmp
  merge3 <- merge(ids_tmp,dat_tmp,by.x='probe_id',by.y='probe')
  merge3
  merge4 <- merge(dat_tmp,ids_tmp,by.x='probe',by.y='probe_id')
  merge4
  #從merge3腻贰、merge4相同的那列順序沒有改變忿晕，這是最后的矩陣中dat和ids順序有差別
}



dim(ids_dat_merge) #注：檢查數(shù)據(jù)，dim看大小
ids_dat_merge[1:4,1:4] #注：取子集看數(shù)據(jù)情況
ids_dat_merge <- na.omit(ids_dat_merge) #注：檢查下是否有NA值银受，注意理解merge函數(shù)践盼，通過merge函數(shù)獲得ids_dat_merge是ids和dat中*都都都*存在的元素。思考宾巍，如果想得到merge(x,y)后咕幻，在x中有不存在與y中的元素也出現(xiàn)在merge后的矩陣中，用什么參數(shù)顶霞？
dim(ids_dat_merge)
ids_dat_merge####用$符號看一下肄程，現(xiàn)在的表達(dá)矩陣有哪些列（列名）按完$健，按tab健

##接下來進(jìn)行從ids中按照merge后的ids_dat_merge表達(dá)矩陣中按照探針名字选浑，從ids中取子集
ids_tmp <- ids[ids_dat_merge$probe_id,] #注：我已經(jīng)得到了dat和ids共有列探針（probe）merge后的表達(dá)矩陣了蓝厌，想按照merge后的ids_dat_merge中的probe_id這一列提取探針名字，然后將ids按照這個探針名字提取子集從而想生成一個ids_tmp古徒，從而獲得和表達(dá)表達(dá)矩陣的探針名相同的ids注釋拓提！
dim(ids_tmp)
ids_tmp[1:4,1:2]#注：一看數(shù)據(jù)具體情況，竟然都是NA!

#踩坑放松6：從數(shù)據(jù)框中提取子集隧膘，可以按rownames也就是行名代态，也可以按照坐標(biāo)，可是為什么都是NA呢疹吃？


dim(ids)
ids[1:4,1:2]
#答：ids是行名為1蹦疑、2、3萨驶、4...歉摧，并不是探針名！從ids中取子集時(shí)，要么按照探針名(行名)叁温，要么按照坐標(biāo)位置,要么按照邏輯值豆挽，就這三種置济，三種音羞，三種！重要的事情說三遍哦馅扣！而ids_dat_merge$probe_id后得到的是探針名probe_id锰镀，而ids的行名又不是探針名probe_id娘侍，而是1、2泳炉、3憾筏、4...,所以此處取出的都是NA!

#既然不是，那我把ids的行名變成探針名probe_id不就好了嘛花鹅！
rownames(ids_tmp) <- ids_dat_merge$probe_id 
##看下大小
dim(ids_tmp)
ids_tmp[1:4,1:2] #注：雖然重新復(fù)制了探針名氧腰，但是需要注意此處的ids_tmp由于上一步取子集的錯誤操作，現(xiàn)在都是NA,所以需要重新獲得一個ids_tmp.


####################  重磅嘉賓-match函數(shù)第二次出場  ####################

ids_tmp <- ids[match(ids_dat_merge$probe_id,ids$probe_id),] 
#注：上面的match函數(shù)刨肃，從ids$probe_id的探針名中按照ids_dat_merge$probe_id在ids$probe_id中的排名位置古拴，也就是坐標(biāo)位置，從ids中按照這個坐標(biāo)位置取子集真友，得到ids_tmp
######下面是match函數(shù)小理解
if(T){
  use_id <- c('A','B','C','D')
  use_id
  m<-data.frame(u=c('B','C','D','A'))
  match(use_id,m$u)
  a<-m[match(use_id,m$u),]
  a
}

dim(ids_tmp)
ids_tmp[1:4,1:2]
ids_tmp <- na.omit(ids_tmp)#注：被NA弄怕了黄痪，再用na.omit檢查下，match函數(shù)用熟練或用對了盔然，就可省略
ids <- ids_tmp
save(ids,file='step1.3_ids.Rdata')

######################    此處小貼士    ###################
#其實(shí)前面的到的ids_dat_merge矩陣中已經(jīng)含有基因symbol這一列了桅打，是不是可以就把矩陣中的symbol變成行名就可以了呢？不就已經(jīng)完成id轉(zhuǎn)換了嘛愈案？那試一下

ids_dat_merge[1:4,1:4]
rownames(ids_dat_merge) <- ids_dat_merge$symbol

##踩坑放松8:rownames(ids_dat_merge) <- ids_dat_merge$symbol 報(bào)錯挺尾，不允許有重復(fù)的‘row.names’,因此基因名symbol需要去除重復(fù)！


####################  （3）用%in%函數(shù)    ##################

rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata')
dim(dat)
dim(ids)
ids <- ids[ids$probe_id %in% rownames(dat),]#注：想從dat的rownames中提取出ids$symbol后存在于rownames(dat)中的元素站绪，%in%是進(jìn)行判斷遭铺，邏輯值為TRUE的，從ids中提取出來崇众。
dim(ids)
ids[1:4,1:2]#注：此時(shí)為什么是ids[1:4,1:2]掂僵，如果用ids[1:4,1:4]會報(bào)錯undefined columns selected航厚，是因?yàn)閕ds本身就只有兩列顷歌。
##踩坑放松9:此處是溫馨小提示：同前面踩坑放松6是類似的道理，雖然現(xiàn)在ids的行名是1幔睬、2眯漩、3、4，但是%in%是進(jìn)行判斷赦抖，得到的是邏輯值舱卡，所以不用管行名是什么，我們就是按照TRUE的提榷佑轮锥！三種提取子集的方法：名字，坐標(biāo)要尔，邏輯值舍杜！


#下面是%in%函數(shù)小理解,新建模擬ids_tmp和dat_tmp.
if(F){
  ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
                  symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
  dat_tmp <- data.frame(row.names = c('1053_at','117_at','121_at','1316_at','1320_at'),
                        GSM835268=c(1,1,1,1,1),
                        GSM835269=c(2,2,2,2,2))
  ids_tmp
  dat_tmp
  m <- ids_tmp[ids_tmp$probe_id %in% rownames(dat_tmp),]
  m
  n <- dat_tmp[rownames(dat_tmp) %in% ids_tmp$probe_id,]
  n
  #如果探針名順序不一樣
  ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
                        symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
  dat_tmp <- data.frame(row.names = c('121_at','1316_at','1053_at','117_at','1320_at'),
                        GSM835268=c(1,1,1,1,1),
                        GSM835269=c(2,2,2,2,2))
  #注意上面dat_tmp中的row.names的順序。
  ids_tmp
  dat_tmp
  m <- ids_tmp[ids_tmp$probe_id %in% rownames(dat_tmp),]
  m
  n <- dat_tmp[rownames(dat_tmp) %in% ids_tmp$probe_id,]
  n
  #上面探針名順序不一樣赵辕，可以幫助理解取子集既绩，從哪個矩陣取，就按哪個矩陣中的邏輯值取
}


identical(ids$probe_id,rownames(dat))
ids[1:4,1:2]
dat[1:4,1:4]

####################  重磅嘉賓-match函數(shù)第三次出場  ####################

dat <- dat[match(ids$probe_id,rownames(dat)),]#注：ids的probe_id是全部存在于dat的rownames里的,因此用match函數(shù)從dat中按照ids的順序提取子集还惠，一定是能得到探針完全對應(yīng)的饲握。此處就是按照坐標(biāo)在提取子集，這個坐標(biāo)指的就是a在b中的順序蚕键，注意救欧，如果數(shù)據(jù)框的行名是1、2锣光、3颜矿、4，與這個1嫉晶、2骑疆、3、4無關(guān)替废，僅僅就是a中的元素一個一個試驗(yàn)在b中排第幾位箍铭，比如a中的元素依次在b中的元素進(jìn)行一一對應(yīng)，排第1位記個1椎镣，排第5位記個5诈火，排第10位記個10，然后從b中按照這個1状答、5冷守、10提取b中的元素。所以從b中提取的元素就是和a完全相同的惊科。

save(ids,file = 'step1.3_ids.Rdata')#注拍摇，上面的dat取子集是按照坐標(biāo)位置提取的，此處并沒有保存dat馆截，但是要知道此出dat用match函數(shù)提取子集的方法充活。


#######################  4.將探針名轉(zhuǎn)換為基因名  ################
######################    基因名symbol去重復(fù) ######################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load(file='step1.3_ids.Rdata')
dim(dat)
dat[1:4,1:4]
dim(ids)
ids[1:4,1:2]
dat <- dat[ids$probe_id,]#注：從dat表達(dá)矩陣（其實(shí)是數(shù)據(jù)框蜂莉，就是這么叫一下表達(dá)矩陣）中按照ids$probe_id也就是探針名字，提取子集混卵。

######下面是數(shù)據(jù)框提取子集的小理解映穗，按照rownames的名字或者坐標(biāo)位置均可以提取子集，示例如下：
  if(F){ 
    x <- data.frame(row.names = c('a','b','c'),xm = c('xh','xl','xh'),nl=c(15,18,21))
    x
    x['a',]
    x[1,]
    identical(x['a',],x[1,])
  }

dim(dat)
dim(ids)
ids[1:4,1:2]
#同理幕随，看看dat和ids能不能一一對應(yīng)
identical(rownames(dat),ids$probe_id)
setequal(rownames(dat),ids$probe_id)
rownames(dat) <- ids$symbol
  #再踩坑一次蚁滋！這里rownames(dat) <- ids$symbol #這個時(shí)候由于rownames的symbol有重復(fù)的，會報(bào)錯赘淮，所以先需要進(jìn)行去重復(fù),還是和上面相同的報(bào)錯枢赔，再一次提醒我們：不允許有重復(fù)的'row.names'，因此需要基因名去重復(fù)


##################基因名symbol去重復(fù)如下：
if(T){
  ids$median=apply(dat,1,median) #ids新建median這一列拥知，列名為median踏拜，同時(shí)對dat這個矩陣按行操作，取每一行的中位數(shù)低剔，將結(jié)果給到median這一列的每一行
  ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序速梗，將對應(yīng)的行賦值為一個新的ids
  ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項(xiàng)，'!'為否襟齿，即取出不重復(fù)的項(xiàng)姻锁，去除重復(fù)的gene ，保留每個基因最大表達(dá)量結(jié)果s
  dim(ids)
  ids[1:3,1:3]
  dat=dat[ids$probe_id,] #新的ids取出probe_id這一列猜欺，將dat按照取出的這一列中probe_id所在的每一行組成一個新的dat位隶，這時(shí)獲得的dat的探針名就是和ids的probe_id這一列是完全對應(yīng)的
  rownames(dat)=ids$symbol#把ids的symbol這一列中的每一行給dat，作為dat的rowname
  dim(dat)
  dat[1:4,1:4]
  #identical(rownames(dat),ids$symbol) #可驗(yàn)證开皿，也可不驗(yàn)證涧黄，只要能確保 dat=dat[ids$probe_id,]這步是正確的
  boxplot(dat)
  dat=log(dat+1)
  save(dat,ids,group_list,file = 'step1.4_dat_ids_group_list.Rdata')

}