基因组的 mappability 是指测序产生的短读长能够被唯一且准确比对到参考基因组特定区域的能力。它反映了基因组不同区域因序列特征（如重复序列、低复杂度区域等）对比对结果可靠性的影响。例如，高度重复区域（如 Alu 元件）的 mappability 通常较低，因为 reads 可能匹配到多个位置。

在实际分析中，有时我们会过滤低 mappability 区域的突变，避免因比对错误导致的假阳性误判。本文将介绍如何使用生物信息学工具 GenMap 获得基因组的 mappability 信息。

GenMap

GenMap github page：

https://github.com/cpockrandt/genmap

GenMap paper：

GenMap: ultra-fast computation of genome mappability

https://doi.org/10.1093/bioinformatics/btaa222

相较于先前的 mappability 计算工具（例如 GEM），GenMap 的搜索策略有所优化并且算法不依赖于启发式方法，因此具有更高的运行效率和准确性。此外 GenMap 支持多种输出格式，便于下游灵活分析。

安装

推荐通过 conda 安装 GenMap ：

shell
1
conda install -c bioconda genmap

也可以通过下载源码进行编译：

shell
1
2
3
4
git clone --recursive https://github.com/cpockrandt/genmap.git
mkdir genmap-build && cd genmap-build
cmake ../genmap -DCMAKE_BUILD_TYPE=Release
make genmap

运行以下命令以检查安装是否完成：

shell
1
genmap --help

运行

构建索引

在计算 mappability 前，需要对目标基因组文件建立索引：

shell
1
genmap index -F /path/to/fasta.fasta -I /path/to/index/folder

-F 指定参考基因组序列路径，-I 指定索引输出目录位置，请注意此处需要的内存较大：

It needs about 6n space in main memory (or 10n for fasta files >2GB). n is the number of bases in your fasta file(s). It might be more or less depending on the number and length of the individual sequences.

如果你分析的物种为常见的模式生物，可以在以下表格下载 GenMap 提前构建好的索引文件：

Genome Index size (compressed) Download
Human GRCh38 [1] 5.4 GB GRCh38 index
Human hs37-1kg [2] 5.4 GB hs37-1kg index
Mouse GRCm38 4.9 GB GRCm38 index
D. melanogaster dm6 0.2 GB dm6 index
C. elegans ce11 0.1 GB ce11 index
Wheat T. aestivum ta45 [3] 21.9 GB ta45 index
Source: https://github.com/cpockrandt/genmap#pre-built-indices

如果你拟分析的物种不在该列表中且内存量不足，可以考虑调整 genmap index 的 -S 参数至 20，不过这会导致后续计算 mappability 的速度减慢。

计算 mappability

在构建好索引后，可以运行以下命令计算基因组的 mappability：

shell
1
genmap map -K 100 -E 2 -I /path/to/index/folder -O /path/to/output/folder -w

各参数的含义如下：

  • -K：使用的 k-mer 长度，此处应当依据自己使用的测序数据读长决定，例如如果你使用的是 150bp 的 reads 数据，那么该参数应当调整为 150 以对应实际比对时的情况。
  • -E：允许的错配数量，此处错配数量的选择需根据应用场景决定，例如 CRISPR 设计需要特定的错配容忍度、在变异调用中映射 reads 时允许一定的碱基错配数量等。
  • -I：索引文件所处路径 。
  • -O：结果文件输出路径，可包含文件前缀，例如 ./hg38.mappability.k30_e2（此时输出路径为当前所在目录）。
  • -w：输出 wig 文件格式，此外可选的还有 -t / -bg / -d 等参数，对应的文件具体信息可见 genmap map --help。

在输出 wig 文件后，可通过 wigToBigWig 将 wig 文件转为 bigwig 文件（可选）：

shell
1
2
3
4
# 安装
conda install bioconda::ucsc-wigtobigwig
# 使用
wigToBigWig [input wig] [chrom size file] [output bigwig]

各参数的含义如下：

  • [input wig]：输入的 wig 文件。
  • [chrom size file]：染色体大小信息文件，可以通过 samtools faidx 获得。
  • [output bigwig]：输出的 bigwig 文件。

拓展

需要注意的是，GenMap 是为了短 reads 和最多 4 个错配开发的，并且目前并不支持 InDel（插入/缺失）的计算。

此外，你可以在 genmap map 的运行中通过指定 -S 参数，以计算特定区域内的 mappability：

shell
1
genmap map -S exome.bed -K 100 -E 2 -I /path/to/index/folder -O /path/to/output/folder -w

以上述命令为例，你可以将 exome.bed 替换为自己的外显子组坐标文件，以使 genmap 仅对外显子组中的区域进行计算。

参考材料

  1. GenMap: ultra-fast computation of genome mappability

  2. https://github.com/cpockrandt/genmap/issues/19

  3. https://github.com/cpockrandt/genmap/issues/5