seurat提取表达矩阵_GPL17586、GPL19251和GPL16686平
台芯⽚ID转换
芯⽚分析中经常会遇到Affymetrix Human Transcriptome Array 2.0芯⽚，由于⽬前还没有现成的R包可以⽤，因此分析⽅法也不统⼀。见⽣信技能树Jimmy⽼师HTA2.0芯⽚⽐较⿇烦，其实这类常见的有3个平台，3种类型：
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]
GPL16686 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array [transcript (gene) version]
对于这三种平台可以去Affymetrix的官⽹去查看其区别，也可以去NCBI去查看。
⼀、获得芯⽚平台信息⽂件
通常基因芯⽚分析，⼀般要下载其平台信息。⼀般来说我们下载GPL是为了得到芯⽚的探针对应基因ID的关系列表。详情可以了解：解读GEO数据存放规律及下载，⼀⽂就够⼀⽂就够-从GEO数据库下载得 到表达矩阵⾸先是GPL17586平台的芯⽚，下载其对应的平台⽂件，这类⽂件通常都⽐较⼤，加上国内下载速度太慢，通常都是等了很久，还是下载不了。
查看GPL17586平台下单个的GSE对应的⽂件，发现其信息与相同；下载单个GSE对应的soft⽂件后，同样可以做id转换。
下⾯以GPL17586平台下的GSE110359为例，进⾏id转换，⾸先是下载GSE110359对应的GSE110359_⽂件
1、读⼊下载好的soft⽂件
##
rm(list = ls())
options(stringsAsFactors = F)
#加载R包
library(GEOquery)
gse
str(gse)
length(gse)
2、提取探针、探针对应的基因及其位于染⾊体上的位置等信息
id_probe
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
GPL17586平台信息
提取需要的第1列(ID)或者第2列(probeset_id)和第8列(gene_assignment)。当然也可以不提取，每⼀列都保留。
probe2gene
3、提取第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
View(head(ids2))
ids2[1:20,1:2]#含有缺失值
table(table(unique(ids2$symbol)))#30907 ,30906个基因，⼀个空字符
save(ids2,probe2gene,file='gse-probe2gene.Rdata')
以上为id转换的第⼀种⽅法，下⾯看第⼆种⽅法：
4、使⽤biomaRt包进⾏id转换
biomaRt包其实也依赖于⽹速
load("gse-probe2gene.Rdata")
dim(probe2gene)
View(head(probe2gene))
# 加载biomaRt包
library(biomaRt)
value
attr
ensembl
ids
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
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
综上，可以看出两种⽅法得到的基因数量差别不⼤都是从7万多个探针中，获得了差不多3万个基因。
5、表达矩阵和基因id的合并
下⾯就是表达矩阵和id的合并了；下载表达矩阵，推荐使⽤Jimmy的万能包GEOmirror⼀⾏命令搞定。library(GEOmirror)
geoChina(gse = "GSE110359", mirror = "tercent")
library(GEOmirror)
gset
gset
a=exprs(gset[[1]])
a[1:4,1:4]
gset[[1]]@annotation
#过滤表达矩阵
exprSet
library(dplyr)
exprSet
dim(exprSet)
exprSet[1:5,1:5]
#ids过滤探针
ids
dim(ids)
ids[1:5,1:2]
#ids2[1:5,1:2]
#合并表达矩阵和ids
idcombine
tmp
ids$symbol,
function(x) rownames(x)[which.max(rowMeans(x))])
probes
print(dim(exprSet))
exprSet
print(dim(exprSet))
rownames(exprSet)
return(exprSet)
}
new_exprSet
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平台信息
GPL19251平台信息
转换id之后总会有很多探针得不到对应的基因或者很多探针对应⼀个基因。其实，基因和探针的关系是多对多的关系。
友情推荐:“ID转换和⽣存分析”的钉钉号：35371384，对这⼏个细节知识点感兴趣的可以加⼊，我们这个⽉18号(下周六晚上)⼋点准时授课。ID转换靠的是深厚的背景知识加上⼀点代码技巧
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议