2

我有一个fasta文件。从该文件中，我需要获取在序列的末尾和/或开头包含GTACAGTAGGand的唯一序列，并将它们放入新的 fasta 文件中。CAACGGTTTTGCC所以这里有一个例子：

>m121012_054644_42133_c100390582550000001523038311021245_s1_p0/7/2516_3269
***GTACAGTAGG***GTACACACAGAACGCGACAAGGCCAGGCGCTGGAGGAACTCCAGCAGCTAGATGCAAGCGACTA
TCAGAGCGTTGGGTCCAGAACGAAGAACAGTCACTCAAGACTGCTTT***CAACGGTTTTGCC***

>m121012_054644_42133_c100390582550000001523038311021245_s1_p0/7/3312_3597
CGCGGCATCGAATTAATACGACTCACTATAGGTTTTTTTATTGGTATTTTCAGTTAGATTCTTTCTTCTTAGAGGGTACA
GAGAAAGGGAGAAAATAGCTACAGACATGGGAGTGAAAGGTAGGAAGAAGAGCGAAGCAGACATTATTCA

>m121012_054644_42133_c100390582550000001523038311021245_s1_p0/7/3708_4657
***CAACGGTTTTGCC***ACAAGATCAGGAACATAAGTCACCAGACTCAATTCATCCCCATAAGACCTCGGACCTCTCA
ATCCTCGAATTAGGATGTTCTCCCCATGGCGTACGGTCTATCAGTATATAAACCTGACATACTATAAAAAAGTATACCAT
TCTTATCATGTACAGTAGG***GTACAGTAGG***

>m121012_054644_42133_c100390582550000001523038311021245_s1_p0/7/4704_5021
***GTACAGTAGG***GTGGGAGAGATGGCAGAAAGGCAGAAAGGAGAAAGATTCAGGATAACTCTCCTGGAGGGGCGAG
GTGCCATTCCCTGTGGTCACTTATTCTAAAGGCCCCAACCCTTCAAC***CAACGGTTTTGCC***

>m121012_054644_42133_c100390582550000001523038311021245_s1_p0/8/4223_4358
AAATATTGGGTCAAAGAACCGTTACTTTTCTTATATATGCGGCGCGAGGTTTTATATACTGATAAGAACCTACGCCATGG
GACATCTAATTCAGAGGGAAGAAGGTCCATGTCTGTTTGGATGAAATTGAGTCTG

（*添加用于突出显示）

我需要一些方法来获取在序列的末尾和/或开头包含 GTACAGTAGG 和 CAACGGTTTTGCC 的唯一序列，并将它们放在新的 fasta 文件中。我对此很陌生。我什至不确定是否可以做到。提前感谢您提供的任何帮助。

4

4 回答 4

1

虽然检查字符串在结尾和/或开头是否具有特定序列对于 Python 来说确实是一项非常简单的任务（请参阅str.startswith和str.endswith），但这并不能解决从 FASTA 文件中提取序列作为字符串的问题。这里有一些问题，例如序列必须与其注释分开，而且它们可以跨越多行。因此，将字符串方法直接应用于文件的行不会产生所需的结果。

这就是为什么您需要一个实际的 FASTA 格式解析器。解析器将处理一个 FASTA 文件，并将注释和序列作为 Python 字符串提供给您。BioPython 确实提供了一个，你可以这样做：

from Bio import SeqIO
def filter_records(source, substrings):
    for rec in source:
        if any(rec.seq.startswith(sub) or rec.seq.endswith(sub)
                 for sub in substrings):
             yield rec

substrings = ('GTACAGTAGG', 'CAACGGTTTTGCC')
SeqIO.write(filter_recors(SeqIO.parse('my.fasta', 'fasta'),
            'filtered.fasta', 'fasta')

我还可以推荐使用pyteomics（我参与开发的基于 Python 的蛋白质组学微框架）来操作 FASTA 文件：

from pyteomics import fasta
def filter_fasta(source, substrings):
    for descr, seq in source:
        if any(seq.startswith(sub) or seq.endswith(sub) for sub in substrings):
            yield descr, seq

substrings = ('GTACAGTAGG', 'CAACGGTTTTGCC')
fasta.write(filter_fasta(fasta.read('my.fasta'), substrings), 'filtered.fasta')
于 2012-12-14T21:29:00.213 回答
1

这是在 Biopython 中执行此操作的一种方法：

from Bio import SeqIO

source = 'fasta_file_name.fa'
outfile = 'filtered.fa'
sub1 ='GTACAGTAGG'
sub2 = 'CAACGGTTTTGCC'

def seq_check(seq, sub1, sub2):
    # basically a function to check whether seq starts or ends with sub1 or sub2
    return seq.startswith(sub1) or seq.startswith(sub2) or \
        seq.endswith(sub1) or seq.endswith(sub2)

seqs = SeqIO.parse(source, 'fasta')
filtered = (seq for seq in seqs if seq_check(seq.seq, sub1, sub2))
SeqIO.write(filtered, outfile, 'fasta')

我们正在使用生成器表达式（以 'filtered' 开头的行），因此过滤是在程序读取源文件时完成的。这具有节省内存的优点。

我们还创建了一个新函数来检查序列的开始和结束，以使程序更具可读性。理论上，您可以在生成器表达式中进行检查，但这会使行变得不必要地长。

希望有帮助:)。

于 2012-12-15T03:37:14.913 回答
1

可能不是最好的方法，具体取决于序列的大小，但这将完成工作。

import re
data_file ="location_of_fasta_file"
sequence = ''
Valid = False
for line in open(data_file):
    line = line.rstrip()
    if re.match("^>",line):
        if re.findall('^GTACAGTAGG',sequence) or re.findall('GTACAGTAGG$',sequence) or re.findall('^CAACGGTTTTGCC',sequence) or re.findall('CAACGGTTTTGCC$',sequence):
            print header_line
            print sequence
        header_line=line
        sequence = ''
        continue
    else:
        sequence += line
# below is needed to allow printing of final sequence which is not followed by a new fasta entry
if re.findall('^GTACAGTAGG',sequence) or re.findall('GTACAGTAGG$',sequence) or re.findall('^CAACGGTTTTGCC',sequence) or re.findall('CAACGGTTTTGCC$',sequence):
    print header_line
    print sequence
于 2012-12-15T15:32:43.030 回答
1

Python 有一个内置于字符串的方法，称为startswith()，还有一个称为endswith()。因此，您可以测试它是否以一个开头并以另一个结尾，反之亦然。

于 2012-12-14T21:07:12.897 回答