基本思路都是从不同软件的结果通过手动操作转为minimap2的paf格式
再转为chain文件，基于chain文件便能获得坐标转换
paf文件格式解读:https://wap.sciencenet.cn/blog-994715-1341509.html
1.基于Minimap2比对
minimap2 -cx asm5 --cs ref.fasta query.fasta > query2ref.paf
transanno minimap2chain query2ref.paf --output query2ref.chain
有了chain文件之后很多软件就都能做坐标转换, 比如vcf文件
transanno liftvcf --chain query2ref.chain --output query2ref.vcf --fail fail.vcf --new-assembly query.fa --original-assembly ref.fa --vcf old.vcf
2.基于nucmer比对
nucmer --mum --mincluster 500 -t 30 ref.fasta query.fasta -p ref_query
delta-filter -m -i 90 -l 100 ref_query.delta > ref_query.filter.delta
show-coords -c -r ref_query.filter.delta > ref_query.filter.coords
sed '1,5d' galiliv3_golaniv4.filter.Chr.coords | awk 'BEGIN{ OFS="\t" }{if($4<$5) print $16,"q.length",$4-1,$5,"+",$15,"r.length",$1-1,$2,$5-$4+1,$5-$4+1,'60',"cg:Z:"$5-$4+1"M","cs:Z::"$5-$4+1; else print $16,"q.length",$4-1,$5,"-",$15,"r.length",$1-1,$2,($5-$4+1)*-1,($5-$4+1)*-1,'60',"cg:Z:"($5-$4+1)*-1"M","cs:Z::"($5-$4+1)*-1}' > tmp1
#tmp1总共14列, 此时还有2列信息(第2列、、第7列)需要补充
bioawk -c fastx '{print $name, length($seq)}' ref.fasta > ref.length
bioawk -c fastx '{print $name, length($seq)}' query.fasta > query.length
#填充了第2列
python 1.py query.length tmp1 | awk '{$3="";$4="";print $0}' | sed 's/\s\+/ /g' | tr " " "\t" > tmp2
#填充了第7列
python 2.py ref.length tmp2 | awk 'BEGIN{ FS="\t";OFS="\t" }{print $3,$4,$5,$6,$7,$1,$2,$10,$11,$12,$13,$14,$15,$16}' > nucmer.paf
transanno minimap2chain nucmer.paf --output query2ref.nucmer.chain
2.5 1.py的内容
import sys
import re
list1 = {}
#f = open("tmp.sam",'r')
f = open(sys.argv[1],'r')
for line in f:
line = line.strip()
content = line.split('\t') #具体分隔符，具体第几列
name = content[0]
list1[name] = line
f = open(sys.argv[2],'r') #用传入参数就这样
for line in f:
line = line.strip()
content1 = line.split()
name1 = content1[0]
#print(content)
if name1 in list1:
print(list1[name1] + "\t" + line)
2.6 2.py的内容，和1.py是一样的
import sys
import re
list1 = {}
f = open(sys.argv[1],'r') #文件多的话就用传入参数
for line in f:
line = line.strip()
content = line.split('\t') #具体分隔符，具体第几列
name = content[0]
list1[name] = line
f = open(sys.argv[2],'r') #用传入参数就这样
for line in f:
line = line.strip()
content1 = line.split()
name1 = content1[5]
#print(content)
if name1 in list1:
print(list1[name1] + "\t" + line)
3.基于挂载后的contig基因组和染色体级别基因组的坐标对应文件
3D-DNA或者yahs生成的坐标对应文件.png
#这个脚本处理3DDNA输出的结果，可以得到 *.agp，也就是contig级别和染色体级别基因组的位置对应关系
#python juicer_assembly2agp_fa.py review.assembly ctg.fa Chr 32 0
#.agp文件转换vcf
grep -w -v 'N' Galili.Chr.agp | grep -v '#' > agp.1
python 1.py galili.ctg.length agp.1 > agp.2 #contig级别基因组的长度信息, 每一行是这样的: ptg00001 101241
python 2.py Galili.Chr.length agp.2 > agp.3 #染色体级别基因组的长度信息, 这样制取bioawk -c fastx '{print $name, length($seq)}' galili.ctg.fa > galili.ctg.length
awk 'BEGIN{ FS="\t";OFS="\t" }{print $6,$11,$7-1,$8,$9,$1,$13,$2-1,$3,$8-$7+1,$8-$7+1,'60',"cg:Z:"$8-$7+1"M","cs:Z::"$8-$7+1}' agp.3 | sort > out.paf
rm agp.1 agp.2 agp.3
/data/01/user164/software/transanno-x86_64-unknown-linux-musl-v0.3.0/transanno minimap2chain out.paf --output out.chain
transanno liftvcf --chain out.chain --output final.vcf --fail fail.vcf --new-assembly query.fa --original-assembly ref.fa --vcf old.vcf
CrossMap.py gff carmeli.chain carmeli.ctg.gff | awk -F "->" '{print $2}' | sed 's/^[\t ]\+//' > carmeli.chr.tmp
bedtools sort -i carmeli.chr.tmp > carmeli.chr.gff
1.py
#!/usr/bin/python3
import sys
import re
list1 = {}
f = open(sys.argv[1],'r') #文件多的话就用传入参数
for line in f:
line = line.strip()
content = line.split('\t') #具体分隔符，具体第几列
name = content[0]
list1[name] = line
f = open(sys.argv[2],'r') #用传入参数就这样
for line in f:
line = line.strip()
content1 = line.split('\t')
name1 = content1[5]
if name1 in list1:
print(line + "\t" + list1[name1])
#print(list1[name1])
2.py
#!/usr/bin/python3
import sys
import re
list1 = {}
f = open(sys.argv[1],'r') #文件多的话就用传入参数
for line in f:
line = line.strip()
content = line.split('\t') #具体分隔符，具体第几列
name = content[0]
list1[name] = line
f = open(sys.argv[2],'r') #用传入参数就这样
for line in f:
line = line.strip()
content1 = line.split('\t')
name1 = content1[0]
if name1 in list1:
print(line + "\t" + list1[name1])
https://github.com/hillerlab/make_lastz_chains
非常方便快捷的一个软件, 输入两个软屏蔽的fasta, 一键即可运行完毕
实测运行时间大概在1-3天, 当然是因为CPU给的比较多
make_chains.py target query target.fa query.fa --project_dir target2query -f --chaining_memory 70
make_chains.py contig Chr contig.fa Chr.fa --project_dir contig2Chr -f --chaining_memory 70
一行即可: 把target的坐标转换到query上
