如何使用Biopython将CSV文件中的snp/indels添加到FASTA文件中？

Question

我想修改一个FASTA文件中的序列。这个FASTA文件包含了人类基因组（每条染色体的序列），它的ID是 >1, >2,... >22, >X, >Y, >MT 和 >GL000207.1。

我想在每条染色体序列中引入的修改（突变）都在一个CSV文件里。这里有个例子：

chrom;position;ref;var;Gene;VAR
1;21424;C;T;WASH7P;snp.LOH
1;33252099;CACATGCATGACTATTCCTAGCC;-;YARS;indel_somatic
5;107061668;-;GT;EFNA5(dist=55072),FBXL17(dist=133066);indel_somatic
22;22677038;G;C;BMS1P20;snp_somatic
MT;16093;T;C;NONE(dist=NONE),NONE(dist=NONE);snp.LOH
X;22012649;-;T;SMS;indel_somatic

每一行描述了染色体的编号，以及在染色体上找到的snp/indel的位置。接下来的两列表示参考的核苷酸和需要插入到FASTA文件中的突变。这种修改可以是替换、删除（多个核苷酸）或插入（多个核苷酸）。最后两列不重要。输出应该是包含这些突变的新FASTA文件。

我已经创建了以下脚本。我知道我离想要的结果还有很远……我会努力改进，但在此期间，如果有人能给点建议，那就太好了。

from bisect import bisect_right
from collections import defaultdict
from Bio import SeqIO
from Bio.Seq import MutableSeq
from Bio.Alphabet import IUPAC
import csv

def line_to_snp(line):
    row = line.split(";")
    return row[0], int(row[1]), row[2], row[3], row[4], row[5]


with open('Modified_build.fasta', 'w') as f1:
   reference = SeqIO.read("human_g1k_v37.fasta", "fasta")
   datafile = open('snp_all.csv', 'r')
   snp = line_to_snp(line)
   for seq_record in SeqIO.parse(reference):
    mutable_seq = MutableSeq (reference, IUPACUnambigousDNA())
    if snp[0] == seq_record.id:
           mutable_seq[snp[1]] = snp[3]
           f1.write(seq_id)
           f1.write(seq_record)

Answer 1

你现在的做法是个不错的开始，但你的代码并没有对打开的CSV文件做任何处理（datafile没有被使用，line也是未定义的）。

我建议你从你的SNP文件中建立一个以访问键为基础的数据字典。你可以使用csv模块来读取用;分隔的文件。在遍历SeqRecord对象时，你可以从这个字典中获取突变数据。

一些杂项错误：

• MutableSeq不能修改记录，只能修改字符串或Seq。
• Python中的字符串（因此Seq对象）使用的是从零开始的索引，而序列索引是从一开始的。
• 未定义的名称IUPACUnambiguousDNA（这是的一部分）。

---

我的建议是使用下面这样的方式：

from Bio import SeqIO
from Bio.Alphabet import IUPAC
from Bio.Seq import MutableSeq
from Bio.SeqRecord import SeqRecord

import csv


# Build a dictionary of position & mutation by accession
snp_dict = {}
with open('snp_all.csv') as snp_all:
    csv_data = csv.reader(snp_all, delimiter=';')
    header = next(csv_data)
    for chrom, position, ref, var, _, _ in csv_data:
        snp_dict[chrom] = (int(position), ref, var)

# Iterate over (missing) FASTA file, create mutations
mutated_records = []
reference = SeqIO.read("human_g1k_v37.fasta", "fasta")
for record in SeqIO.parse("this_file_missing.fasta", "fasta"):
    mutable_seq = MutableSeq(reference.seq, IUPAC.IUPACUnambigousDNA())
    if record.id in snp_dict:  # Are all sequences mutated in snp_all?
        position, ref, var = snp_dict[record.id]
        mutable_seq[position - 1] = var
        record.seq = mutable_seq.toseq()  # Re-use the record- preserves id, etc.
        mutated_records.append(record)

with open('Modified_build.fasta', 'w') as f1:
    SeqIO.write(mutated_records, f1, "fasta")