0

此代码用于从 fasta 文件中提取和拆分序列

outfile=open('outf','w')
for line in open('input'):
      if line[0]==">":
         outfile.write('\n')
     else:
    outfile.write(line.strip())

   outfile.close()
 all_codons=[]
 for line in open('outf', 'r'):
     seq=line.strip()
     codons = [seq[i:i+3] for i in xrange(0, len(seq), 3) if len(seq[i:i+3])==3]
     all_codons.append(codons)

然后，从分割的序列中，我想取三个序列，它的长度是 9（9 个碱基）例如：

    CGTAACAAG 
    AATCCGGAG 
    CCGCCTCGG

我将第一个序列分成三个碱基的 3 个子序列，因此，从一个序列中我获得 3 个子序列，我对另外两个序列做同样的事情。

像这样：

    CGT    AAC     AAG 
    AAT    CCG     GAG 
    CCG    CCT     CGG

例子：

identical_segment('CGT')

我想将此函数应用于三个序列的每个子序列，然后在所有 fasta 文件上应用相同的东西。所以，目的是获取矩阵，例如我取第一个子序列'CGT'并应用函数 same_segment() ，它返回 28，其余 8 个子序列相同。所以我得到一个矩阵（3,3）：

28         2             3
4          23            35
23         4             27

我能做些什么？

4

2 回答 2

0

您的代码中存在一些问题。

首先，您只需要文件中的某些行，丢弃其他行，然后将所需的行输出到文件中。我不确定为什么需要最后一步。线的直接处理更有效。

def processLines(inputname):
    all_codons=[]
    for line in open(inputname):
        if line[0]==">":
            seq=line.strip()
            codons = [seq[i:i+3] for i in xrange(0, len(seq), 3) if
                      len(seq[i:i+3])==3]
        all_codons.append(codons)
    return all_codons

此外，每次调用 same_segment 都会生成一个字典，用作从 str 到 score 的映射。当呼叫数量增加时，它可能会变得昂贵。为避免这种情况，您可以尝试两种方法：

code={"a":0,"c":1,"g":2,"t":3} 
def identical_segment(input_string):
   .... # what you have written

或创建一个实例包含字典的类。

为了处理多个文件，请执行以下操作：

output = [processLines(filename) for filename in filenames]
# filenames is an iterable

或者如果您想将输入名称映射到输出：

outputDict = {filename: processLines(filename) for 
              filename in filenames}

毕竟，在每个输出上调用您的分析函数并将它们写入输出文件。

总结一下你应该在这篇文章中学到的东西：

  1. 输出文件可能不是最佳选择，因为文件 IO 很昂贵。如果将其写入某个文件，则意味着您必须再次将其读入，这会加倍昂贵。

  2. 不应一遍又一遍地创建同一个对象。校对您的代码以确保不会发生这种情况。

  3. 将你的主要任务分成几个小任务，然后为每个任务想一个简单直观的方法开始。在这个例子中，我们有 processfiles-> analysis-> output_result

  4. 理解是在 Python 中迭代事物的一种有用方法，而且它更具可读性。您可以搜索列表理解和字典理解以了解更多信息。

自己尝试一下。我很乐意在这里阅读您改进的代码。

于 2013-07-21T03:04:39.203 回答
0

尝试使用BioPython从您的 fasta 文件中提取核苷酸序列。使用这个包，

from Bio import AlignIO

for record in AlignIO.parse('filename.fasta', 'fasta'):
    print record.id, record.seq

# or store in a new file
seqs = []
for record in AlignIO.parse('filename.fasta', 'fasta'):
    seqs.append(record.seq + '\n')

with open(outfile, 'w') as out:
    out.writelines(seqs)
于 2013-09-23T16:37:30.697 回答