芯片分析中常常會遇到Affymetrix Human Transcriptome Array 2.0芯片,因爲目前尚未現成的R包能夠用，所以分析方法也不統一。見生信技能樹Jimmy老師HTA2.0芯片比較麻煩，其實這類常見的有3個平臺，3種類型：web
GPL17586 [HTA-2_0] Affymetrix Human Transcriptome Array 2.0 [transcript (gene) version]數據庫GPL19251 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array [probe set (exon) version]apiGPL16686 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array [transcript (gene) version]微信
對於這三種平臺能夠去Affymetrix的官網去查看其區別，也能夠去NCBI去查看。編輯器
1、得到芯片平臺信息文件
一般基因芯片分析，通常要下載其平臺信息。通常來講咱們下載GPL是爲了獲得芯片的探針對應基因ID的關係列表。詳情能夠了解：解讀GEO數據存放規律及下載，一文就夠一文就夠-從GEO數據庫下載獲得表達矩陣首先是GPL17586平臺的芯片，下載其對應的平臺文件GPL17586.soft.gz，這類文件一般都比較大，加上國內下載速度太慢，一般都是等了好久，仍是下載不了。學習
查看GPL17586平臺下單個的GSE對應的GSEXXX.soft.gz文件，發現其信息與GPL17586.soft.gz相同；下載單個GSE對應的soft文件後，一樣能夠作id轉換。flex
下面以GPL17586平臺下的GSE110359爲例，進行id轉換，首先是下載GSE110359對應的GSE110359_family.soft.gz文件this
一、讀入下載好的soft文件
##
rm(list = ls())
options(stringsAsFactors = F)
#加載R包
library(GEOquery)
gse <- getGEO(filename = "GSE110359_family.soft.gz",destdir = ".")
str(gse)
length(gse)
二、提取探針、探針對應的基因及其位於染色體上的位置等信息
id_probe <- gse@gpls$GPL17586@dataTable@table
dim(id_probe)
head(id_probe)
id_probe[1:4,1:15]
View(head(id_probe))## you need to check this , which column do you need
提取須要的第1列（ID）或者第2列（probeset_id）和第8列（gene_assignment）。固然也能夠不提取，每一列都保留。
probe2gene <- id_probe[,c(2,8)]
三、提取第8列 gene_assignment中的基因名稱,並添加到probe2gene
library(stringr)
probe2gene$symbol=trimws(str_split(probe2gene$gene_assignment,'//',simplify = T)[,2])
plot(table(table(probe2gene$symbol)),xlim=c(1,50))
head(probe2gene)
dim(probe2gene)
View(head(probe2gene))
ids2 <- probe2gene[,c(1,3)]
View(head(ids2))
ids2[1:20,1:2]#含有缺失值
table(table(unique(ids2$symbol)))#30907 ,30906個基因，一個空字符
save(ids2,probe2gene,file='gse-probe2gene.Rdata')
以上爲id轉換的第一種方法，下面看第二種方法：url
四、使用biomaRt包進行id轉換
biomaRt包其實也依賴於網速spa
load("gse-probe2gene.Rdata")
dim(probe2gene)
View(head(probe2gene))
# 加載biomaRt包
library(biomaRt)
value <- probe2gene$probeset_id
attr <- c("affy_hta_2_0","hgnc_symbol")
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
ids <- getBM(attributes = attr,
filters = "affy_hta_2_0",
values = value,
mart = ensembl,
useCache = F)
dim(ids)#[1] 1041 2
View(head(ids))
save(ids,file = "GPL17586_ids.Rdata")
#去重以後
attributes <- listAttributes(ensembl)
View(attributes) # 查看支持的芯片轉換格式
save(ids,ensembl,y,file = "ensembl.Rdata")
plot(table(table(ids$hgnc_symbol)),xlim=c(1,50))
table(table(unique(ids$hgnc_symbol)))#去重以後有29262，丟失了一不少
# 去重複
ids3 <- ids[!duplicated(ids$hgnc_symbol),]
綜上，能夠看出兩種方法獲得的基因數量差異不大都是從7萬多個探針中，得到了差很少3萬個基因。
五、表達矩陣和基因id的合併
下面就是表達矩陣和id的合併了；下載表達矩陣，推薦使用Jimmy的萬能包GEOmirror一行命令搞定。
library(GEOmirror)
geoChina(gse = "GSE110359", mirror = "tercent")
library(GEOmirror)
gset <- geoChina(gse = "GSE110359", mirror = "tercent")
gset
a=exprs(gset[[1]])
a[1:4,1:4]
gset[[1]]@annotation
#過濾表達矩陣
exprSet <- a
library(dplyr)
exprSet <- exprSet[rownames(exprSet) %in% ids2$probeset_id,]
dim(exprSet)
exprSet[1:5,1:5]
#ids過濾探針
ids <- ids2[match(rownames(exprSet),ids2$probeset_id),]
dim(ids)
ids[1:5,1:2]
#ids2[1:5,1:2]
#合併表達矩陣和ids
idcombine <- function(exprSet, ids){
tmp <- by(exprSet,
ids$symbol,
function(x) rownames(x)[which.max(rowMeans(x))])
probes <- as.character(tmp)
print(dim(exprSet))
exprSet <- exprSet[rownames(exprSet) %in% probes,]
print(dim(exprSet))
rownames(exprSet) <- ids[match(rownames(exprSet), ids$probeset_id),2]
return(exprSet)
}
new_exprSet <- idcombine(exprSet,ids)
new_exprSet[1:4,1:6]
dim(new_exprSet)
rownames(new_exprSet)
save(new_exprSet,file = "new_exprSet.Rdata")
接下來就是質控、差別分析，火山圖、GO和KEEG分析，生信技能樹上已經有不少這類的推文了，這裏就不作演示了。
該方法也適用於GPL16686與GPL19251平臺的芯片。只是GPL16686的信息這樣的能夠用Y叔叔的包進行轉換id。
GPL16686平臺信息
GPL19251平臺信息
轉換id以後總會有不少探針得不到對應的基因或者不少探針對應一個基因。其實，基因和探針的關係是多對多的關係。
友情推薦:「ID轉換和生存分析」羣的釘釘羣號：35371384，對這幾個細節知識點感興趣的能夠加入，咱們這個月18號(下週六晚上)八點準時授課。ID轉換靠的是深厚的背景知識加上一點代碼技巧
寫在後面
由於這個投稿是兩個星期前的，可是生信技能樹學習筆記太多了，因此排班到了今天，上面的ID轉換公開課已經結束了，可是今晚有生存分析公開課，感興趣仍然是能夠進去哈，見：生存分析的10年和20年時間點
本文分享自微信公衆號 - 生信技能樹（biotrainee）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。