最近接到粉丝求助，有一篇文献写到：We found that 16 differentially expressed genes (Table 2) represented by specific probe sets (‘_at’ suffix) mapped to previously reported linkage peaks on chromosomes 1p34, 5q12, 9q22, 9q34, 13q32, 14q32, and 20q13.

也就是说，差异基因这个时候是由affymetrix的芯片探针表示，然后作者注释到了基因组区段，也就是 1p34, 5q12这样的标识，行话叫做：cytoBand。

ChromHeatMap 包存放有 cytoBand坐标信息

早在教程：染色体全局可视化 ，我就提到过ChromHeatMap 包 存放有 cytoBand坐标信息，查看的代码如下：

# BiocManager::install('ChromHeatMap')
library('ChromHeatMap')
data("cytobands")
hc=cytobands[[1]]
head(hc)

可以看到是如下所示的数据框：

> head(hc)
 chr start end band stain arm
1 chr1 0 2300000 36.33 gneg p
2 chr1 2300000 5300000 36.32 gpos25 p
3 chr1 5300000 7100000 36.31 gneg p
4 chr1 7100000 9200000 36.23 gpos25 p
5 chr1 9200000 12600000 36.22 gneg p
6 chr1 12600000 16100000 36.21 gpos50 p

其6个字段分别是：染色体编号(chr)、在染色体中的起始位置(start)、终止位置(end)、cM (band)、染色标识(stain)，长臂短臂(arm， short (p) and long (q) )。

探针的坐标在各个芯片包也可以获得

比如 hgu133plus2.db 芯片包，如下：

library(hgu133plus2.db)
probe2pos=toTable(hgu133plus2CHRLOC)
head(probe2pos)

坐标如下：

> head(probe2pos)
 probe_id start_location Chromosome
1 1053_at -74231501 7
2 1053_at -74231501 7
3 117_at 161524539 1
4 121_at -113215996 2
5 1255_g_at 42155376 6
6 1316_at 40062964 17

两个坐标在R里面使用grange进行overlap

初学者需要自己去读文档，了解 https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html 包的方方面面，才能把它用好，而且用的妙。

需要首先把它们两个坐标转为 GRanges 对象，然后 findOverlaps 函数即可

library(GenomicRanges)
gr_probes= GRanges(
 seqnames = paste0('chr',probe2pos$Chromosome),
 ranges = IRanges( probe2pos$start_location, probe2pos$start_location+1),
 strand = ifelse(probe2pos$start_location>1,'+','-'),
 probe=probe2gene$probe_id
 )
gr_probes
gr_cytobands= GRanges(
 seqnames = hc$chr,
 ranges = IRanges( hc$start, hc$end)
)
gr_cytobands
findOverlaps(gr_probes,gr_cytobands)

可以看到，两个坐标系统就对应起来了。

> findOverlaps(gr_probes,gr_cytobands)
Hits object with 39409 hits and 0 metadata columns:
 queryHits subjectHits
 <integer> <integer>
 [1] 3 43
 [2] 5 651
 [3] 6 319
 [4] 7 319
 [5] 8 319
 ... ... ...
 [39405] 85850 144
 [39406] 85851 144
 [39407] 85852 144
 [39408] 85853 144
 [39409] 85854 144
 -------
 queryLength: 85862 / subjectLength: 862
>

最后根据坐标检索即可，代码如下：

tmp=findOverlaps(gr_probes,gr_cytobands)
comb=cbind(probe2gene[tmp@from,],hc[tmp@to,])
head(comb)

结果如下：

> head(comb)
 probe_id start_location Chromosome chr start end band stain arm
3 117_at 161524539 1 chr1 158800000 163800000 23.3 gneg q
5 1255_g_at 42155376 6 chr6 40600000 45200000 21.1 gneg p
6 1316_at 40062964 17 chr17 37800000 41900000 21.31 gneg q
7 1316_at 40062192 17 chr17 37800000 41900000 21.31 gneg q
8 1316_at 40062964 17 chr17 37800000 41900000 21.31 gneg q
12 1431_at 133527362 10 chr10 130500000 135374737 26.3 gneg q

但是可能有一个问题，我没有去深究这两个包里面的坐标信息的参考基因组版本问题。大家如果要实战，需要额外注意，我这里仅仅是教程哈。

文末友情推荐

要想真正入门生物信息学建议务必购买全套书籍，一点一滴攻克计算机基础知识，书单在：什么，生信入门全套书籍仅需160 。
如果大家没有时间自行慢慢摸索着学习，可以考虑我们生信技能树官方举办的学习班：

• 数据挖掘学习班第5期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。
• 生信爆款入门-第7期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。

如果你课题涉及到转录组，欢迎添加一对一客服：详见：你还在花三五万做一个单细胞转录组吗？

号外：生信技能树知识整理实习生招募，长期招募，也可以简单参与软件测评笔记撰写，开启你的分享人生！另外：绝大部分生信技能树粉丝都没有机会加我微信，已经多次满了5000好友，所以我开通了一个微信好友，前100名添加我，仅需150元即可，3折优惠期机会不容错过哈。我的微信小号二维码在：0元，10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》