随着基因组数据的不断完善,序列信息的不断更新,一个物种的基因组会有多个不同的组装版本。以人为例,目前就有hg16、hg17、hg18、hg19、hg38这几个组装版本。那么如何在这些不同版本之间进行坐标转化,以适应特定的分析需求呢?下面介绍CrossMap软件的使用。
1 基因组坐标转换常用软件
- NCBI的Remap:支持BED、VCF、GFF、GTF等。
- UCSC的LiftOver:支持BED,本地版本还支持GFF。
- CrossMap:BAM, CRAM, SAM, BED, Wiggle, BigWig, GFF, GTF and VCF等。
2 CrossMap安装与使用
2.1 安装
#安装(pip or conda)
conda install CrossMap
pip3 install CrossMap
2.2 使用
在使用之前先介绍下自己的需求。我有一部分小鼠的变异数据(BED格式),但是用的是mm9的坐标,需要把其转化为最新的mm10坐标。
2.2.1 BED格式文件介绍
BED文件格式提供了一种灵活的方式来定义的数据行,以用来描述注释的信息。BED行有3个必须的列和9个额外可选的列,每行的数据格式要求一致。
- chrom:染色体或scafflold的名字。
- chromStart:染色体或scaffold的起始位置(从零开始)。
- chromEnd:染色体或scaffold的结束位置,染色体的末端位置没有包含到显示信息里面。
2.2.2 chain文件下载
chain文件作为转换坐标时候的输入文件,可以根据你的物种名在UCSC上下载(这里以小鼠mm9转换为mm10为例,下载mm9ToMm10.over.chain.gz文件即可)。
3.2.3 坐标转换
CrossMap.py
#使用相对简单
CrossMap.py xx(什么文件格式:如bed) xx(chain文件) xx(原始坐标文件) xx(输出文件名)
#Note:由于我的bed文件莫名多出一列,这里先去掉一列,然后进行坐标转换。(注:bed文件里面第一行不能为列名,如果是请去掉)
awk '{$1="";print $0}' gain_1K_sign_6.txt >gain_1K_sign_6.bed
CrossMap.py bed mm9ToMm10.over.chain.gz gain_1K_sign_6.bed gain_1K_sign_6_mm10.bed
#生成的结果map到的结果为xx.bed,同时也生成unmap的结果。