利用已知序列从fasta文件中提取序列和头文件

Question

我试图做比较两个文件和提取序列有其他子集。我还想提取标识符。但是，我能做的是能够提取包括子集的序列。示例文件包括：

text.fa
>header1
ETTTHAASCISATTVQEQ*TLFRLLP
>header2
SKSPCSDSDY**AAA
>header3
SSGAVAAAPTTA

而且

^{pr2}$

当我运行代码时，我得到以下输出：

ETTTHAASCISATTVQEQ*TLFRLLP
SSGAVAAAPTTA

但是，我的预期输出是带标题的：

>header1
ETTTHAASCISATTVQEQ*TLFRLLP
>header3
SSGAVAAAPTTA

我的代码分为两部分，首先用这些序列创建文件，然后尝试从原始fasta文件中提取它们的头：

def get_nucl(filename):
    with open(filename,'r') as fd:
        nucl = []
        for line in fd:
            if line[0]!='>':
                nucl.append(line.strip())
        return nucl
def finding(filename,reffile):
        nucl = get_nucl(filename)
        with open(reffile,'r') as reffile2:
            for line in reffile2:
                for element in nucl:
                    if line.strip() in element:
                            yield(element)



    with open('sequencesmatched.txt','w') as output:
            results = finding('text.fa','textref.fa',)
            for res in results:
                print(res)
                output.write(res + '\n')

因此，在这个sequencesmatched.txt中，我得到了text.fa的序列，它们的子串是{}。作为：

ETTTHAASCISATTVQEQ*TLFRLLP
SSGAVAAAPTTA

所以在另一部分中，要检索相应的标题和这些序列：

    def finding(filename,seqfile):
        with open(filename,'r') as fastafile:
                with open(seqfile,'r') as sequf:
                        alls=[]
                        for line in fastafile:
                                alls.append(line.strip())
                        print(alls)
                        sequfs = []
                        for line2 in sequf:
                                sequfs.append(line2.strip())
                                if str(line.strip()) == str(line2.strip()):
                                        num = alls.index(line.strip())
                                        print(alls[num-1] + line)


print(finding('text.fa','sequencesmatched.txt'))

但是，作为输出，我只能检索一个序列，这是第一个匹配：

>header1
ETTTHAASCISATTVQEQ*TLFRLLP

也许我可以不使用第二个文件，但是我不能进行正确的循环来获取序列和它们各自的头。所以，我走了很远的路。。在

如果你能帮忙我会很高兴的！在

Answer 1

如果文件的结构始终相同，则可以更轻松地执行以下操作：

def get_nucl(filename):
    with open(filename, 'r') as fd:
        headers = {}
        key = ''
        for line in fd.readlines():    
            if '>' in line:
                key = line.strip()[1:] # to remove the '>'
            else:
                headers[key] = line.strip()

    return headers 

这里我假设你的文件以“>；headern”开头，如果不是的话，你必须添加一些测试。现在你有一个像headers['header1'] = 'ETTTHAASCISATTVQEQ*TLFRLLP'这样的单词。在

所以现在要找到匹配项，你只需使用这句话：

^{pr2}$

因此，当您有一个头与它们的值匹配的dict时，如果您有一个子字符串，并且已经将头值作为键，那么就可以签入dict。在

刚才看到你做了print(finding(....)，你的函数已经打印出来了，所以就调用它。在