2

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

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

我有两个文本文件,大小完全相同,行数相同。我需要检查第 2 行、第 6 行(每次 +4)文本,查看它的开头文本,检查它是否类似于某些预定义的文本,如果是,则将该行与相应文件中的 4 块一起写入,并写入相同另一个相应文件中的行。(对于那些听起来很熟悉的人,我正在尝试从 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
4

2 回答 2

1

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

from itertools import izip

for line,line2 in izip(infile, infile2):
于 2012-09-26T19:03:13.927 回答
0
>>> def dist(w1,w2):
...     return len(w1)-sum(map(lambda x:int(x[0]==x[1]),zip(w1,w2)))
...
>>> dist("DOG","FOG")
1
>>> dist("DOG","FOF")
2
>>> words = ["DOG","FOG","DAG","CAT","RAT","AOG","AAG"]
>>> print filter(lambda x:dist(target,x)<2,words)
['DOG', 'FOG', 'DAG', 'AOG']

然后做你想做的事

>>> import itertools
>>> my_alphabet = ["A","T","C","G"]
>>> target = "ATG"
>>> print filter(lambda x:dist(x,target)<2,itertools.permutations(my_alphabet,len(target)))
[('A', 'T', 'C'), ('A', 'T', 'G'), ('A', 'C', 'G'), ('C', 'T', 'G')]
于 2012-09-25T16:08:41.463 回答