SYRI Synteny and Rearrangement Identifier

使用SYRI鉴定基因组变异及可视化

1.调整两条比对基因组的染色体方向

SYRI对比对的基因组有非常严格的要求：

1. 两个基因组的染色体ID必须一一对应(同源染色体ID必须一致，染色体数量也必须一致，染色体ID不能是数字)；
2. 同源染色体strand方向必须一致。

  我们根据之前已有的共线性关系，手动生成一个染色体操作表格，然后用pythpn代码进行ID替换与strand方向的改变。

  这里RAW_ID是其中一个样基因组染色体ID，NEW_ID是新的染色体ID，这个与另外一个样本的染色体ID是一一对应的(另外一个样本的染色体ID就是这种1,2,3,4..11), STRAND为链的方向，+代表两个样本这条染色体方向一致，-代表方向相反，需手动生成反向互补链。
替换表格:

RAW_ID	NEW_ID	STRAND
ctg000130	Chr1	+
ctg000050	Chr2	-
ctg000120	Chr3	+
ctg000020	Chr4	-
ctg000100	Chr5	+
ctg000040	Chr6	-
ctg000030	Chr7	+
ctg000110	Chr8	-
ctg000010	Chr9	+
ctg000000	Chr10	-
ctg000090	Chr11	-

from sys import argv, exit


def read_fa(fa_file):
    tmp_dic = []
    for line in open(fa_file):
        if line[0] == '\n':
            continue
        if line[0] == '>':
            tmp_dic.append([line[1:].split()[0], ''])
        else:
            tmp_dic[-1][1] += line.strip()
    tmp_dic = {t_lst[0]:t_lst[1] for t_lst in tmp_dic}

    return tmp_dic

def reverse_complement(dna_sequence):
    complement_dict = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C', 'N':'N'}
    reverse_sequence = dna_sequence[::-1]
    complement_sequence = ''.join(complement_dict[base] for base in reverse_sequence)
    return complement_sequence


if __name__ == '__main__':
    if len(argv) != 4:
        exit(f"python {argv[0]} [raw_fa] [info_tab] [paired_fa]")
    
    info_tab = {}
    for line in open(argv[2]):
        if line.startswith("RAW_ID"):
            continue
        tmp = line.strip().split()
        info_tab[tmp[0]] = tmp[1:]
 
    raw_fas = read_fa(argv[1])
    ouf_w = open(argv[3], 'w')
    for seq_id in info_tab:
        if info_tab[seq_id][1] == '+':
            tmp_seq = raw_fas[seq_id].upper()
        else:
            # 输出反向互补链
            tmp_seq = reverse_complement(raw_fas[seq_id].upper())
        ouf_w.write(">"+info_tab[seq_id][0]+"\n" + tmp_seq+"\n")
            
    ouf_w.close()
    

# 直接运行
python pair_chr.py ../nextDevono.fa info.tab gaisen.fa

2.SYRI鉴定SV并可视化

# nucmer比对
nucmer -c 100 -l 50 -t 8 -p Z7.gaisen Z7.fa gaisen.fa

# 比对结果过滤
delta-filter -m -i 90 -l 100 Z7.gaisen.delta > Z7.gaisen.filtered.delta

# 比对格式转换
show-coords -THrd Z7.gaisen.filtered.delta > Z7.gaisen.filtered.coords

# SYRI鉴定基因组变异
python3 /data/chaofan/software/miniconda/envs/syri_env/bin/syri -s /data/chaofan/software/miniconda/envs/syri_env/bin/show-snps -c Z7.gaisen.filtered.coords -d Z7.gaisen.filtered.delta -r Z7.fa -q gaisen.fa --prefix Z7.gaisen.

# 生成基因组信息文件，tab分隔(详细内容见代码框下)
vi genome.txt

# plotsr可视化 -H：输出图形高 -W： 输出图形高
plotsr \
    --sr Z7.gaisen.syri.out \
    --genomes genome.txt -o genome.pdf -b pdf -H 16 -W 12

genome.txt文件内容：

genome.txt
#file	name	tags
Z7.fa	Z7	lw:2
gaisen.fa	Gaisen	lw:2

  嗯，结果看着还行。