首页 文章

通过已知序列从fasta文件中提取序列和头

提问于
浏览 次
2

我试图比较两个文件并提取具有其他子集的序列 . 而且，我也想提取标识符 . 但是，我能做的是能够提取包括子集的序列 . 示例文件是：

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

和，

textref.fa
>textref.fa
CISA
AAAP
AATP

当我运行代码时，我有这个输出：

ETTTHAASCISATTVQEQ*TLFRLLP
SSGAVAAAPTTA

但是，我的预期输出是 Headers ：

>header1
ETTTHAASCISATTVQEQ*TLFRLLP
>header3
SSGAVAAAPTTA

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

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 的序列具有 textref.fa 的子串 . 如：

ETTTHAASCISATTVQEQ*TLFRLLP
SSGAVAAAPTTA

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

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

也许我可以在没有第二个文件的情况下完成它，但我无法制作正确的循环来获取序列及其各自的 Headers . 因此，我走了很长的路..

如果你能提供帮助，我会很高兴的！

python indexing fasta

1 回答

  • 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' 这样的词典 .

    所以现在找到你刚刚使用那个dict的匹配：

    def finding(filename, reffile):
    headers = get_nucl(filename)
    with open(reffile, 'r') as f:
        matches = {}
        for line in f.readlines():
            for key, value in headers.items():
                if line.stip() in value and key not in matches:
                    matches[key] = value

    return matches

    因此，当您有一个带有与其值匹配的标头的dict时，您可以只检查dict是否有子字符串，并且您已经将标头值作为键 .

    刚刚看到你做了 print(finding(....) ，你的功能已经打印出来，所以就这么说吧 .

    回复于

相关问题