現(xiàn)在已知這個片段在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的序列內(nèi)容
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]}":
# 處理重復(fù)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()
在實際操作過程中竟然出現(xiàn)了重復(fù)三次可款，并且在vcf文件中不連續(xù)的pos。一方面是可以選擇手動把第三個改成dup3末盔，另一方面就是只選取3號染色體上的SV筑舅。因為最終比對結(jié)果中也沒有出現(xiàn)這些少量的重復(fù)位置。
blastn結(jié)果
