不同基因组之间的坐标转换

基本思路都是从不同软件的结果通过手动操作转为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.更多其他比对软件待续