Mafft仅使用Python创建一个文件

nicole 发表于 Dev

妮可

所以我正在做一个项目，以对齐序列ID及其代码。给我一个条形码文件，其中包含DNA序列的标签，即TTAGG。然后有几个标签（ATTAC，ACCAT等）从序列文件中删除并带有seq ID。例子：

sequence file --> SEQ 01  TTAGGAACCCAAA
barcode file --> TTAGG

我想要的输出文件将删除条形码，并使用它来创建新的fasta格式文件。示例：testfile.TTAGG，该文件在打开时应具有

>SEQ01
AACCCAAA

这些文件有几个。我想获取创建并运行它们的每个文件mafft，但是运行脚本时，它只集中于的一个文件mafft。我上面提到的文件可以正常mafft运行，但是在运行时，它仅运行最后创建的文件。

这是我的脚本：

#!/usr/bin/python

import sys
import os

fname = sys.argv[1]
barcodefname = sys.argv[2]

barcodefile = open(barcodefname, "r")
for barcode in barcodefile:
        barcode = barcode.strip()
        outfname = "%s.%s" % (fname, barcode)
        outf = open(outfname, "w+")
        handle = open(fname, "r")
        mafftname = outfname + ".mafft"
        for line in handle:
                newline = line.split()
                seq = newline[0]
                brc = newline[1]
                potential_barcode = brc[:len(barcode)]
                if potential_barcode == barcode:
                        outseq = brc[len(barcode):]
                        barcodeseq = ">%s\n%s\n" % (seq,outseq)
                        outf.write(barcodeseq)
        handle.close()
        outf.close()
cmd = "mafft %s > %s" % (outfname, mafftname)
os.system(cmd)
barcodefile.close()

我希望这很清楚！请帮忙！我尝试更改缩进量，并在关闭文件时进行调整。在大多数情况下，它根本不会制作.mafft文件，有时它会制作文件，但不会放任何东西，但是大多数情况下，它仅适用于最后创建的文件。

示例：代码的开头将创建文件，例如-

testfile.ATTAC
testfile.AGGAC
testfile.TTAGG

然后在运行mafft时仅创建testfile.TTAGG.mafft（使用正确的输入）

我尝试关闭outf文件，然后再次打开它，它告诉我我正在强制它。我已更改为outf仅写文件，没有任何更改。

ODiogoSilva

mafft只对齐最后一个文件的原因是因为其执行不在循环中。

如代码所示，您在循环的每个迭代中都创建一个输入文件名变量（outfname），但此变量始终在下一个迭代中被覆盖。因此，当您的代码最终到达mafft执行命令时，outfname变量将包含循环的最后一个文件名。

要解决此问题，只需在循环内插入mafft执行命令：

#!/usr/bin/python

import sys
import os

fname = sys.argv[1]
barcodefname = sys.argv[2]

barcodefile = open(barcodefname, "r")
for barcode in barcodefile:
        barcode = barcode.strip()
        outfname = "%s.%s" % (fname, barcode)
        outf = open(outfname, "w+")
        handle = open(fname, "r")
        mafftname = outfname + ".mafft"
        for line in handle:
                newline = line.split()
                seq = newline[0]
                brc = newline[1]
                potential_barcode = brc[:len(barcode)]
                if potential_barcode == barcode:
                        outseq = brc[len(barcode):]
                        barcodeseq = ">%s\n%s\n" % (seq,outseq)
                        outf.write(barcodeseq)
        handle.close()
        outf.close()
        cmd = "mafft %s > %s" % (outfname, mafftname)
        os.system(cmd)

barcodefile.close()