现在已知这个片段在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的序列内容
#!/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。因为最终比对结果中也没有出现这些少量的重复位置。