python - 大文本文件和另一个文本文件中相应行的条形码分割(python)

Question

我正在努力让自己在 python 方面变得更好。有一些工具可以做这些事情，但我想自己做有两个原因。

1. 学习一些更好的方法
2. 操作灵活

我有两个文本文件，大小完全相同，行数也相同。我需要检查文本的第 2、6(每次 +4)行，查看它的开头文本，检查它是否与某些预定义文本相似，如果是，则将该行与相应文件中的 4 block 一起写入，并写入相同的内容另一个相应文件中的行。 (对于那些听起来很熟悉的人，我正在尝试从 Illumina 配对的末端序列数据中分离条形码数据)。

我已经有一个工作代码，但问题是它需要几天才能完成。 100,000 行我花了大约 10 分钟，我有 2 亿行。

我把代码和我的想法一起贴在这里。好的，我有 100 个键，它们是 ATCCGG、ACCTGG 等。但是，如果我有一个不匹配，我想认为它是正确的，例如 DOG 可以有 AOG、BOG、DIG、DAG、DOF、DOH....

def makehamming2(text,dist):

    dicthamming=dict()
    rep=["A","T","C","G"]

    if dist==1:
        for i in range(len(text)):

            for j in range(len(rep)):
                chars=list(text)
                if rep[j]<>chars[i]:
                    chars[i]=rep[j]
                    word="".join(chars)
                    dicthamming[word]=text
    return dicthamming

我正在使用 dist=1。

我将此函数用于 100 个条形码，因此，我的字典中大约有 ~100*18 个项目。

count=0
eachline=1
writeflag=0
seqlen=int(seqlen)
cutlen=len(cutsite)
infile=open(inf, "r")
for line in infile:
        count+=1
        if eachline==1:
            writeflag=0
            header=line
            eachline=2
        elif eachline==2:
            eachline=3
            line=line.strip()
            if line[0:6] in searchdict.keys():

            barcode=searchdict[line[0:6]]

            towritefile=outfile+"/"+barcode+".fastq"


            seq=line[6:seqlen+6]
            qualstart=6
            writeflag=1
            seqeach[barcode]=seqeach.get(barcode,0)+1

    elif eachline==3:
        eachline=4
        third=line
    elif eachline==4:

        eachline=1
        line=line.strip()
        if writeflag==1:
            qualline=line[qualstart:qualstart+seqlen]
            addToBuffer=header+seq+"\n"+third+qualline+"\n"
            bufferdict[towritefile]=bufferdict.get(towritefile,"")+addToBuffer


            Fourlinesofpair=getfrompair(inf2,count, seqlen)


            bufferdictpair[towritefile[:-6:]+"_2.fastq"]=\
            bufferdictpair.get(towritefile[:-6:]+"_2.fastq","")+Fourlinesofpair

                if (count/4)%10000==0:
                    print "writing" , str((count/4))
                    for key, val in bufferdict.items():

                        writefile1=open(key,"a")
                        writefile1.write(val)
                        bufferdict=dict()


                    for key, val in bufferdictpair.items():


                        writefile1=open(key,"a")
                        writefile1.write(val)
                        bufferdictpair=dict()


                    end=(time.time()-start)/60.0
                    print "finished writing", str(end) , "minutes"


    print "writing" , str(count/4)                
    for key, val in bufferdict.items():


        writefile1=open(key,"a")
        writefile1.write(val)
        bufferdict=dict()
        writefile1.close()
    for key, val in bufferdictpair.items():

        writefile1=open(key,"a")
        writefile1.write(val)
        bufferdictpair=dict()
        writefile1.close()

    end=(time.time()-start)/60.0
    print "finished writing", str(end) , "minutes"

getfrompair 是一个函数，

def getfrompair(inf2, linenum, length):

    info=open(inf2,"r")
    content=""
    for count, line in enumerate(info):
        #print str(count)

    if count == linenum-4:
        content=line
    if count == linenum-3:
        content=content+line.strip()[:length]+"\n"
    if count == linenum-2:
        content=content+line
    if count == linenum-1:
        content=content+line.strip()[:length]+"\n"
        #print str(count), content



        return content

所以，我的主要问题是如何优化它。在大多数情况下，我会假设此代码在至少 8GB 内存和 >4 核处理器中运行。我可以使用多处理器吗？我在这里的另一个线程中使用了建议的缓冲区，因为这比在每行之后写入磁盘更快。

预先感谢您对我的教导。

编辑 1根据 Ignacio 的建议，我进行了分析，“getfrompair”函数占用了一半以上的运行时间？有没有更好的方法从文件中获取特定行而无需在某个时间逐行检查。

剖析结果来自分数(10000 行，而不是原始的 8 亿行)

     68719 function calls in 2.902 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       66    0.000    0.000    0.000    0.000 :0(append)
       32    0.003    0.000    0.003    0.000 :0(close)
     2199    0.007    0.000    0.007    0.000 :0(get)
        8    0.002    0.000    0.002    0.000 :0(items)
        3    0.000    0.000    0.000    0.000 :0(iteritems)
      750    0.001    0.000    0.001    0.000 :0(join)
     7193    0.349    0.000    0.349    0.000 :0(keys)
    39977    0.028    0.000    0.028    0.000 :0(len)
        1    0.000    0.000    0.000    0.000 :0(mkdir)
      767    0.045    0.000    0.045    0.000 :0(open)
      300    0.000    0.000    0.000    0.000 :0(range)
        1    0.005    0.005    0.005    0.005 :0(setprofile)
       96    0.000    0.000    0.000    0.000 :0(split)
        1    0.000    0.000    0.000    0.000 :0(startswith)
        1    0.000    0.000    0.000    0.000 :0(stat)
     6562    0.016    0.000    0.016    0.000 :0(strip)
        4    0.000    0.000    0.000    0.000 :0(time)
       48    0.000    0.000    0.000    0.000 :0(update)
       46    0.004    0.000    0.004    0.000 :0(write)
      733    1.735    0.002    1.776    0.002 RC14100~.PY:273(getfrompair)
        1    0.653    0.653    2.889    2.889 RC14100~.PY:31(split)
        1    0.000    0.000    0.000    0.000 RC14100~.PY:313(makehamming)
        1    0.000    0.000    0.005    0.005 RC14100~.PY:329(processbc2)
       48    0.003    0.000    0.005    0.000 RC14100~.PY:344(makehamming2)
        1    0.006    0.006    2.896    2.896 RC14100~.PY:4(<module>)
     4553    0.015    0.000    0.025    0.000 RC14100~.PY:74(<genexpr>)
     2659    0.014    0.000    0.023    0.000 RC14100~.PY:75(<genexpr>)
     2659    0.013    0.000    0.023    0.000 RC14100~.PY:76(<genexpr>)
        1    0.001    0.001    2.890    2.890 RC14100~.PY:8(main)
        1    0.000    0.000    0.000    0.000 cProfile.py:5(<module>)
        1    0.000    0.000    0.000    0.000 cProfile.py:66(Profile)
        1    0.000    0.000    0.000    0.000 genericpath.py:15(exists)
        1    0.000    0.000    0.000    0.000 ntpath.py:122(splitdrive)
        1    0.000    0.000    0.000    0.000 ntpath.py:164(split)
        1    0.000    0.000    0.000    0.000 os.py:136(makedirs)
        1    0.000    0.000    2.902    2.902 profile:0(<code object <module> at 000000000211A9B0, file "RC14100~.PY", line 4>)
        0    0.000             0.000          profile:0(profiler)



Process "Profile" terminated, ExitCode: 00000000

Answer 1

您的 getfrompair 函数使这成为一个经典的 O(n^2) 问题，因为您每次获得匹配项时都会通读第二个文件。相反，您想做的是同时从两个文件中读取，这样您只需要浏览一次。 izip是做到这一点的方法。

from itertools import izip

for line,line2 in izip(infile, infile2):