使用python根据SNP变异检测vcf文件中的碱基突变信息并替换参考基因组中相应位置的单碱基

001、

(base) root@PC1:/home/test2# ls
a.fasta  b.vcf  test.py
(base) root@PC1:/home/test2# cat a.fasta                          ## 测试fasta文件
>NC_000964.3 Bacillus subtilis subsp. subtilis str. 168 chromosome, complete genome
ATCGTTTTCGGCTTTTTTTAGTATCCACAGAGGTTATCGACAACATTTTCACATTACCAACCCCTGTGGACAAGGTTTTT
TCAACAGGTTGTCCGCTTTGTGGATAAGATTGTGACAACCATTGCAAGCTCTCGTTTATTTTGGTATTATATTTGTGTTT
TAACTCTTGATTACTAATCCTACCTTTCCTCTTTATCCACAAAGTGTGGATAAGTTGTGGATTGATTTCACACAGCTTGT
GTAGAAGGTTGTCCACAAGTTGTGAAATTTGTCGAAAAGCTATTTATCTACTATATTATATGTTTTCAACATTTAATGTG
TACGAATGGTAAGCGCCATTTGCTCTTTTTTTGTGTTCTATAACAGAGAAAGACGCCATTTTCTAAGAAAAGGAGGGACG
TGCCGGAAGATGGAAAATATATTAGACCTGTGGAACCAAGCCCTTGCTCAAATCGAAAAAAAGTTGAGCAAACCGAGTTT
TGAGACTTGGATGAAGTCAACCAAAGCCCACTCACTGCAAGGCGATACATTAACAATCACGGCTCCCAATGAATTTGCCA
GAGACTGGCTGGAGTCCAGATACTTGCATCTGATTGCAGATACTATATATGAATTAACCGGGGAAGAATTGAGCATTAAG
TTTGTCATTCCTCAAAATCAAGATGTTGAGGACTTTATGCCGAAACCGCAAGTCAAAAAAGCGGTCAAAGAAGATACATC
TGATTTTCCTCAAAATATGCTCAATCCAAAATATACTTTTGATACTTTTGTCATCGGATCTGGAAACCGATTTGCACATG
(base) root@PC1:/home/test2# cat b.vcf                                   ## 测试vcf文件
##reference=file://Bacillus_subtilis.str168.fasta
##source=SelectVariants
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Bacillus_subtilis
NC_000964.3     1       .       A       C       3487    PASS    AC=1;AF=1.00;AN=1;DP=87;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=34.24;SOR=0.812      GT:AD:DP:GQ:PL  1:0,87:87:99:3517,0
NC_000964.3     2       .       T       A       3529    PASS    AC=1;AF=1.00;AN=1;DP=89;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=30.63;SOR=0.809      GT:AD:DP:GQ:PL  1:0,89:89:99:3559,0
NC_000964.3     3       .       C       T       3621    PASS    AC=1;AF=1.00;AN=1;DP=93;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=29.09;SOR=0.714      GT:AD:DP:GQ:PL  1:0,93:93:99:3651,0
NC_000964.3     4       .       G       T       3272    PASS    AC=1;AF=1.00;AN=1;DP=86;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=32.93;SOR=0.765      GT:AD:DP:GQ:PL  1:0,85:85:99:3302,0
NC_000964.3     111     .       A       G       3062    PASS    AC=1;AF=1.00;AN=1;DP=75;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=30.83;SOR=0.832      GT:AD:DP:GQ:PL  1:0,75:75:99:3092,0
NC_000964.3     115     .       G       A       2962    PASS    AC=1;AF=1.00;AN=1;DP=72;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=29.67;SOR=0.807      GT:AD:DP:GQ:PL  1:0,72:72:99:2992,0
NC_000964.3     118     .       A       G       2805    PASS    AC=1;AF=1.00;AN=1;DP=71;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=60.00;QD=34.04;SOR=0.779      GT:AD:DP:GQ:PL  1:0,71:71:99:2835,0
(base) root@PC1:/home/test2# cat test.py                                ## 测试脚本
#!/usr/bin/python

in_file1 = open("a.fasta", "r")
in_file2 = open("b.vcf", "r")
out_file = open("result.txt", "w")

seq_all = {}
fasta_id_all = []
scaff_id_all = []

for i in in_file1:
    i = i.strip()
    if i[0] == ">":
        fasta_id_all.append(i)
        seq = i.split()[0].split(">")[1]
        scaff_id_all.append(seq)
        seq_all[seq] = []
    else:
        for j in i:
            seq_all[seq].append(j)

for i in in_file2:
    i = i.strip()
    if i[0] == "#":
        continue
    else:
        i = i.split("\t")
        seq_all[i[0]][int(i[1]) - 1] = i[4]

for i in scaff_id_all:
    print(fasta_id_all[scaff_id_all.index(i)], file = out_file)
    print("".join(seq_all[i]), file = out_file)

in_file1.close()
in_file2.close()
out_file.close()
(base) root@PC1:/home/test2# python test.py                        ## 执行程序
(base) root@PC1:/home/test2# ls
a.fasta  b.vcf  result.txt  test.py
(base) root@PC1:/home/test2# cat result.txt                        ## 程序运行 结果
>NC_000964.3 Bacillus subtilis subsp. subtilis str. 168 chromosome, complete genome
CATTTTTTCGGCTTTTTTTAGTATCCACAGAGGTTATCGACAACATTTTCACATTACCAACCCCTGTGGACAAGGTTTTTTCAACAGGTTGTCCGCTTTGTGGATAAGATGGTGACAGCCATTGCAAGCTCTCGTTTATTTTGGTATTATATTTGTGTTTTAACTCTTGATTACTAATCCTACCTTTCCTCTTTATCCACAAAGTGTGGATAAGTTGTGGATTGATTTCACACAGCTTGTGTAGAAGGTTGTCCACAAGTTGTGAAATTTGTCGAAAAGCTATTTATCTACTATATTATATGTTTTCAACATTTAATGTGTACGAATGGTAAGCGCCATTTGCTCTTTTTTTGTGTTCTATAACAGAGAAAGACGCCATTTTCTAAGAAAAGGAGGGACGTGCCGGAAGATGGAAAATATATTAGACCTGTGGAACCAAGCCCTTGCTCAAATCGAAAAAAAGTTGAGCAAACCGAGTTTTGAGACTTGGATGAAGTCAACCAAAGCCCACTCACTGCAAGGCGATACATTAACAATCACGGCTCCCAATGAATTTGCCAGAGACTGGCTGGAGTCCAGATACTTGCATCTGATTGCAGATACTATATATGAATTAACCGGGGAAGAATTGAGCATTAAGTTTGTCATTCCTCAAAATCAAGATGTTGAGGACTTTATGCCGAAACCGCAAGTCAAAAAAGCGGTCAAAGAAGATACATCTGATTTTCCTCAAAATATGCTCAATCCAAAATATACTTTTGATACTTTTGTCATCGGATCTGGAAACCGATTTGCACATG

002、验证

参考：https://mp.weixin.qq.com/s?__biz=MzIxNzc1Mzk3NQ==&mid=2247491511&idx=1&sn=0d329538b90decbd8f7c9d6164cd2ae9&chksm=97f5afafa08226b9a38ddfaad521285d07fe5bae99ade4327ba6c39c274f807cf6ac6fc122e5&scene=178&cur_album_id=2403674812188688386#rd