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

【问题标题】：How to add snp/indels from a CSV file to a FASTA file using Biopython?如何使用 Biopython 将 CSV 文件中的 snp/indels 添加到 FASTA 文件中？
【发布时间】：2014-04-29 09:14:36
【问题描述】：

我想修改 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)

【问题讨论】：

你有没有尝试解决这个问题？你有任何代码要显示吗？ Biopython tutorial 应该是您熟悉读/写 FASTA 记录的良好起点。
是的，我已经添加了到目前为止我所做的...
好的，太好了！现在，您具体遇到了什么问题？
我不知道如何专门为 csv 文件中明确标识的一个核苷酸更改 fasta 文件的核苷酸。如果我们看上面的例子，我想改变fasta文件1号染色体中T的核苷酸21424

标签： python csv biopython fasta mutation

【解决方案1】：

您当前的方法是一个好的开始，但是您的代码对打开的 CSV 文件没有任何作用（datafile 未被触及，line 未定义）。

我会从您的 SNP 文件中构建一个加入键控的数据字典。您可以使用csv 模块读取;-delimited 文件。在遍历SeqRecord 对象时，您可以从该字典中获取变异数据。

其他错误：

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

我的建议是使用如下方法：

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")

【讨论】：