0

帮我，

我有一个fasta文件我想对它应用一些操作。我假设我的文件包含500个序列，for i=1 to 500我想取三个序列并应用一些功能，所以我会在166次里做同样的操作，每次我取3个序列

for i=1 to 500 do (500 number of sequences in fasta file)

take 3 sequences

    apply some functions

示例： 我的文件包含 9 个序列

   1-tatctattaccc

   2-gctgcgataagc

   3-tcctacttttgt

   4-caggaaaagaaa

   5-actgaatccctt

   6-ctgaagttgact

   7-aggtttgaagtg

   8-aacttccaactc

   9-gaaaagcaccct

我采用前 3 个序列

       seq1-tatctattaccc

       seq2-gctgcgataagc

       seq3-tcctacttttgt

我应用了一些函数，然后，我取序列号 4、5、6，我做与序列号 1、2、3 相同的事情，然后我对 7、8、9 做同样的事情 这是我的函数：

def identical(input_string):
code={"a":0,"c":1,"g":2,"t":3}
p=[code[i] for i in input_string]
n=len(input_string)
c=0
for i, n in enumerate(range(n, 0, -1)):
    c +=p[i]*(4**(n-1))
    return c+1

我的函数必须只使用这个：{"a","c","g","t"}，但在 fasta 文件序列中以'>'开头，如下所示：

>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG
AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG
CCGCCTCGGGAGCGTCCATGGCGGGTTTGAACCTCTAGCCCGGCGCAGTTTGGGCGCCAAGCCATATGAA
AGCATCACCGGCGAATGGCATTGTCTTCCCCAAAACCCGGAGCGGCGGCGTGCTGTCGCGTGCCCAATGA
ATTTTGATGACTCTCGCAAACGGGAATCTTGGCTCTTTGCATCGGATGGAAGGACGCAGCGAAATGCGAT
AAGTGGTGTGAATTGCAAGATCCCGTGAACCATCGAGTCTTTTGAACGCAAGTTGCGCCCGAGGCCATCA
GGCTAAGGGCACGCCTGCTTGGGCGTCGCGCTTCGTCTCTCTCCTGCCAATGCTTGCCCGGCATACAGCC
AGGCCGGCGTGGTGCGGATGTGAAAGATTGGCCCCTTGTGCCTAGGTGCGGCGGGTCCAAGAGCTGGTGT
TTTGATGGCCCGGAACCCGGCAAGAGGTGGACGGATGCTGGCAGCAGCTGCCGTGCGAATCCCCCATGTT
GTCGTGCTTGTCGGACAGGCAGGAGAACCCTTCCGAACCCCAATGGAGGGCGGTTGACCGCCATTCGGAT
GTGACCCCAGGTCAGGCGGGGGCACCCGCTGAGTTTACGC

>gi|2765657|emb|Z78532.1|CCZ78532 C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAGAATATATGATCGAGTG
AATCTGGAGGACCTGTGGTAACTCAGCTCGTCGTGGCACTGCTTTTGTCGTGACCCTGCTTTGTTGTTGG
GCCTCCTCAAGAGCTTTCATGGCAGGTTTGAACTTTAGTACGGTGCAGTTTGCGCCAAGTCATATAAAGC
ATCACTGATGAATGACATTATTGTCAGAAAAAATCAGAGGGGCAGTATGCTACTGAGCATGCCAGTGAAT
TTTTATGACTCTCGCAACGGATATCTTGGCTCTAACATCGATGAAGAACGCAGCTAAATGCGATAAGTGG
TGTGAATTGCAGAATCCCGTGAACCATCGAGTCTTTGAACGCAAGTTGCGCTCGAGGCCATCAGGCTAAG
GGCACGCCTGCCTGGGCGTCGTGTGTTGCGTCTCTCCTACCAATGCTTGCTTGGCATATCGCTAAGCTGG
CATTATACGGATGTGAATGATTGGCCCCTTGTGCCTAGGTGCGGTGGGTCTAAGGATTGTTGCTTTGATG
GGTAGGAATGTGGCACGAGGTGGAGAATGCTAACAGTCATAAGGCTGCTATTTGAATCCCCCATGTTGTT
GTATTTTTTCGAACCTACACAAGAACCTAATTGAACCCCAATGGAGCTAAAATAACCATTGGGCAGTTGA
TTTCCATTCAGATGCGACCCCAGGTCAGGCGGGGCCACCCGCTGAGTTGAGGC

所以，

if line=='>' then  pass

即：>gi|2765657|emb|Z78532.1|CCZ78532 C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA必须忽略像这样的行当我使用我的函数时，它给出了这个错误：

KeyError: '>'

或错误如：KeyError: 'CGT'

我该怎么做？

4

4 回答 4

1

在 python 2.x 上：

import itertools
with open('path/to/file') as infile:
    for seq1, seq2, seq3 in itertools.izip(infile, infile, infile):
        do_something(seq1, seq2, seq3)

我itertools.izip通常使用的原因zip是因为在 python 2.x 中，zip返回一个列表，然后您对其进行迭代。itertools.izip另一方面，迭代适当的元素而不首先创建列表。这节省了内存使用量，以及分配内存和创建列表所需的时间。

在 python 3.x 上，内置函数与 python 2.x 中的zip功能完全相同itertools.izip（要在 python 3.x 中获得 python 2.x 的zip功能，您需要这样做L = list(zip(…))）：

with open('path/to/file') as infile:
    for seq1, seq2, seq3 in zip(infile, infile, infile):
        do_something(seq1, seq2, seq3)

希望这可以帮助

于 2013-07-27T22:59:55.820 回答
1 
with open(file, 'r') as f:
    while True:
        line1 = f.readline()
        line2 = f.readline()
        line3 = f.readline()
        if not line3:
            break
        doSomething(line1, line2, line3)
于 2013-07-27T22:47:55.423 回答
1

你可以这样做：

import itertools

with open('/tmp/lines.txt', 'rU') as fasta:
    data=itertools.izip_longest(*[fasta]*3)

即使行数可能不是 3 的倍数，这也具有传递所有行的优势。最后一个元组的其余部分（最后一个元组的 0、1 或 2）用 填充None。你会希望something()预料到这一点。

或者，如果您想扩展为三个变量：

with open('/tmp/lines.txt', 'rU') as fasta:
    for a,b,c in itertools.izip_longest(*[fasta]*3):
        something(a,b,c)

请注意，如果文件中的行数不是 3 的整数倍，则 izip_longest 会完成元组，因此None您可能拥有c或b, cas None。如果不完整，只需测试或丢弃最后一个元组。

如果您知道如果不完整则不会使用最后一个元组，请使用以下命令：

with open('/tmp/lines.txt', 'rU') as fasta:
    for a,b,c in itertools.izip(*[fasta]*3):
        something(a,b,c)
于 2013-07-27T23:11:07.770 回答
0

我认为在 python 中处理 fasta 序列的更好方法总是使用 Biopython。Ii 一开始可能有点棘手，但值得努力。

首先你应该安装 biopython。只需在终端中输入：

sudo apt-get install python-biopython

现在你在 python 中拥有了一个令人难以置信的 fasta 解析器！

>>> from Bio import SeqIO

因此，您现在可以遍历记录（序列）并独立检索它们的 id、描述、序列、字母表、获取反向补码等。

with open("your_fasta_file.fasta", "r") as infh:
    parser = SeqIO.parse(infh, "fasta")

现在您有了 fasta 文件的“解析器”，这是一个可以在循环中使用的迭代器：

for sequence in parser:
    do stuff with sequence.id 
    do stuff with str(sequence.seq) # See below why str(seq)
    do stuff with sequence.description

对于像 ">gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA 基因和 ITS1 和 ITS2 的记录"

sequence.id --> 类似于 "gi|2765658|emb|Z78533.1|"

sequence.description --> 将是整个标题 ">gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA 基因和 ITS1 和 ITS2"

sequence.seq --> 将是序列。此时您必须小心，因为 sequence.seq 返回一个 Seq 对象： Seq('AGCTAGTAGCTG... ATGAC', Alphabet())

如果要将序列作为普通字符串，只需使用：

str(Seq_object)
于 2015-03-17T08:21:26.427 回答