cMAP在线分析——旧版build2中Tag list的制作？

之前分享过cMap分析专题：

• Connectivity Map(cMap)的探索应用（一）

• Connectivity Map(cMap)的探索应用（二）

• cMAP在线分析——旧版build2的使用

• cMAP新版clue的使用——List Marker

• cMAP新版clue的使用——query

• cMAP新版clue的使用——Touchstone

• cMAP新版clue的使用——Data Library

• cMap新版clue的使用——Morpheus

• cMAP新版clue的使用——CELL APP

• Connectivity Map(cMap)的探索应用（三）

然后，有部分小伙伴咨询在“cMAP在线分析——旧版build2的使用”这篇教程中tag list的制作问题，今天就这部分做简要介绍！

首先来看下官方说明：

tag list: A list of probe sets from the feature set representing either the 'up' genes in a signature (the 'up tag list') or the 'down' genes in a signature (the 'down tag list'). One up tag list and one down tag list together represent a signature and both are required to execute a query using the quick query page or to upload a signature using the load signature page. Note that the total number of probe sets in the up tag list and down tag lists may not together exceed 1000. Text files with the .grp extension are used to contain tag lists. A .grp file has one probe set on each line.
Tag lists may be conveniently constructed from a signature populated with UniGene identifiers, RefSeq identifiers, Genbank identifiers, gene symbols, etc using the 'Batch Query' tool on NetAffx and selecting "HG-U133A" from the 'Select a GeneChip Array' menu. A variety of tools are available for the generation of tag lists from signatures from non-human species. These include GeneCruiser (http://www.genecruiser.org). Our preferred method uses UniGene and HomoloGene from NBCI (http://www.ncbi.nlm.nih.gov/) to obtain human UniGene identifiers for signature genes for submission to NetAffx, as described above.

关键大意：tag list通常分上调基因list和下调基因list，两者的基因总和不能大于1000个，同时要以.grp格式保存。基因名称可以是 UniGene ID, RefSeq ID, Genbank ID, gene symbols等，然后通过NetAffx的Batch Query功能进行转换，得到HG-U133A芯片中基因对应的探针ID。如何是其他物种，需要先通过同源性比较转化为人的基因，然后再提交给NetAffx完成探针ID的转化。至于如何进行同源性比较，官网提及 GeneCruiser（这个网址没打开）和NCBI HomoloGene等方式。

probe set: The collection of match and mismatch oligonucleotides on an Affymetrix GeneChip microarray designed against a given transcript which together allow the relative level of that transcript to be estimated. Probe sets are uniquely identified with a code number that, by convention, ends with "_at" ( eg 200800_s_at). Tag lists and instances are populated with probe sets from the feature set. Detailed descriptions of individual probe sets can be found through NetAffx.
feature set: The universe of probe sets recognized by cmap, defined as the 22,283 probe sets on the HG-U133A array.
HG-U133A: Affymetrix GeneChip Human Genome U133A Array (part number 510681). The probe sets on this array define the feature set.

关键大意：taglist的探针ID的来源主要是HG-U133A这款芯片，共计22283条探针信息，每个探针ID都有独一无二的数字编码，并以_at结尾为特征。

从GEO Database中检索GPL96，即可获取该芯片的注释信息（如下图），这时候可以通过对应关系，将人的基因名字（可以是官方gene symbols/RefSeq ID ）转化为芯片的探针ID。

image

目前待解决的问题：

1. 常规情况下，如何将基因转化为HG-U133A芯片的探针ID?

方式一：根据官方说明，拿到人的基因名字（UniGene ID/RefSeq ID/ Genbank ID/ gene symbols）后，这里采用的是 gene symbols，将基因名字提交给NetAffx（http://www.affymetrix.com/analysis/netaffx/index.affx）的Batch Query，完成基因名字—探针ID的转换（如下图）。

image

方式二：本地下载好HG-U133A芯片的注释文件（即上文中的GPL96平台信息），然后通过R语言直接将提供的基因名称与注释文件匹配，得到基因对应的探针ID。

rm(list = ls())
GPL96=read.table("GPL96-57554.txt",  header = TRUE,fill = T,sep = "\t",
                 comment.char = "#",
                 stringsAsFactors = FALSE,
                 quote = "")
colnames(GPL96)#1，11列为探针ID和gene symbol
data=read.table("基因名称.txt")
colnames(data)<- "Gene Symbol"
##获取基因名称和探针ID之间的关系
gene_tran_ID <- GPL96[GPL96$Gene.Symbol %in% data$`Gene Symbol`,c(1,11)]
# ID Gene.Symbol
# 1087  201559_s_at       CLIC4
# 1088    201560_at       CLIC4
# 1706    202178_at       PRKCZ
# 3406    203879_at      PIK3CD
# 3702    204175_at      ZNF593
# 5154    205627_at         CDA
# 8728    209234_at       KIF1B
# 10668 211230_s_at      PIK3CD
# 11733 212348_s_at       KDM1A
# 12649   213268_at      CAMTA1
# 17927 218562_s_at      TMEM57
# 18495   219131_at      UBIAD1
# 19228 219864_s_at       RCAN3
# 21150 221790_s_at     LDLRAP1
# 21241 221881_s_at       CLIC4
# 22116    57082_at     LDLRAP1

2. 如何将其它物种的基因转化为HG-U133A芯片的探针ID?

关键步骤在于将其他物种基因通过同源性映射到人的基因，然后再通过NetAffx的Batch Query，完成基因名字—探针ID的转换。那如何进行不同物种之间的同源基因转化呢？

• 在线同源基因转化：http://refdic.rcai.riken.jp/tools/matchom.cgi。如下图，提交了一份小鼠的基因名称（gene symbol）,想将其转化为同源的人基因名称，最后结果：提交17个小鼠的基因，成功匹配到14个人的同源基因。

image

• 本地转化—homologene包，该包基于NCBI的homologene数据库开发，旨在用于匹配不同物种之间的同源基因。该包的中包含多个函数，主要使用的是homologene，它的使用方式如下：
  Usage：homologene(genes, inTax, outTax, db = homologene::homologeneData)
  Arguments（参数）
• genes：A vector of gene symbols or NCBI ids
• inTax：taxid of the species that the input genes are coming from
• outTax：taxid of the species that you are seeking homology
• db：Homologene database to use.
  Examples：homologene(c('Eno2','17441'), inTax = 10090, outTax = 9606)

从上述使用方式来看，只需要准备gene symbols or NCBI ids，然后说清楚，由什么物种（inTax）转化到什么物种（outTax）即可。这里用NCBI的taxid来表示不同物种。NCBI人为的给自然界所有物种都给了一个编号，这个编号就是taxid，比如其中人类的编号是 9606，10090是小鼠。

10090   Mus musculus
10116   Rattus norvegicus
28985   Kluyveromyces lactis
318829  Magnaporthe oryzae
33169   Eremothecium gossypii
3702    Arabidopsis thaliana
4530    Oryza sativa
4896    Schizosaccharomyces pombe
4932    Saccharomyces cerevisiae
5141    Neurospora crassa
6239    Caenorhabditis elegans
7165    Anopheles gambiae
7227    Drosophila melanogaster
7955    Danio rerio
8364    Xenopus (Silurana) tropicalis
9031    Gallus gallus
9544    Macaca mulatta
9598    Pan troglodytes
9606    Homo sapiens
9615    Canis lupus familiaris
9913    Bos taurus

通过homologene()函数可以非常快速的进行不同物种间同源基因的转化，且比在线工具的匹配度更高。

rm(list = ls())
options(download.file.method = 'libcurl')
options(url.method='libcurl')
install.packages("homologene")
library("homologene")
data=read.table("小鼠基因名称.txt")
dat=as.character(t(data))
genelist <- dat
# [1] "Prkcz"   "Camta1"  "Pik3cd"  "Kif1b"   "Ubiad1"  "Fbxo6"   "Agtrap"  "Rsc1a1"  "Cda"     "Kdm1a"   "Pithd1"  "Rcan3"   "Clic4"  
# [14] "Tmem57"  "Ldlrap1" "Pdik1l"  "Zfp593" 
output <- homologene(genelist, inTax = 10090, outTax = 9606)
# 10090    9606 10090_ID 9606_ID
# 1    Prkcz   PRKCZ    18762    5590
# 2   Camta1  CAMTA1   100072   23261
# 3   Pik3cd  PIK3CD    18707    5293
# 4    Kif1b   KIF1B    16561   23095
# 5   Ubiad1  UBIAD1    71707   29914
# 6    Fbxo6   FBXO6    50762   26270
# 7   Agtrap  AGTRAP    11610   57085
# 8   Rsc1a1  RSC1A1    69994    6248
# 9      Cda     CDA    72269     978
# 10   Kdm1a   KDM1A    99982   23028
# 11  Pithd1  PITHD1    66193   57095
# 12   Rcan3   RCAN3    53902   11123
# 13   Clic4   CLIC4    29876   25932
# 14  Tmem57  TMEM57    66146   55219
# 15 Ldlrap1 LDLRAP1   100017   26119
# 16  Pdik1l  PDIK1L   230809  149420
# 17  Zfp593  ZNF593    68040   51042

遇到的问题？

问：有500个gene（都是gene名），在转换探针ID后只剩下不到一半了，匹配率很低是为什么？
答：可能的原因：

1. 这500个基因可能是其它物种的，如果是其它物种需要转化成人的同源基因，再行转化探针ID。
2. 基因名称有很多别名，可能正好你提供的名字和平台注释提供的名字不同而被过滤。但实际操作过程中发现，这个影响因素很小。撇开匹配度的问题，本身HG-U133A芯片收录的基因本身就是不全面的，总共收录的基因也就两万多个。现在很多基因在HG-U133A芯片中就没有涉及相应探针。
3. 若是由芯片本身基因设计不足带来的低匹配度，这该怎么办？个人觉得，技术的革新是正常的，基因检测数量的增加也很正常。但cmap分析基于的小分子表达数据库的数据大多是基于HG-U133A这款芯片，是将这个数据库作为比较参照的，所以这就导致我们要做这个分析，那么就要保证咱们的数据最好也是HG-U133A芯片数据，但显然不可能。因此，建议咱们的差异数据结果出来以后，先与HG-U133A芯片中收录的基因取交集，在交集的差异基因基础上，在分别取上调500基因，下调500基因去进行cmap分析。这样可以有效避免由于技术差异带来的低匹配基因，从而避免由于技术革新对cmap分析结果带来的影响。

往期回顾
Connectivity Map(cMap)的探索应用（一）
Connectivity Map(cMap)的探索应用（二）
cMAP在线分析——旧版build2的使用
Connectivity Map(cMap)的探索应用（三）
cMAP新版clue的使用——List Marker
cMAP新版clue的使用——query
cMAP新版clue的使用——Touchstone
cMAP新版clue的使用——Data Library
cMap新版clue的使用——Morpheus
cMAP新版clue的使用——CELL APP

今天的内容就到这里，更多内容可关注公共号“YJY技能修炼”~~~