Mercurial > repos > martasampaio > phagepromoter
changeset 23:a455e43048cd draft
Uploaded
| author | martasampaio |
|---|---|
| date | Sat, 20 Apr 2019 11:09:34 -0400 |
| parents | 5acc4fa8b62d |
| children | 6459e7054353 |
| files | README.rst auxiliar.py model1600.sav model2400.sav phagepromoter.py phagepromoter.xml pssm10_6.txt pssm10_8.txt pssm35_6.txt pssm35_9.txt pssm35_cbb.txt pssm35_lb.txt pssm35_mu.txt pssm35_t4.txt pssm_21.txt pssm_23.txt pssm_27.txt pssm_32.txt scaler1600.sav scaler2400.sav tool_dependencies.xml |
| diffstat | 21 files changed, 894 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.rst Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,29 @@ +=============== +PhagePromoter +=============== + +Get promoter of phage genomes + +PhagePromoter is a python script that predicts promoter sequences in phage genomes, using a machine learning SVM model. This model was built from a train dataset with 19 features and 3200 examples (800 positives and 2400 negatives), each representing a 65 bp sequence of a phage genome. The positive cases represent the phage sequences that are already identified as promoters. + +**Inputs:** + +* genome format: fasta vs genbank; +* genome file: acepts both genbank and fasta formats; +* both strands (yes or no): allows the search in both DNA strands; +* threshold: represents the probability of the test sequence be a promoter (float between 0 and 1)" +* family: The family of the testing phage - Podoviridae, Siphoviridae or Myoviridae; +* Bacteria: The host of the phage. The train dataset include the following hosts: Bacillus, EColi, Salmonella, Pseudomonas, Yersinia, Klebsiella, Pectobacterium, Morganella, Cronobacter, Staphylococcus, Streptococcus, Streptomyces, Lactococcus. If the testing phage has a different host, select the option 'other', and it is recommended the use of a higher threshold value for more accurate results. +* phage type: The type of the phage, according to its lifecycle: virulent or temperate; + +**Outputs:** +This tool outputs two files: a FASTA file and a table in HTML, with the locations, sequence, score and type (recognized by host or phage RNAP) of the predicted promoters. + +**Requirements:** + +* Biopython +* Sklearn +* Numpy +* Pandas + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/auxiliar.py Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun May 27 17:37:09 2018 + +@author: Marta +""" + + +#get the phage host from the file 'bacteria.xlsx' +def get_bacteria(file): + import pandas as pd + df = pd.read_excel(file,header=0,index_col=0) + bacteria = {} + for ind,row in df.iterrows(): + bac = row['Bacteria'] + bacteria[ind] = bac + return bacteria + +#get the phage family from the file 'family.xlsx' +def get_families(file): + import pandas as pd + df = pd.read_excel(file,header=0,index_col=0) + families = {} + for ind,row in df.iterrows(): + fam = row['Family'] + families[ind] = fam + return families + +#get phage lifecycle from the file 'lifecycle.xlsx' +def get_lifecycle(file): + import pandas as pd + df = pd.read_excel(file,header=0,index_col=0) + types = {} + for ind,row in df.iterrows(): + lc = row['lifecycle'] + types[ind] = lc + return types + +#reads a file with a PSSM and return the max possible score of that PSSM +def get_max_pssm(file_pssm): + from Bio.Alphabet import IUPAC + from Bio.motifs import matrix + m = [] + fic = open(file_pssm,'r') + rf = fic.readline() + while rf: + new_l = [] + l = rf.strip().split('\t') + for val in l: + x = float(val) + new_l.append(x) + m.append(new_l) + rf = fic.readline() + a = IUPAC.unambiguous_dna + dic = {'A':m[0],'C':m[1], 'G':m[2], 'T':m[3]} + pssm = matrix.PositionSpecificScoringMatrix(a,dic) + return pssm.max + +#reads a file with a PSSM and returns a list of scores in all positions of the sequence +#returns the score divided by the maximum possible value +def get_scores(file_pssm, seq): + from Bio.Alphabet import IUPAC + from Bio.motifs import matrix + maxi = get_max_pssm(file_pssm) + m = [] + fic = open(file_pssm,'r') + rf = fic.readline() + while rf: + new_l = [] + l = rf.strip().split('\t') + for val in l: + x = float(val) + new_l.append(x) + m.append(new_l) + rf = fic.readline() + a = IUPAC.unambiguous_dna + dic = {'A':m[0],'C':m[1], 'G':m[2], 'T':m[3]} + pssm = matrix.PositionSpecificScoringMatrix(a,dic) + scores = [] + positions = [] + a = IUPAC.unambiguous_dna + seq.alphabet = a + for pos, score in pssm.search(seq, both=False,threshold=-50): + scores.append(score/maxi) + positions.append(pos) + return scores,positions + +#returns the frequencia of A and T bases in a sequence +def freq_base(seq): + A = seq.count('A') + T = seq.count('T') + AT = A+T + return AT + +#returns the free energy value of that sequence +def free_energy(seq): + dic1 = {'AA':-1.00, + 'TT':-1.00, + 'AT':-0.88, + 'TA':-0.58, + 'CA':-1.45, + 'AC':-1.44, + 'GG':-1.84, + 'CC':-1.84, + 'GA':-1.30, + 'AG':-1.28, + 'TC':-1.30, + 'CT':-1.28, + 'TG':-1.45, + 'GT':-1.44, + 'GC':-2.24, + 'CG':-2.17} + total = 0 + i = 0 + j = 1 + while i < len(seq)-1: + dint = seq[i]+seq[j] + total += dic1[dint] + i += 1 + j += 1 + return total \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/phagepromoter.py Sat Apr 20 11:09:34 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): + from sklearn.externals import joblib + scaler = joblib.load(scaler_file) + model = joblib.load(model_file) + 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) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/phagepromoter.xml Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,120 @@ +<tool id="get_proms" name="PhagePromoter" version="0.1.0"> + <description> +Get promoters of phage genomes + </description> + <requirements> + <requirement type="package">biopython</requirement> + <requirement type="package">scikit-learn</requirement> + <requirement type="package">numpy</requirement> + <requirement type="package">pandas</requirement> + </requirements> + <command detect_errors="exit_code" interpreter="python3"><![CDATA[ + phagepromoter.py "$input_type.genome_format" "$genome" "$both" "$threshold" "$family" "$bacteria" "$lifecycle" + "$adv.model" ]]> + </command> + <inputs> + <conditional name="input_type"> + <param type="select" name="genome_format" label='file format'> + <option value="genbank" selected="yes">genbank</option> + <option value="fasta">fasta</option> + </param> + <when value="genbank"> + <param type="data" name="genome" format="genbank" label='genome'/> + </when> + <when value="fasta"> + <param type="data" name="genome" format="fasta" label='genome'/> + </when> + </conditional> + <param type="boolean" name="both" label='Search both strands' checked="false" truevalue="-both" falsevalue="" /> + <param name="threshold" type="float" value="0.50" label="Threshold" help="Probabilty of being a promoter (float between 0 and 1)" /> + <param type="select" name="family" label='Phage family'> + <option value="Podoviridae" selected="yes">Podoviridae</option> + <option value="Siphoviridae">Siphoviridae</option> + <option value="Myoviridae">Myoviridae</option> + </param> + <param type="select" name="bacteria" label='Host bacteria Genus'> + <option value="Escherichia coli" selected="yes">Escherichia coli</option> + <option value="Salmonella">Salmonella</option> + <option value="Pseudomonas">Pseudomonas</option> + <option value="Yersinia">Yersinia</option> + <option value="Morganella">Morganella</option> + <option value="Cronobacter">Cronobacter</option> + <option value="Staphylococcus">Staphylococcus</option> + <option value="Streptococcus">Streptococcus</option> + <option value="Lactococcus">Lactococcus</option> + <option value="Streptomyces">Streptomyces</option> + <option value="Klebsiella">Klebsiella</option> + <option value="Bacillus">Bacillus</option> + <option value="Pectobacterium">Pectobacterium</option> + <option value="other">other</option> + </param> + <param type="select" name="lifecycle" label='Phage type'> + <option value="virulent" selected="yes">virulent</option> + <option value="temperate">temperate</option> + </param> + <section name = 'adv' title= 'Advanced Options' expanded = 'False'> + <param type = "select" name="model" label="Model"> + <option value="SVM2400" selected="yes">SVM2400</option> + <option value="ANN1600">ANN1600</option> + </param> + </section> + </inputs> + <outputs> + <data name="output1" format="html" from_work_dir="output.html" /> + <data name="output2" format="fasta" from_work_dir="output.fasta" /> + <data name="output3" format="genbank" from_work_dir="output.gb" /> + </outputs> + <tests> + <test> + <param name="genome_format" value="genbank"/> + <param name="genome" value="NC_015264.gb"/> + <param name="both" value="False"/> + <param name="threshold" value="0.50"/> + <param name="family" value="Podoviridae"/> + <param name="bacteria" value="Pseudomonas"/> + <param name="lifecycle" value="virulent"/> + <param name="model" value="SVM2400"/> + <output name="output1" file="output.html"/> + <output name="output2" file="output.fasta"/> + <output name="output3" file="output.gb"/> + </test> + </tests> + <help><![CDATA[ + +=============== +PhagePromoter +=============== + +Get promoters of phage genomes + +PhagePromoter is a python script that predicts promoter sequences in phage genomes, using machine learning models. Two different datasets were used to developed two models: the ANN model was built using a dataset with 26 features and 2400 examples (800 positives and 1600 negatives) and the SVM model was created using a dataset with 19 features and 3200 examples (800 positives and 2400 negatives). +Each example represents a sequence of 65 base pairs of a phage genome. The positive examples correspond to phage sequences already identified as promoters. + +**Inputs:** + +* genome format: fasta vs genbank (default); +* genome file: acepts both GenBank and FASTA formats; +* both strands: yes or no (default). Allows the search only in the direct strand or in both DNA strands; +* threshold: represents the probability of the test sequence being a promoter (a float between 0 and 1, default=0.50). For example, if threshold=0.90, the model will only return predicted sequences with more than 90% probability of being a promoter. The larger the genome, the higher the threshold should be. +* Family: The family of the testing phage - Podoviridae (default), Siphoviridae or Myoviridae; +* Host: The host of the phage. The training dataset include the following hosts: Bacillus, Escherichia coli (default), Salmonella, Pseudomonas, Yersinia, Klebsiella, Pectobacterium, Morganella, Cronobacter, Staphylococcus, Streptococcus, Streptomyces, Lactococcus. If the testing phage has a different host, select the option 'other'. +* Phage type: The type of the phage, according to its lifecycle: virulent or temperate; + +**Advanced options:** + +* Model: the user can choose which model to run: the SVM model (default) or the ANN model. The SVM model uses more negative data, so it will return less promoters but with a higher probability of being real promoters. However, it can fail to detect some of the real promoters. On the other hand, the ANN model will predict more promoters, so it can identify more real promoters, but it is expected to predict more false negatives. + +**Outputs:** + +This tool outputs two files: a FASTA file and a table in HTML, with the locations, sequence, score and type (recognized by host or phage RNAP) of the predicted promoters. +In addition, the tool will output a GenBank file with the predicted promoters as features. + +**Requirements:** + +* Biopython +* Sklearn +* Numpy +* Pandas + + ]]> </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm10_6.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-3.24 1.93 -0.34 1.38 1.43 -3.05 +-2.14 -4.24 -1.03 -1.44 -1.19 -4.05 +-2.29 -4.46 -1.12 -1.44 -1.53 -3.59 +1.79 -3.59 1.17 -0.61 -0.96 1.9
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm10_8.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +0.09 1.72 -4.49 1.92 1.95 1.95 -4.49 1.92 +-1.68 -1.32 -4.49 -4.49 -4.49 -4.49 -3.49 -4.49 +-0.79 -2.49 -4.49 -2.91 -4.49 -4.49 -2.49 -4.49 +1.03 -2.91 1.95 -4.49 -4.49 -4.49 1.88 -2.91
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_6.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-2.84 -2.12 -2.4 1.49 -0.74 1.16 +-2.95 -4.65 -2.56 -0.26 1.13 -1.56 +-4.33 -2.48 1.73 -3.33 -1.69 -1.65 +1.88 1.83 -1.65 -1.95 -0.14 0.15
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_9.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +0.93 0.79 -0.65 -1.87 -0.14 -1.46 -1.14 -0.14 1.13 +-1.46 -1.87 -1.46 -2.46 1.35 -3.46 -2.46 -0.14 -0.65 +-0.87 -1.87 -1.87 1.79 -1.87 -1.87 -3.46 -1.46 -1.87 +0.24 0.79 1.45 -3.46 -1.87 1.71 1.71 0.86 -0.14
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_cbb.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-2.46 -2.46 1.45 1.79 1.79 -2.46 0.13 +-2.46 -2.46 -2.46 -2.46 -2.46 1.54 -2.46 +-2.46 1.79 -0.14 -2.46 -2.46 -2.46 1.24 +1.79 -2.46 -2.46 -2.46 -2.46 -0.46 -1.46
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_lb.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-0.81 -0.81 -0.81 -0.81 -0.81 0.19 -0.81 -0.81 0.78 0.19 -0.81 -0.81 -0.81 -0.81 +-0.81 -0.81 -0.81 1.19 -0.81 -0.81 -0.81 0.19 -0.81 0.19 -0.81 -0.81 -0.81 1.19 +-0.81 -0.81 1.19 -0.81 1.19 -0.81 0.78 -0.81 0.19 -0.81 -0.81 -0.81 1.19 -0.81 +1.19 1.19 -0.81 -0.81 -0.81 0.78 0.19 0.78 -0.81 0.19 1.19 1.19 -0.81 -0.81
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_mu.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-1.17 -0.17 1.15 -1.17 1.15 1.15 -1.17 -0.17 0.42 -1.17 0.42 -1.17 0.83 -1.17 +0.83 0.83 -1.17 -0.17 -1.17 -1.17 0.83 1.15 0.83 1.42 -1.17 0.42 -1.17 0.83 +-0.17 -1.17 -0.17 -1.17 -1.17 -0.17 -1.17 -1.17 -1.17 -1.17 0.42 0.42 -0.17 -1.17 +-0.17 -0.17 -1.17 1.15 -0.17 -1.17 0.42 -1.17 -1.17 -1.17 -0.17 -0.17 -0.17 0.42
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm35_t4.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-2.43 -3.43 -2.43 -2.43 1.86 -2.43 1.33 +-3.43 -3.43 -3.43 -3.43 -3.43 1.86 -3.43 +1.82 -3.43 -2.43 -3.43 -3.43 -3.43 -3.43 +-2.43 1.9 1.82 1.86 -2.43 -3.43 0.38
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm_21.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-2.0 -1.0 -0.42 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 1.7 -2.0 -1.0 -1.0 -2.0 -2.0 -2.0 -2.0 0.81 0.58 +1.17 1.46 -2.0 -2.0 -1.0 -2.0 -1.0 -2.0 1.7 1.7 -0.42 -2.0 1.7 -2.0 0.81 1.7 1.7 -2.0 -2.0 0.81 0.81 +-0.42 -2.0 1.46 1.7 1.58 -2.0 -2.0 1.7 -2.0 -2.0 1.46 -2.0 -2.0 -1.0 -0.42 -2.0 -2.0 -2.0 -2.0 -2.0 -1.0 +-0.42 -1.0 -2.0 -2.0 -2.0 1.7 1.58 -2.0 -2.0 -2.0 -2.0 -2.0 -2.0 1.46 0.0 -2.0 -2.0 1.7 1.7 -2.0 -2.0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm_23.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +0.07 1.71 1.18 -0.18 1.65 0.6 -1.86 0.82 -2.86 -0.79 -5.18 1.88 -1.86 -4.18 1.69 0.97 1.26 -0.05 -0.79 0.07 1.11 0.28 1.25 +-0.32 -2.86 -2.86 -4.18 -2.86 0.49 0.6 0.52 1.68 -2.37 1.88 -1.86 1.84 -1.09 -2.86 -3.18 -1.18 -5.18 -0.54 -1.86 -1.86 -1.48 -1.09 +-3.59 -2.86 -3.59 -4.18 -3.18 -0.32 0.95 -0.86 -1.37 -1.59 -1.86 -5.18 -3.18 -5.18 -2.59 -3.59 -1.18 1.56 1.41 1.3 0.6 1.16 -0.18 +1.05 -1.18 0.6 1.59 -0.72 -1.86 -1.86 -2.01 -1.86 1.53 -5.18 -5.18 -5.18 1.79 -1.09 0.89 -0.48 -4.18 -3.59 -2.18 -4.18 -2.37 -1.86
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm_27.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +-3.09 -2.09 1.67 -1.5 -1.09 -0.28 -3.09 -3.09 -2.09 -3.09 -3.09 -3.09 1.87 -3.09 -3.09 -3.09 -3.09 -3.09 -1.5 -3.09 -3.09 -2.09 -3.09 1.72 1.56 1.5 -1.5 +1.82 -3.09 -3.09 1.67 1.56 1.3 -1.09 -0.5 1.16 1.67 1.67 -1.09 -3.09 1.87 -3.09 1.87 1.67 -0.77 -3.09 0.82 1.16 -0.28 0.91 -2.09 -2.09 -0.5 0.08 +-3.09 1.82 -2.09 -2.09 -3.09 -0.77 0.08 1.16 0.5 -3.09 -0.77 -2.09 -3.09 -3.09 -2.09 -3.09 -1.09 1.67 1.77 1.0 0.61 -3.09 -3.09 -1.5 -3.09 -3.09 -0.5 +-2.09 -3.09 -1.09 -2.09 -1.09 -3.09 1.23 -0.09 -3.09 -0.77 -3.09 1.67 -3.09 -3.09 1.82 -3.09 -2.09 -3.09 -3.09 -3.09 -3.09 1.5 0.91 -3.09 -0.5 -1.5 0.91
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pssm_32.txt Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,4 @@ +1.65 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 1.65 -1.81 -1.81 0.78 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 -1.81 -1.81 -1.81 1.65 1.65 +-1.81 -1.81 1.65 1.65 -1.81 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 -1.81 1.65 1.65 -1.81 -1.81 1.19 -1.81 -1.81 -1.81 -1.81 0.78 -1.81 1.65 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 +-1.81 1.65 -1.81 -1.81 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 0.78 -1.81 -1.81 -1.81 -0.81 -1.81 -1.81 1.51 1.65 1.65 1.65 -1.81 1.65 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 1.65 -1.81 -1.81 +-1.81 -1.81 -1.81 -1.81 1.65 -1.81 1.65 -1.81 -1.81 -1.81 -1.81 1.65 -1.81 -1.81 1.51 -1.81 0.19 -0.81 -1.81 -1.81 -1.81 0.78 -1.81 -1.81 1.65 -1.81 1.65 -1.81 1.65 -1.81 -1.81 -1.81
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Sat Apr 20 11:09:34 2019 -0400 @@ -0,0 +1,6 @@ +<tool_dependency> + <package name="biopython"></package> + <package name="numpy" ></package> + <package name="pandas"></package> + <package name="scikit-learn"></package> +</tool_dependency>