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

时间:2022-07-22
本文章向大家介绍芯片探针到基因组区段坐标的映射，主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项，具有一定的参考价值，需要的朋友可以参考一下。

最近接到粉丝求助，有一篇文献写到：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

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

随机文章
本站知识点必读
最近更新