将条形码与序列python匹配？答案

【问题标题】：Matching Barcodes to sequences python?将条形码与序列python匹配？
【发布时间】：2015-04-26 04:20:56
【问题描述】：

我有序列文件和条形码文件。条形码文件可能具有任何长度的条形码，例如“ATTG、AGCT、ACGT”。例如，序列文件看起来像“ATTGCCCCCCCGGGGG, ATTGTTTTTTT, AGCTAAAA”。我需要将条形码与开头包含它们的序列进行匹配。然后对于具有相同条形码的每组序列，我必须使用程序的其余部分（已经编写）对它们进行计算。我只是不知道如何让他们匹配。我已经使用了打印语句，它被搞砸的部分是“potential_barcode = line(:len(barcode)”行。此外，它说#simple to fasta 那是我应该在匹配序列中阅读的地方。我对此很陌生，所以我可能犯了很多错误。感谢您的帮助！

bcodefname = sys.argv[1]
infname = sys.argv[2]
barcodefile = open(bcodefname, "r")
for barcode in barcodefile:
        barcode = barcode.strip()
        print "barcode: %s" % barcode
        outfname = "%s.%s" % (bcodefname,barcode)
#           print outfname
        outf = open("outfname", "w")
        handle = open(infname, "r")
        for line in handle:
                potential_barcode = line[:len(barcode)]
                print potential_barcode
                if potential_barcode == barcode:
                        outseq = line[len(barcode):]
                        sys.stdout.write(outseq)
                        outf.write(outseq)
                        fastafname = infname + ".fasta"
                        print fastafname
                        mafftfname = fastafname + ".mafft"
                        stfname = mafftfname + ".stock"
                        print stfname
#simp to fasta#
#                       handle = open(infname, "r")
                        outf2 = open(fastafname, "w")
                        for line in handle:
                                linearr = line.split()
                                seqid = linearr[0]
                                seq = linearr[1]
                                outf2.write(">%s\n%s\n" % (seqid,seq))
#                       handle.close()
#                       outf.close()
#mafft#
                        cmd = "mafft %s > %s" % (fastafname,mafftfname)
                        sys.stderr.write("command: %s\n" % cmd)
                        os.system(cmd)
                        sys.stderr.write("command done\n")

【问题讨论】：

我没有看到potential_barcode = line[:len(barcode)] 行有问题，它适用于我所有的简单测试用例。您能否提供一些产生错误结果的示例条形码/行对？
贾斯汀，由于您是新用户并且尚未发布更新，我会提一下，如果下面的答案对您有所帮助，您应该考虑将其标记为“已接受”。

标签： python regex sequence

【解决方案1】：

我不知道为什么这段代码不适合你。但这里有一些改进代码的技巧：

您正在读取每个条形码的整个序列文件。如果有 100 个条形码，则您将读取序列文件 100 次。你应该做的是读取一次条形码文件并创建一个条形码列表。

首先定义一个我们将用来检查匹配的函数：

def matches_barcode(sequence, barcodes):
  for bc in barcodes:
    if sequence.startswith(bc):
      return True
  return False

（请注意，我使用startswith，而不是构造一个新字符串并进行比较；startswith 应该更快。）

现在从文件中读取条形码：

barcodes = []
with open(bcodefname, "r") as barcodefile:
  for barcode in barcodefile:
    barcodes.append(barcode.strip())

（请注意，我使用了with open...；您的代码会泄漏打开的文件，如果您有很多条形码，这可能会阻止程序运行。）

然后将序列文件读一遍，检查每个序列是否与条形码匹配：

with open(infname, "r") as seqfile:
  for sequence in seqfile:
    if matches_barcode(sequence.strip(), barcodes):
      # Got a match, continue processing.

这会快很多，因为它执行的 I/O 少得多：它读取 2 个文件而不是 N + 1 个文件，其中 N 是条形码的数量。不过，它仍然是一个非常幼稚的算法：如果它太慢，您将不得不研究更复杂的算法来检查匹配。

如果您仍然没有得到您期望的匹配，您需要调试：打印出正在比较的确切字符串，这样您就可以看到发生了什么。在这些情况下使用repr 来真正查看数据是个好主意，包括空格和所有内容。

【讨论】：