寻找SV.vcf文件中是否含有某个片段的插入

现在已知这个片段在3号染色体上，基本思路是先将vcf文件变成fa文件，然后用blastn处理。

#!/bin/bash
python SVvcf2fa.py INS_SViper_out.vcf> All_INS_Mo17.fa
makeblastdb -in All_INS_Mo17.fa -dbtype nucl -title All -input_type fasta -max_file_sz 1GB -parse_seqids -out blastdb
blastn -query ins.fa -db blastdb -outfmt 6 -num_threads 1 >blast_map_SV2all_result

第四列记录的是INS的序列内容

vcf文件

#!/bin/python
import sys
vcf_fi = open(sys.argv[1], "rt")

vcf = vcf_fi.readline()
pos = ""

while vcf:
    if vcf.startswith("##"):
        pass
    # FAM存了也没用
    elif vcf.startswith("#CHROM"):
        header = vcf.rstrip().split("\t")[9::]
    # 记录SV_seq和pos
    elif vcf.startswith("chr3"):
        vcf_line = vcf.rstrip().split("\t")

        if pos == f">{vcf_line[0]}_{vcf_line[1]}":
            #  处理重复header
            print(f"{pos}_dup2")
        else:
            pos = f">{vcf_line[0]}_{vcf_line[1]}"
            print(pos)
        seq = vcf_line[4]
        print(seq)
    vcf = vcf_fi.readline()

vcf_fi.close()

在实际操作过程中竟然出现了重复三次，并且在vcf文件中不连续的pos。一方面是可以选择手动把第三个改成dup3，另一方面就是只选取3号染色体上的SV。因为最终比对结果中也没有出现这些少量的重复位置。

blastn结果