在一个 fasta 中查找 contig 的长度，使用另一个 fasta 的标头作为 python 中的查询答案

【问题标题】：Find length of a contig in one fasta, using the header of another fasta as query in python在一个 fasta 中查找 contig 的长度，使用另一个 fasta 的标头作为 python 中的查询
【发布时间】：2020-07-26 22:29:24
【问题描述】：

我正在尝试找到一个 python 解决方案，以使用序列的完整标题作为查询来提取 fasta 文件中特定序列的长度。完整的标头存储为管道中较早的变量（即“CONTIG”）。我想将此脚本的输出保存为变量，以便稍后在同一管道中使用。

以下是使用 Lucía Balestrazzi 提供的代码的脚本的更新版本。

附加信息：以下 with 语句嵌套在一个较大的 for 循环中，该循环遍历原始基因组的子样本。我目录中的第一个子样本 fasta 有一个长度为 40129801 的单个序列“>chr1:0-40129801”。我正在尝试写出一个文本文件“OUTPUT”，其中包含有关每个子样本 fasta 的一些基本信息。此文本文件将用作下游另一个程序的输入。

原始 fasta 文件中的标题名称是 chr1、chr2 等...而子示例 fastas 中的标题名称类似于：

batch1.fa >chr1:0-40k

batch2.fa >chr1:40k-80k

...等等...

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse(ORIGINAL_GENOME, "fasta")) #not the subsample
with open(GENOME_SUBSAMPLE, 'r') as FIN: 
    for LINE in FIN:
        if LINE.startswith('>'):
                #Example of "LINE"... >chr1:0-40129801
                HEADER = re.sub('>','',LINE)
                #HEADER = chr1:0-40129801
                HEADER2 = re.sub('\n','',HEADER)
                #HEADER2 = chr1:0-40129801 (no return character on the end)
                CONTIG = HEADER2.split(":")[0]
                #CONTIG = chr1 
                PART2_HEADER = HEADER2.split(":")[1]
                #PART2_HEADER = 0-40129801 
                START = int(PART2_HEADER.split("-")[0])
                #START = 0
                END = int(PART2_HEADER.split("-")[1])
                #END = 40129801
                LENGTH = END-START
                #LENGTH = 40129801 minus 0 = 40129801
            
                #This is where I'm stuck...
                ORIGINAL_CONTIG_LENGTH = len(record_dict[CONTIG]) #This returns "KeyError: 1"
                #ORIGINAL_CONTIG_LENGTH = 223705999 (this is from the full genome, not the subsample). 

                OUTPUT.write(str(START) + '\t' + str(HEADER2) + '\t' + str(LENGTH) + '\t' + str(CONTIG) + '\t' + str(ORIGINAL_CONTIG_LENGTH) + '\n')
                #OUTPUT = 0    chr1:0-40129801    40129801    chr1    223705999
OUTPUT.close()

我对生物信息学比较陌生。我知道我搞砸了我如何使用字典，但我不太确定如何解决它。

任何建议将不胜感激。谢谢！

【问题讨论】：

标签： python bioinformatics biopython fasta

【解决方案1】：

你可以这样做：

import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("genome.fa", "fasta"))
print(len(record_dict["chr1"]))

或

import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("genome.fa", "fasta"))
seq = record_dict["chr1"]
print(len(seq))

编辑：替代代码

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse("genome.fa", "fasta")
names = record_dict.keys()
for HEADER in names:
    #HEADER = chr1:0-40129801
    ORIGINAL_CONTIG_LENGTH = len(record_dict[HEADER])
    CONTIG = HEADER.split(":")[0]
    #CONTIG = chr1 
    PART2_HEADER = HEADER.split(":")[1]
    #PART2_HEADER = 0-40129801 
    START = int(PART2_HEADER.split("-")[0])
    END = int(PART2_HEADER.split("-")[1])
    LENGTH = END-START

这个想法是您定义一次字典，获取其键的值（所有 contigs 标头）并将它们存储为变量，然后循环通过标头提取您需要的信息。无需循环文件。

干杯

【讨论】：

这太棒了！无论如何我可以使用一个变量作为标题名称，所以我不必使用“chr1”......CONTIG = chr1...... print(len(record_dict[CONTIG]))
@Gunther 我不确定我是否理解正确你想要什么。如果要打印所有重叠群的长度而不一一命名，则在定义字典后，您可以将所有名称放在一个带有names = record_dict.keys()的变量中，然后循环提取长度的名称：for n in names: print(len(record_dict[n])) .您还可以在打印函数中添加 contig 的名称：for n in names: print(n, len(record_dict[n]))。希望对您有所帮助！
嗨@Lucía！到目前为止，您在引导我找到正确答案方面提供了巨大的帮助。我的脚本仍然有问题。我在原始问题中添加了更多信息。如果您有时间提供任何帮助，我将不胜感激。谢谢！
@Gunther 我一直在查看您的代码。 de dict 确实存在问题，您只需定义一次，而不是在每个循环中。让我发布一个替代代码，看看它是否满足您的需求。
我非常感谢你愿意帮助我，Lucía。对不起，我最初的问题不清楚。我在上面添加了更多上下文。我正在尝试根据原始基因组的较小子样本的标题提取一个基因组（原始）中的序列长度。

【解决方案2】：

这可行，只是将“CONTIG”变量更改为字符串。感谢 Lucía 在过去几天提供的所有帮助！

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse(ORIGINAL_GENOME, "fasta")) #not the subsample
with open(GENOME_SUBSAMPLE, 'r') as FIN: 
    for LINE in FIN:
        if LINE.startswith('>'):
                #Example of "LINE"... >chr1:0-40129801
                HEADER = re.sub('>','',LINE)
                #HEADER = chr1:0-40129801
                HEADER2 = re.sub('\n','',HEADER)
                #HEADER2 = chr1:0-40129801 (no return character on the end)
                CONTIG = HEADER2.split(":")[0]
                #CONTIG = chr1 
                PART2_HEADER = HEADER2.split(":")[1]
                #PART2_HEADER = 0-40129801 
                START = int(PART2_HEADER.split("-")[0])
                #START = 0
                END = int(PART2_HEADER.split("-")[1])
                #END = 40129801
                LENGTH = END-START
                #LENGTH = 40129801 minus 0 = 40129801
            
                #This is where I'm stuck...
                ORIGINAL_CONTIG_LENGTH = len(record_dict[str(CONTIG)]) 
                #ORIGINAL_CONTIG_LENGTH = 223705999 (this is from the full genome, not the subsample). 

                OUTPUT.write(str(START) + '\t' + str(HEADER2) + '\t' + str(LENGTH) + '\t' + str(CONTIG) + '\t' + str(ORIGINAL_CONTIG_LENGTH) + '\n')
                #OUTPUT = 0    chr1:0-40129801    40129801    chr1    223705999
OUTPUT.close()

【讨论】：