crossmap转基因组版本

crossmap 软件转基因组版本。

安装

可以用 conda 或 pip 安装

conda install -c bioconda crossmap

从v0.7.0版本开始，主程序 CrossMap.py 已更名为 CrossMap，原因是pyproject.toml下脚本名称使用“.”存在限制。为保持向后兼容性，请在您的~/.bashrc文件中添加以下行：

alias CrossMap.py='CrossMap'

输入文件

Chain文件

格式见之前博客

用户输入文件

CrossMap 支持下面的格式，一般用 BED 。

• BAM, CRAM, or SAM
• BED or BED-like formats (The BED file must contain at least the fields: chrom, start, and end.)
• Wiggle (variableStep, fixedStep, and bedGraph formats are supported.)
• BigWig
• GFF or GTF
• VCF
• GVCF
• MAF

输出文件

输出文件格式取决于输入文件

Input Format	Output Format
BED	BED (genome coordinates will be updated)
BAM	BAM (genome coordinates, header section, all SAM flags, and insert sizes will be updated)
CRAM	BAM (requires pysam >= 0.8.2)
SAM	SAM (genome coordinates, header section, all SAM flags, and insert sizes will be updated)
Wiggle	BigWig
BigWig	BigWig
GFF	GFF (genome coordinates will be updated to the target assembly)
GTF	GTF (genome coordinates will be updated to the target assembly)
VCF	VCF (header section, genome coordinates, and reference alleles will be updated)
GVCF	GVCF (header section, genome coordinates, and reference alleles will be updated)
MAF	MAF (genome coordinates and reference alleles will be updated)

用法

转换 BED 格式文件

BED（Browser Extensible Data，浏览器可扩展数据）文件是一个制表符分隔的文本文件，用于描述基因组区域或基因注释信息。每一行代表一个特征，可包含 3 至 12 列。

CrossMap 通过更新染色体和坐标信息，同时保留其他字段不变，实现 BED 文件在不同参考基因组版本之间的转换。

• 对于少于 12 列的 BED 文件，仅更新染色体和坐标字段（不对吧，下面说了会更新链信息）。
• 对于12 列的 BED 文件，除以下列外，其余所有列均会更新（一般提供不了完整的12列）：
  • 第 4 列：特征名称（feature name）
  • 第 5 列：得分（score）
  • 第 9 列：RGB 显示颜色（RGB display color）

源基因组区域若映射到目标基因组的多个位置，则会被拆分为多个独立的条目。

输入输出

标准 BED 文件有 12列，但是 CrossMap 也支持 BED-like 文件（就是不是标准12列），如下

BED3: first three columns — chrom, start, end

BED6: first six columns — chrom, start, end, name, score, strand

Other BED-like formats: between 3 and 12 columns, where extra columns are arbitrary

注意事项：

1. 对于类似BED的格式，CrossMap仅更新chrom、start、end和strand字段，其余所有列原样保留。
2. 以#、browser或track开头的行将被忽略。
3. 列数少于三列的行将被跳过。
4. 第2列和第3列必须是整数。
5. 若未提供链(strand)信息，则默认假定为+链。
6. 对于标准的12列BED文件，如果任一外显子块无法唯一比对，则整条记录将被跳过。
7. 输入文件input_chain_file和input_bed_file可以是普通文件或压缩文件（扩展名如.gz、.Z、.z、.bz、.bz2、.bzip2），也可以是本地或远程URL。
8. 若未指定output_file，CrossMap会将结果输出到控制台；此时，原始（未转换）条目也会一并输出。
9. 如果某个输入区域无法连续地比对到目标参考基因组，则该区域会被拆分。
10. *.unmap文件包含无法明确转换的条目。

按照我个人的理解，如果采用 BED6 格式，则最后2列应该这么设置

• score列 ：全部设为0，CrossMap不会使用或修改这一列，没用。
• strand列 ：按照实际情况填写，重测序下机数据统一使用 + 链（因为肯定和其使用的参考基因组的 forward链一致）

CrossMap 软件会基于 chain 文件的 qStrand (查询链方向) 字段来决定是否翻转 BED6 格式中的链信息。如果 qStrand 为 - ，则链信息（+/-）会翻转。理由如下

CrossMap 处理 BED6 格式时，会根据比对情况自动调整链信息。其核心逻辑是：

• 如果 chain 文件中对应区域的 qStrand 为 -，则表示查询序列（即你的 BED 文件所在的原始基因组）在比对到目标基因组时发生了反向互补（reverse complement），也就是发生了倒位（inversion）。
• 此时，为了准确反映该区域在目标基因组上的真实方向，CrossMap 会翻转 BED6 记录中的链信息（+ 变 -，- 变 +）。

如果 qStrand 为 +，则表示方向一致，链信息保持不变。

这个机制确保了转换后 BED 文件中的链信息能准确代表该区域在新基因组上的实际方向。

命令

转换 bed 文件使用 CrossMap bed 命令，使用帮助见 CrossMap bed -h，如下。这基本很清晰了

usage: CrossMap bed [-h] [--chromid {a,s,l}] [--unmap-file UNMAP_FILE]
                    input.chain input.bed [out_bed]

positional arguments:
  input.chain           Chain file
                        (https://genome.ucsc.edu/goldenPath/help/chain.html)
                        Describes pairwise alignments between two genomes.
                        The chain file may be plain text or compressed
                        (.gz, .Z, .z, .bz, .bz2, .bzip2).

  input.bed             Input BED file. The first 3 columns must be
                        “chrom”, “start”, and “end”.
                        Supports plain text, compressed (.gz, .bz2, etc.),
                        or URL sources (http://, https://, ftp://).
                        Compressed remote files are not supported.

  out_bed               Output BED file.
                        If omitted, CrossMap writes to STDOUT.

options:
  -h, --help            Show this help message and exit
  --chromid {a,s,l}     Style of chromosome IDs:
                          a = as-is
                          l = long style (e.g., "chr1", "chrX")
                          s = short style (e.g., "1", "X")
  --unmap-file UNMAP_FILE
                        File to save unmapped entries.
                        Ignored if [out_bed] is not provided.

简单举例

CrossMap bed GRCh37_to_GRCh38.chain.gz Test1.hg19.bed output.hg38

如果输入序列不能匹配到目标序列的一段连续序列，那么输出结果会自动拆分，举例如下

$ cat Test2.hg19.bed
chr20  21106623   21227258
chr22  30792929   30821291

$ CrossMap bed GRCh37_to_GRCh38.chain.gz Test2.hg19.bed
2024-01-12 08:53:10 [INFO]  Read the chain file "GRCh37_to_GRCh38.chain.gz"
chr20  21106623   21227258  (split.1:chr20:21106623:21144549:+)  chr20  21125982  21163908
chr20  21106623   21227258  (split.2:chr20:21144549:21176014:+)  chr20  21163909  21195374
chr20  21106623   21227258  (split.3:chr20:21176014:21186161:+)  chr20  21195375  21205522
chr20  21106623   21227258  (split.4:chr20:21186161:21227258:+)  chr20  21205523  21246620
chr22  30792929   30821291  (split.1:chr22:30792929:30819612:+)  chr22  30396940  30423623
chr22  30792929   30821291  (split.2:chr22:30819615:30821291:+)  chr22  30423627  30425303

参考文献

1. https://crossmap.readthedocs.io/en/latest/#