简介
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