# HG changeset patch
# User martasampaio
# Date 1572280744 14400
# Node ID a555e95b2066b0e1b5bb4f55bd75d6eeeca505b3
# Parent  4289d1b76d668586a47f184486556cfcf4c16d48
Update phagepromoter.py -> pickle instead joblib

diff -r 4289d1b76d66 -r a555e95b2066 phagepromoter.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/phagepromoter.py	Mon Oct 28 12:39:04 2019 -0400
@@ -0,0 +1,570 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jul 19 13:45:05 2018
+
+@author: Marta
+"""
+
+from Bio import SeqIO
+import numpy as np
+import pandas as pd
+from auxiliar import free_energy,freq_base
+from Bio.Seq import Seq
+from Bio.SeqRecord import SeqRecord
+from Bio.Alphabet import IUPAC
+from auxiliar import get_bacteria, get_families, get_max_pssm, get_scores, get_lifecycle
+
+#division of the test genome in sequences of 65 bp
+def get_testseqs65(form,fic,both=False):
+    ALL = []
+    indexes = []
+    a = 0
+    rec = SeqIO.read(fic,form)
+    genome = rec.seq
+    i = 0
+    j = 65
+    while j < len(genome):
+        s = genome[i:j]
+        ALL.append([1,i,j,s])
+        i += 20
+        j += 20
+        a += 1
+        indexes.append(rec.name+":"+str(a))
+    if both:
+        comp = genome.reverse_complement()
+        size = len(rec.seq)
+        i = 0
+        j = 65
+        while j < len(comp):
+            s = comp[i:j]
+            ALL.append([-1,size-j,size-i,s])
+            i += 20
+            j += 20
+            a += 1
+            indexes.append(rec.name+":"+str(a))
+    df = pd.DataFrame(ALL, index=indexes, columns=['strand','iniprom','endprom','seq'])
+    return df
+
+#calculate the scores of all sequences (similar to get_posScores and get_negScores)
+def get_testScores(loc,test):
+    scores = []
+    posis = []
+    sizes = []
+    dic = {}
+    for ind,row in test.iterrows():
+        _,window = ind.split(':')
+        strand = row['strand']
+        ini = row['iniprom']
+        end = row['endprom']
+        seq = row['seq']
+        pos = [ini,end,strand]
+        dic[window] = pos
+        s = seq
+        score10_6,pos10_6 = get_scores(os.path.join(loc,'pssm10_6.txt'), s)
+        maxi10_6 = get_max_pssm(os.path.join(loc,'pssm10_6.txt'))
+        score10_8,pos10_8 = get_scores(os.path.join(loc,'pssm10_8.txt'), s)
+        maxi10_8 = get_max_pssm(os.path.join(loc,'pssm10_8.txt'))
+        scores23,pos23 = get_scores(os.path.join(loc,'pssm_23.txt'), s)
+        maxi23 = get_max_pssm(os.path.join(loc,'pssm_23.txt'))
+        scores21,pos21 = get_scores(os.path.join(loc,'pssm_21.txt'), s)
+        maxi21 = get_max_pssm(os.path.join(loc,'pssm_21.txt'))
+        scores27,pos27 = get_scores(os.path.join(loc,'pssm_27.txt'), s)
+        maxi27 = get_max_pssm(os.path.join(loc,'pssm_27.txt'))
+        scores32,pos32 = get_scores(os.path.join(loc,'pssm_32.txt'), s)
+        maxi32 = get_max_pssm(os.path.join(loc,'pssm_32.txt'))
+        score23 = max(scores23)
+        score21 = max(scores21)
+        score27 = max(scores27)
+        score32 = max(scores32)
+        maxiphage = max(score23,score21,score27,score32)
+        if maxiphage == score23: phage_max = score23*maxi23
+        elif maxiphage == score21: phage_max = score21*maxi21
+        elif maxiphage == score27: phage_max = score27*maxi27
+        elif maxiphage == score32: phage_max = score32*maxi32
+        score35_6,pos35_6 = get_scores(os.path.join(loc,'pssm35_6.txt'), s)
+        maxi35_6 = get_max_pssm(os.path.join(loc,'pssm35_6.txt'))
+        score35_9,pos35_9 = get_scores(os.path.join(loc,'pssm35_9.txt'), s)
+        maxi35_9 = get_max_pssm(os.path.join(loc,'pssm35_9.txt'))
+        score35_t4,pos35_t4 = get_scores(os.path.join(loc,'pssm35_t4.txt'), s)
+        maxi35_t4 = get_max_pssm(os.path.join(loc,'pssm35_t4.txt'))
+        score35_cbb,pos35_cbb = get_scores(os.path.join(loc,'pssm35_cbb.txt'), s)
+        maxi35_cbb = get_max_pssm(os.path.join(loc,'pssm35_cbb.txt'))
+        score35_lb,pos35_lb = get_scores(os.path.join(loc,'pssm35_lb.txt'),s)
+        maxi35_lb = get_max_pssm(os.path.join(loc,'pssm35_lb.txt'))
+        score35_mu, pos35_mu = get_scores(os.path.join(loc,'pssm35_mu.txt'),s)
+        maxi35_mu = get_max_pssm(os.path.join(loc,'pssm35_mu.txt'))
+        
+        dists6 = []
+        score6 = []
+        for p in pos10_6:
+            for a in range(14,22):
+                d = p-a-6
+                if d >= 0: 
+                    s10 = score10_6[p]
+                    s35_6 = score35_6[d]
+                    score6.append([s35_6,s10])
+                    dists6.append([d,p])
+        dists9 = []            
+        score9 = []
+        for p in pos10_6:
+            for a in range(11,14):
+                d = p-a-9
+                if d >= 0: 
+                    s10 = score10_6[p]
+                    s35_9 = score35_9[d]
+                    score9.append([s35_9,s10])
+                    dists9.append([d,p])
+        distst4 = []
+        scoret4 = []
+        distscbb = []
+        scorecbb = []
+        for p in pos10_6:
+            for a in range(16,18):
+                d = p-a-7
+                if d >= 0: 
+                    s10 = score10_6[p]
+                    s35_t4 = score35_t4[d]
+                    s35_cbb = score35_cbb[d]
+                    scoret4.append([s35_t4,s10])
+                    distst4.append([d,p])
+                    scorecbb.append([s35_cbb,s10])
+                    distscbb.append([d,p])
+                    
+        
+        distsmu = []
+        scoremu = []
+        for p in pos10_6:
+            d = p-16-14
+            if d >= 0: 
+                s10 = score10_6[p]
+                s35_mu = score35_mu[d]
+                scoremu.append([s35_mu,s10])
+                distsmu.append([d,p])
+        
+        distslb = []
+        scorelb = []
+        for p in pos10_6:
+            d = p-13-14
+            if d >= 0: 
+                s10 = score10_6[p]
+                s35_lb = score35_lb[d]
+                scorelb.append([s35_lb,s10])
+                distslb.append([d,p])
+    
+        soma6 = [sum(x) for x in score6]
+        soma9 = [sum(x) for x in score9]
+        somat4 = [sum(x) for x in scoret4]
+        somacbb = [sum(x) for x in scorecbb]
+        somamu = [sum(x) for x in scoremu]
+        somalb = [sum(x) for x in scorelb]
+        
+        maxi6 = max(soma6)
+        maxi9 = max(soma9)
+        maxit4 = max(somat4)
+        maxicbb = max(somacbb)
+        maximu = max(somamu)
+        maxilb = max(somalb)
+        maxi_elems = max(maxi6,maxi9,maxit4,maxicbb,maxilb,maximu)
+
+        if maxi_elems == maxilb:
+            indmax = somalb.index(maxilb)
+            sc35 = scorelb[indmax][0]
+            sc10 = scorelb[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = distslb[indmax]
+            size35 = 14
+            elems_maxi = sc35*maxi35_lb+sc10*maxi10_6
+        elif maxi_elems == maximu:
+            indmax = somamu.index(maximu)
+            sc35 = scoremu[indmax][0]
+            sc10 = scoremu[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = distsmu[indmax]
+            size35 = 14
+            elems_maxi = sc35*maxi35_mu+sc10*maxi10_6
+        elif maxi_elems == maxi9:
+            indmax = soma9.index(maxi9)
+            sc35 = score9[indmax][0]
+            sc10 = score9[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = dists9[indmax]
+            size35 = 9
+            elems_maxi = sc35*maxi35_9+sc10*maxi10_6
+        elif maxi_elems == maxit4:
+            indmax = somat4.index(maxit4)
+            sc35 = scoret4[indmax][0]
+            sc10 = scoret4[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = distst4[indmax]
+            size35 = 7
+            elems_maxi = sc35*maxi35_t4+sc10*maxi10_6
+        elif maxi_elems == maxicbb:
+            indmax = somacbb.index(maxicbb)
+            sc35 = scorecbb[indmax][0]
+            sc10 = scorecbb[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = distscbb[indmax]
+            size35 = 7
+            elems_maxi = sc35*maxi35_cbb+sc10*maxi10_6
+        else:
+            indmax = soma6.index(maxi6)
+            sc35 = score6[indmax][0]
+            sc10 = score6[indmax][1]
+            score_elems = [sc35,sc10]
+            posel = dists6[indmax]
+            size35 = 6
+            elems_maxi = sc35*maxi35_6+sc10*maxi10_6
+        
+        if score23 == maxiphage:
+            phage_score = score23
+            posiphage = scores23.index(score23)
+            sizephage = 23
+        elif score21 == maxiphage:
+            phage_score = score21
+            posiphage = scores21.index(score21)
+            sizephage = 21
+        elif score27 == maxiphage:
+            phage_score = score27
+            posiphage = scores27.index(score27)
+            sizephage = 27
+        else:
+            phage_score = score32
+            posiphage = scores32.index(score32)
+            sizephage = 32
+            
+        if elems_maxi >= max(score10_8)*maxi10_8:
+            i = posel[1]
+            ext = s[i-3:i-1]
+            if ext == 'TG': tg = 1
+            else: tg = 0
+            if elems_maxi > phage_max: host = 1
+            else: 
+                host = 0
+                tg = 0
+            sc = max(score10_8)
+            end35 = posel[0]+size35
+            dist = posel[1]-end35
+            scores.append([host, score_elems[1],sc,score_elems[0],phage_score,tg,dist,str(seq)])
+            posis.append([posel[1],posel[0],posiphage])
+            sizes.append([6,size35,sizephage])
+        else:
+            host = 1
+            sc = max(score10_8)
+            i = score10_8.index(sc)
+            ext = s[i-3:i-1]
+            if ext == 'TG': tg = 1
+            else: tg = 0
+            if max(score10_8)*maxi10_8 > phage_max: host = 1
+            else: 
+                host = 0
+                tg = 0
+            end35 = posel[0]+size35
+            dist = posel[1]-end35
+            scores.append([host,score_elems[1],sc,score_elems[0],phage_score,tg,dist,str(seq)])
+            posis.append([i,posel[0],posiphage])
+            sizes.append([8,size35,sizephage])
+    score = pd.DataFrame(scores, index=test.index, columns=['host','score10','score10_8','score35','score_phage','tg','dist','seq'])
+    posis = pd.DataFrame(posis, index=test.index, columns=['pos10','pos35','posphage'])
+    sizes = pd.DataFrame(sizes, index=test.index, columns=['size10','size35','size_phage'])
+    df_all = pd.concat([score,posis,sizes],axis=1)
+    return df_all,dic
+
+def create_dftest(scores_test,dic_window,family,bacteria,lifecycle):
+    tudo = []
+    tudo2 = []
+    for ind,row in scores_test.iterrows():
+        _,window = ind.split(':')
+        posis = dic_window[window]
+        strand=posis[2]
+        if strand == 1: ini=posis[0]
+        else: ini=posis[1]
+        seqprom = row['seq']
+        score10 = row['score10']
+        score10_8 = row['score10_8']
+        score35 = row['score35']
+        scorephage = row['score_phage']
+        size10 = row['size10']
+        size35 = row['size35']
+        sizephage = row['size_phage']
+        ini10 = row['pos10']
+        tg = row['tg']
+        host = row['host']
+        ini35 = row['pos35']
+        dist = row['dist']
+        end10=ini10+size10
+        iniphage = row['posphage']
+        endphage = iniphage+sizephage
+        if strand == 1:
+            if host == 0: 
+                new_seq = seqprom[iniphage:endphage]
+                new_ini = ini+iniphage+1
+                new_end = ini+endphage
+            else:
+                if size10 == 6: 
+                    new_seq = seqprom[ini35:end10]
+                    new_ini = ini+ini35+1
+                    new_end = ini+end10
+                else: 
+                    new_seq = seqprom[ini10:end10]
+                    new_ini = ini+ini10+1
+                    new_end = ini+end10
+            new_pos = '('+str(new_ini)+'..'+str(new_end)+')'
+        else:
+            if host == 0: 
+                new_seq = seqprom[iniphage:endphage]
+                new_ini = ini-endphage+1
+                new_end = ini-iniphage
+            else:
+                if size10 == 6: 
+                    new_seq = seqprom[ini35:end10]
+                    new_ini = ini-end10+1
+                    new_end = ini-ini35
+                else: 
+                    new_seq = seqprom[ini10:end10]
+                    new_ini = ini-end10+1
+                    new_end = ini-ini10
+            new_pos = 'complement('+str(new_ini)+'..'+str(new_end)+')'
+        if size10 == 6:
+            size10_6 = 1
+        else: 
+            size10_6 = 0
+        if size35 == 6:
+            size35_6 = 1
+            size35_7 = 0
+            size35_9 = 0
+            size35_14 = 0
+        elif size35 == 7:
+            size35_6 = 0
+            size35_7 = 1
+            size35_9 = 0
+            size35_14 = 0
+        elif size35 == 9:
+            size35_6 = 0
+            size35_7 = 0
+            size35_9 = 1
+            size35_14 = 0
+        else:
+            size35_6 = 0
+            size35_7 = 0
+            size35_9 = 0
+            size35_14 = 1
+        if sizephage == 23:
+            sizephage_23 = 1
+            sizephage_21 = 0
+            sizephage_32 = 0
+        elif sizephage == 21:
+            sizephage_23 = 0
+            sizephage_21 = 1
+            sizephage_32 = 0
+        elif sizephage == 32:
+            sizephage_23 = 0
+            sizephage_21 = 0
+            sizephage_32 = 1
+        else:
+            sizephage_23 = 0
+            sizephage_21 = 0
+            sizephage_32 = 0
+        if family == 'Podoviridae': 
+            Podo = 1
+            Sipho = 0
+            Myo = 0
+        elif family == 'Siphoviridae':
+            Podo = 0
+            Sipho = 1
+            Myo = 0
+        elif family == 'Myoviridae':
+            Podo = 0
+            Sipho = 0
+            Myo = 1
+        else: 
+            Podo = 0
+            Sipho = 0
+            Myo = 0
+        if bacteria == 'Bacillus':
+            bac = [1,0,0,0,0]
+        elif bacteria == 'Escherichia coli':
+            bac = [0,1,0,0,0]
+        elif bacteria == 'Klebsiella':
+            bac = [0,0,1,0,0]
+        elif bacteria == 'Pectobacterium':
+            bac = [0,0,0,1,0]
+        elif bacteria == 'Cronobacter':
+            bac = [0,0,0,0,1]
+        else:
+            bac = [0,0,0,0,0]
+        if lifecycle == 'virulent':  tp = 0
+        else: tp = 1
+        fe = free_energy(str(seqprom))
+        AT = freq_base(str(seqprom))
+        linha = [score10, score10_8, score35, dist, scorephage, fe, AT, host,size10_6,
+                 size35_6, size35_7, size35_9, size35_14,
+                 sizephage_23, sizephage_21,
+                 sizephage_32, tg, Podo, Sipho, Myo,tp]
+        linha.extend(bac)
+        tudo.append(linha)
+        linha2 = [new_pos,str(new_seq), host, size10_6,
+                 score10, score10_8, size35_6, size35_7, size35_9,size35_14, 
+                 score35, dist, sizephage_23, sizephage_21,
+                 sizephage_32, scorephage, tg, Podo, Sipho, Myo,tp, fe, AT]
+        linha2.extend(bac)
+        tudo2.append(linha2)
+    df_test = pd.DataFrame(tudo, index=scores_test.index, 
+                          columns = ['score10', 'score10_8','score35', 'dist35_10',
+                                     'scorephage','fe', 'freqAT',
+                                     'host','size10',
+                                     'size35_6', 'size35_7','size35_9','size35_14',
+                                     'sizephage_23', 
+                                     'sizephage_21', 'sizephage_32',
+                                     'TG', 'Podo', 'Sipho', 'Myo', 'tp',
+                                     'Bacillus', 'EColi', 
+                                     'Pectobacterium','Klebsiella', 'Cronobacter'])
+    df_INFO = pd.DataFrame(tudo2, index=scores_test.index, 
+                          columns = ['Positions','Promoter Sequence','host','size10',
+                                     'score10', 'score10_8',
+                                     'size35_6', 'size35_7','size35_9','size35_14',
+                                     'score35', 'dist35_10','sizephage_23', 
+                                     'sizephage_21', 'sizephage_32',
+                                     'scorephage', 'TG', 'Podo', 'Sipho', 'Myo', 'tp','fe', 'freqAT',
+                                     'EColi', 'Salmonella', 
+                                     'Pectobacterium','Cronobacter', 'Streptococcus'])
+    return df_test,df_INFO
+
+
+def get_predictions(scaler_file,model_file,test,df_testinfo,threshold):
+    import pickle
+    scaler = pickle.load(open(scaler_file, 'rb'))
+    model = pickle.load(open(model_file, 'rb'))
+    feat_scaled = pd.DataFrame(scaler.transform(test.iloc[:,:7]),index =test.index, columns=test.columns[:7])
+    TEST_scaled = pd.concat([feat_scaled,test.iloc[:,7:]],axis=1)
+    scores = model.predict_proba(TEST_scaled)
+    pos_scores = np.empty((TEST_scaled.shape[0],0), float)
+    for x in scores: pos_scores = np.append(pos_scores,x[1])
+    try: positive_indexes = np.nonzero(pos_scores>float(threshold))[0] #escolher os positivos, podia ser escolher com score > x
+    except ValueError: return 'The threshold value is not a float'
+    else:
+        if len(positive_indexes) == 0: return None
+        else:
+            positive_windows = TEST_scaled.index[positive_indexes]
+            INFO = df_testinfo.loc[positive_windows,['Positions','Promoter Sequence']]
+            promoter_type = []
+            for x in df_testinfo.loc[positive_windows,'host'].tolist():
+                if x == 0: promoter_type.append('phage')
+                else: promoter_type.append('host')
+            INFO['Type'] = promoter_type
+            INFO['Scores'] = np.around(pos_scores[positive_indexes],decimals=3)
+            INFO.index = positive_windows
+            return INFO
+
+
+def get_finaldf(test,rec):
+    new_df = test.groupby(['Positions'])
+    groups = list(new_df.groups.keys())
+    for i in range(len(groups)-1):
+        for j in range(i, len(groups)):
+            if 'complement' in groups[i]: inii = int(groups[i][11:].split('..')[0])
+            else: inii = int(groups[i][1:].split('..')[0])
+            if 'complement' in groups[j]: inij = int(groups[j][11:].split('..')[0])
+            else: inij = int(groups[j][1:].split('..')[0])
+            if inij < inii:
+                temp = groups[i]
+                groups[i] = groups[j]
+                groups[j] = temp
+    new_inds = []
+    for g in groups:
+        inds = new_df.groups[g]
+        if len(inds) == 1: new_inds.append(inds[0])
+        else:
+            maxi = max(new_df.get_group(g)['Scores'])
+            i = new_df.groups[g][new_df.get_group(g)['Scores']==maxi][0]
+            new_inds.append(i)
+    
+    output = test.loc[new_inds,:]
+    strands = []
+    new_pos = []
+    old_pos = output['Positions'].tolist()
+    
+    from Bio.SeqFeature import SeqFeature, FeatureLocation
+    feats = rec.features
+    for ind, row in output.iterrows():
+        pos = row['Positions']
+        if 'complement' in pos: 
+            strands.append('-')
+            new_pos.append(pos[10:])
+            ini,end= pos[11:-1].split('..')
+            new_loc = FeatureLocation(int(ini),int(end),strand=-1)
+        else: 
+            strands.append('+')
+            new_pos.append(pos)
+            ini,end= pos[1:-1].split('..')
+            new_loc = FeatureLocation(int(ini),int(end),strand=1)
+        feat = SeqFeature(new_loc, type='regulatory',qualifiers={'regulatory_class':['promoter'], 'note=':['predicted by PhagePromoter']})
+        feats.append(feat)            
+            
+    output.insert(loc=0, column='Strand', value=strands)
+    output['Positions'] = new_pos
+    output.to_html('output.html',index=False, justify='center')
+    recs = []
+    i = 0
+    for ind,row in output.iterrows():
+        s = Seq(row['Promoter Sequence'])
+        posis = old_pos[i]
+        typ = row['Type']
+        score = row['Scores']
+        sq = SeqRecord(seq=s, id=ind, description=typ+' '+posis+' score='+str(score))
+        recs.append(sq)
+        i += 1
+    SeqIO.write(recs, 'output.fasta','fasta')
+    new_rec = rec
+    new_rec.seq.alphabet = IUPAC.IUPACAmbiguousDNA()
+    new_feats = sorted(feats, key=lambda x: x.location.start)
+    new_rec.features = new_feats
+    SeqIO.write(new_rec,'output.gb','genbank')
+
+if __name__== "__main__":
+    
+    import sys
+    import os
+    __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
+    
+    gen_format = sys.argv[1]
+    genome_file = sys.argv[2]
+    both = sys.argv[3]
+    threshold = sys.argv[4]
+    family = sys.argv[5]
+    host = sys.argv[6]
+    phage_type = sys.argv[7]
+    model = sys.argv[8]
+    '''
+    gen_format = 'genbank'
+    genome_file = 'test-data/NC_015264.gb'
+    both = False
+    threshold = '0.50'
+    family = 'Podoviridae'
+    host = 'Pseudomonas'
+    phage_type = 'virulent'
+    model = 'SVM2400'
+    #model = 'ANN1600'
+    '''
+    
+    rec = SeqIO.read(genome_file, gen_format)
+    test_windows = get_testseqs65(gen_format, genome_file,both)
+    try: score_test,dic_window = get_testScores(__location__,test_windows)
+    except IndexError: print('Error. Input sequence can only have A,C,G or T')
+    else:
+        df_test,df_testinfo = create_dftest(score_test,dic_window,family,host,phage_type)
+        if model == 'ANN1600':
+            scaler_file = os.path.join(__location__, 'scaler1600.sav')
+            model_file = os.path.join(__location__, 'model1600.sav')
+            preds =  get_predictions(scaler_file, model_file, df_test,df_testinfo,threshold)
+            if preds is None: print('There is no sequence with a score value higher or equal to the threshold '+str(threshold))
+            elif type(preds) == str: print(preds)
+            else: output = get_finaldf(preds,rec)
+        else:
+            scaler_file = os.path.join(__location__, 'scaler2400.sav')
+            model_file = os.path.join(__location__, 'model2400.sav')
+            new_df_test = df_test.iloc[:,[0,1,2,3,4,5,6,7,8,9,13,14,16,17,19,20,22,24,25]]
+            preds =  get_predictions(scaler_file, model_file, new_df_test,df_testinfo,threshold)
+            if preds is None: print('There is no sequence with a score value higher or equal to the threshold '+str(threshold))
+            elif type(preds) == str: print(preds)
+            else: output = get_finaldf(preds,rec)
+