在挖掘GEO数据库的过程中，第一步就是获取并整理数据集的表达矩阵。有部分文章会好心的在GSE数据集页面上传整理好的表达矩阵，但绝大多数文章都仅仅是上传基础的信息，需要自行整理表达矩阵。

（如果原作者上传了整理后的表达矩阵，一般都会放在GSE数据集页面的Supplementaryfile处，但是此处的文件也可能是作者的其他附件，需要自行下载并检查是不是我们需要的表达矩阵）

那么整理表达矩阵的关键一步，就是将探针名（probe）转换为基因名（genesymbol），即注释探针。今天小编会介绍第二种注释探针的方法：使用GEO数据集官方推荐的GEOquery包下载GPL文件进行注释。

一、数据集准备

这次选择的数据集是GSE，脑胶质瘤患者的mRNA表达芯片数据。芯片平台为GPL，[HG-UPlus2]AffymetrixHumanGenomeUPlus2.0Array。GPL是非常经典的芯片，在GEO数据库中就有个数据集，共16万样本使用过该芯片分析表达数据。

二、实际操作

本次操作也很简单，只需使用GEOquery包下载好表达数据和GPL文件即可

library(GEOquery)library(Biobase)gset-getGEO("GSE",GSEMatrix=TRUE,AnnotGPL=FALSE)if(length(gset)1)idx-grep("GPL",attr(gset,"names"))elseidx-1gset-gset[[1]]exprSet-exprs(gset);dim(exprSet)

到这里我们已经成功下载表达矩阵，接下来就是下载GPL文件

library(GEOquery)gpl-getGEO(GPL)GPL=Table(gpl)colnames(GPL)#[1]"ID""GB_ACC"#[3]"SPOT_ID""SpeciesScientificName"#[5]"AnnotationDate""SequenceType"#[7]"SequenceSource""TargetDescription"#[9]"RepresentativePublicID""GeneTitle"#[11]"GeneSymbol""ENTREZ_GENE_ID"#[13]"RefSeqTranscriptID""GeneOntologyBiologicalProcess"#[15]"GeneOntologyCellularComponent""GeneOntologyMolecularFunction"

可以看到第十一列是我们所需要的内容"GeneSymbol"，点进去看看

但是我们发现，有的"GeneSymbol"里面，有两个基因SYMBOL，因此我们需要拆分开后再继续注释

library("plyr")GPL-GPL[GPL$`GeneSymbol`!=,]dim(GPL)split_gene-strsplit(as.character(GPL$`GeneSymbol`),"///")probe2gene-mapply(cbind,GPL[,1],split_gene)ID2gene-as.data.frame(do.call(rbind,probe2gene))head(ID2gene)exprSet-exprSet[ID2gene$V1,]dim(exprSet)exprSet[1:5,1:6]tail(sort(table(ID2gene[,1])),n=12L)

##有一些基因对应多个探针或者一个探针对应多个基因，所以我们这里只保留表达量最大的探针

ID2gene$max-apply(exprSet,1,max)ID2gene-ID2gene[order(ID2gene$V2,ID2gene$max,decreasing=T),]ID2gene-ID2gene[!duplicated(ID2gene$V2),]ID2gene-ID2gene[!duplicated(ID2gene$V1),]head(ID2gene)anno=ID2gene

##最后只需要把行名调成一致即可

tmp=rownames(exprSet)%in%anno[,1]exprSet=exprSet[tmp,]dim(exprSet)match(rownames(exprSet),anno$V1)anno=anno[match(rownames(exprSet),anno$V1),]match(rownames(exprSet),anno$V1)dim(exprSet)dim(anno)rownames(exprSet)=anno$V2dim(exprSet)exprSet[1:5,1:5]#GSMGSMGSMGSMGSM#MIR.....9#HSPA.....6#TTLL.....1#ADAM.....4#PRR.....2

最后再点开看看是否注释成功

实际上，如果最初的函数getGEO下载数据比较完整的话，可能gset文件内就有我们的注释文件，比如这个数据集的gset我们进一步探索就会发现gset里面还有更多内容。

根据英文名称来看，assayData是表达矩阵，phenoData是该数据集的表型文件，exprimentData和protocolData可能是该数据集进行实验的时候的相关实验条件。annotation和featureData内就可能存在注释信息。查看后发现annotation写的是芯片内容。

gset

annotation[1]"GPL"

进一步探索发现，注释信息就藏在featureData内

featureData-gset

featureData

dataView(featureData)

对比后发现从featureData内获得的注释信息，和使用GEOquery包下载的注释信息完全一致。因此以后寻找注释信息时，可以通过这两种途径去获得GPL文件~

总的来说，上一期和这一期都是利用GPL文件进行注释，那么还有没有更方便的方法进行探针注释呢~~请期待下一期内容!

注：此推文未经许可禁止转载！阅读推荐：

生信实操丨一个生信菜鸟的上道经验分享-转录组测序（富集分析绘图篇）

生信实操丨高级转录组分析WGCNA应该这么做

生信实操丨带你复现单细胞转录组纯分析文章（一）

生信实操

转录本编码区预测神器-TransDecoder详细教程

生信实操

注释芯片探针的三种方法（一）

如果您或者科室有科研上的困扰扫码备注：科研思路探讨预览时标签不可点收录于话题#个上一篇下一篇


转载请注明地址:http://www.tanhuaa.com/gyth/7889.html