芯片探针到基因组区段坐标的映射

最近接到粉丝求助，有一篇文献写到：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 arm1 chr1 0 2300000 36.33 gneg p2 chr1 2300000 5300000 36.32 gpos25 p3 chr1 5300000 7100000 36.31 gneg p4 chr1 7100000 9200000 36.23 gpos25 p5 chr1 9200000 12600000 36.22 gneg p6 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 Chromosome1 1053_at -74231501 72 1053_at -74231501 73 117_at 161524539 14 121_at -113215996 25 1255_g_at 42155376 66 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_probesgr_cytobands= GRanges( seqnames = hc$chr, ranges = IRanges( hc$start, hc$end))gr_cytobandsfindOverlaps(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 arm3 117_at 161524539 1 chr1 158800000 163800000 23.3 gneg q5 1255_g_at 42155376 6 chr6 40600000 45200000 21.1 gneg p6 1316_at 40062964 17 chr17 37800000 41900000 21.31 gneg q7 1316_at 40062192 17 chr17 37800000 41900000 21.31 gneg q8 1316_at 40062964 17 chr17 37800000 41900000 21.31 gneg q12 1431_at 133527362 10 chr10 130500000 135374737 26.3 gneg q


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