liftover——不同基因组座位转换

简介
The liftOver tool is useful if you wish to convert a large number of coordinate ranges between assemblies.

在线版
UCSC Genome browser中的Tools中有liftover的在线版工具Lift Genome Annotations
。在线版liftover可以进行相同物种不同版本基因组间的转换,也可以进行不同物种间的转换。

本地版
本地版的liftover只能在linux系统中以命令行的形式运行。

下载和安装必要工具地址与安装

1. 安装 UCSC Kent Utilities(包含 liftOver 等工具)

从源码编译或使用预编译版本

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtChain
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainNet
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/netChainSubset
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainSort
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainSplit
chmod +x liftOver axtChain chainNet netChainSubset chainSort chainSplit

2. 安装 LASTZ(全基因组比对工具)

方法一:从源码编译

wget https://github.com/lastz/lastz/archive/refs/tags/1.04.22.tar.gz

git clone https://github.com/lastz/lastz.git

tar -xzf 1.04.22.tar.gz
cd lastz-1.04.22/
make
make install

方法二:使用conda

conda install -c bioconda lastz ucsc-liftover

下载基因组序列
确保您有源基因组和目标基因组的 FASTA 文件:

创建项目目录结构

mkdir -p liftover_project/{genomes,alignment,chains,results}
cd liftover_project

假设您已经有以下基因组FASTA文件:

源基因组:11.fa, 22.fa, 33.fa, 44.fa

目标基因组:ref.fa

将文件放入 genomes 目录

自行生成链文件
步骤 1:准备基因组文件

为每个基因组创建染色体大小文件

samtools faidx genomes/97103v2.fa
cut -f1,2 genomes/97103v2.fa.fai > genomes/97103v2.chrom.sizes

samtools faidx genomes/USVL246-FR2.fa
cut -f1,2 genomes/USVL246-FR2.fa.fai > genomes/USVL246-FR2.chrom.sizes

samtools faidx genomes/PI296341.fa
cut -f1,2 genomes/PI296341.fa.fai > genomes/PI296341.chrom.sizes

samtools faidx genomes/CharlestonGray.fa
cut -f1,2 genomes/CharlestonGray.fa.fai > genomes/CharlestonGray.chrom.sizes

samtools faidx genomes/97103v1.fa
cut -f1,2 genomes/97103v1.fa.fai > genomes/97103v1.chrom.sizes

我以上都出现问题,后来发现minimap2进行基因组比对更加好用

下载安装minimap2
wget https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2
tar -xvf minimap2-2.30_x64-linux.tar.bz2
cd minimap2-2.30_x64-linux
./minimap2

Minimap2基因组比对
基本命令:
minimap2 -cx asm5 --cs=long -t 20
target_genome.fa
query_genome.fa
> alignment.paf
关键参数详解:
-cx asm5: 使用asm5预设,适用于近缘基因组比对
--cs=long: 输出详细的cs标签,便于后续处理
-t 20: 使用20个线程
--secondary=no: 不输出次要比对(可选)
-o alignment.paf: 指定输出文件

针对不同相似度的参数调整:

高相似度基因组(>95%)

minimap2 -cx asm5 --cs=long -t 20 target.fa query.fa > alignment.paf

中等相似度基因组(85-95%)

minimap2 -cx asm10 --cs=long -t 20 target.fa query.fa > alignment.paf

较低相似度基因组

minimap2 -cx asm20 --cs=long -t 20 target.fa query.fa > alignment.paf

PAF转Chain文件
paf2chainAndreaGuarracino/paf2chain: convert PAF format to CHAIN format

installation

paf2chain is built with rust, and so we install using cargo:

git clone https://github.com/AndreaGuarracino/paf2chain
cd paf2chain
cargo install --force --path .

报错没有cargo

CentOS/RHEL (EPEL required)

sudo yum install epel-release
sudo yum install rust cargo

使用kentUtils工具:

下载kentUtils工具集

git clone https://github.com/ucscGenomeBrowser/kent.git

编译所需工具

cd kent/src
make utils

转换PAF为Chain

pafToChain alignment.paf > alignment.chain

anyway,以上尝试各种出错:
最后的方案:
wget https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2
tar -xvf minimap2-2.30_x64-linux.tar.bz2
cd minimap2-2.30_x64-linux
./minimap2

Minimap2基因组比对
基本命令:
minimap2 -cx asm5 --cs=long -t 20
target_genome.fa
query_genome.fa
> alignment.paf
关键参数详解:
-cx asm5: 使用asm5预设,适用于近缘基因组比对
--cs=long: 输出详细的cs标签,便于后续处理
-t 20: 使用20个线程
--secondary=no: 不输出次要比对(可选)
-o alignment.paf: 指定输出文件

针对不同相似度的参数调整:

高相似度基因组(>95%)

minimap2 -cx asm5 --cs=long -t 20 target.fa query.fa > alignment.paf

中等相似度基因组(85-95%)

minimap2 -cx asm10 --cs=long -t 20 target.fa query.fa > alignment.paf

较低相似度基因组

minimap2 -cx asm20 --cs=long -t 20 target.fa query.fa > alignment.paf

查看不同基因组染色体命名方式
less minimap2_result/order_2_97103v2.paf
Chr01……

自己手动编辑bed文件
cat > order.bed
Chr02 12480191 13649064 qtl2-2
Chr05 13985119 14544500 qtl2-5
Chr09 16756994 17369072 qtl2-9
Chr10 18514495 18923367 qtl2-10
Chr08 10018545 11855293 qtl2-8

然后运行
k8 /root/software/minimap2-2.30_x64-linux/paftools.js liftover minimap2_result/order_2_refv2.paf order.bed > rder_2_refv2.bed

根据得分排序,以最优为锚定位点
sort -k4,4 -k5,5nr rder_2_refv2.bed | awk '!seen[$4]++' > rder_2_refv2_best.bed

©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容