2

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

1 回答 1

0

您当前的方法是一个好的开始，但是您的代码对打开的 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")
于 2014-03-24T20:03:06.840 回答