|
32
|
1 import pickle
|
|
|
2 from Bio import SeqIO
|
|
|
3 import os
|
|
|
4 import pandas as pd
|
|
|
5 import numpy as np
|
|
|
6 from local_ctd import CalculateCTD
|
|
|
7 from local_AAComposition import CalculateDipeptideComposition
|
|
|
8 import sys
|
|
|
9 from Bio.SeqUtils.ProtParam import ProteinAnalysis
|
|
|
10
|
|
25
|
11
|
|
|
12 class PDPOPrediction:
|
|
32
|
13
|
|
|
14 def __init__(self, folder='location', mdl='', seq_file='fasta_file.fasta', ttable=11):
|
|
|
15 """
|
|
|
16 Initialize PhageDPO prediction.
|
|
|
17 :param folder: data path
|
|
|
18 :param mdl: ml model, in this case ANN or SVM
|
|
|
19 :param seq_file: fasta file
|
|
|
20 :param ttable: Translational table. By default, The Bacterial, Archaeal and Plant Plastid Code Table 11
|
|
|
21 """
|
|
|
22 self.records = []
|
|
25
|
23 self.data = {}
|
|
|
24 self.df_output = None
|
|
|
25 self.seqfile = seq_file
|
|
32
|
26 self.__location__ = os.path.realpath(os.path.join(os.getcwd(), folder))
|
|
25
|
27
|
|
32
|
28 with open(os.path.join(self.__location__, mdl), 'rb') as m:
|
|
25
|
29 self.model = pickle.load(m)
|
|
|
30 if mdl == 'SVM4311':
|
|
32
|
31 with open(os.path.join(__location__, 'd4311_SCALER'), 'rb') as sl:
|
|
25
|
32 self.scaler = pickle.load(sl)
|
|
|
33 self.name = mdl
|
|
|
34 elif mdl == 'ANN7185':
|
|
32
|
35 with open(os.path.join(__location__, 'd7185_SCALER'), 'rb') as sc:
|
|
25
|
36 self.scaler = pickle.load(sc)
|
|
|
37 self.name = mdl
|
|
|
38
|
|
32
|
39 for seq in SeqIO.parse(os.path.join(self.__location__, self.seqfile), 'fasta'):
|
|
|
40 record = []
|
|
25
|
41 DNA_seq = seq.seq
|
|
|
42 AA_seq = DNA_seq.translate(table=ttable)
|
|
32
|
43 descr_seq = seq.description.replace(' ', '')
|
|
|
44 self.data[descr_seq] = [DNA_seq._data, AA_seq._data]
|
|
|
45 record.append(seq.description)
|
|
|
46 record.append(DNA_seq._data)
|
|
|
47 record.append(AA_seq._data)
|
|
|
48 self.records.append(record)
|
|
|
49
|
|
|
50 columns = ['ID', 'DNAseq', 'AAseq']
|
|
|
51 self.df = pd.DataFrame(self.records, columns=columns)
|
|
|
52 #self.df = self.df.set_index('ID')
|
|
|
53 self.df.update(self.df.DNAseq[self.df.DNAseq.apply(type) == list].str[0])
|
|
|
54 self.df.update(self.df.AAseq[self.df.AAseq.apply(type) == list].str[0])
|
|
25
|
55
|
|
|
56 def Datastructure(self):
|
|
32
|
57 """
|
|
|
58 Create dataset with all features
|
|
|
59 """
|
|
25
|
60 def count_orf(orf_seq):
|
|
32
|
61 """
|
|
|
62 Function to count open reading frames
|
|
|
63 :param orf_seq: sequence to analyze
|
|
|
64 :return: dictionary with open reading frames
|
|
|
65 """
|
|
25
|
66 dic = {'DNA-A': 0, 'DNA-C': 0, 'DNA-T': 0, 'DNA-G': 0, 'DNA-GC': 0}
|
|
|
67 for letter in range(len(orf_seq)):
|
|
|
68 for k in range(0, 4):
|
|
32
|
69 if str(orf_seq[letter]) in list(dic.keys())[k][-1]:
|
|
25
|
70 dic[list(dic.keys())[k]] += 1
|
|
|
71 dic['DNA-GC'] = ((dic['DNA-C'] + dic['DNA-G']) / (
|
|
|
72 dic['DNA-A'] + dic['DNA-C'] + dic['DNA-T'] + dic['DNA-G'])) * 100
|
|
|
73 return dic
|
|
|
74
|
|
|
75 def count_aa(aa_seq):
|
|
32
|
76 """
|
|
|
77 Function to count amino acids
|
|
|
78 :param aa_seq: sequence to analyze
|
|
|
79 :return: dictionary with amino acid composition
|
|
|
80 """
|
|
25
|
81 dic = {'G': 0, 'A': 0, 'L': 0, 'V': 0, 'I': 0, 'P': 0, 'F': 0, 'S': 0, 'T': 0, 'C': 0,
|
|
|
82 'Y': 0, 'N': 0, 'Q': 0, 'D': 0, 'E': 0, 'R': 0, 'K': 0, 'H': 0, 'W': 0, 'M': 0}
|
|
|
83 for letter in range(len(aa_seq)):
|
|
|
84 if aa_seq[letter] in dic.keys():
|
|
|
85 dic[aa_seq[letter]] += 1
|
|
|
86 return dic
|
|
|
87
|
|
|
88 def sec_st_fr(aa_seq):
|
|
32
|
89 """
|
|
|
90 Function to analyze secondary structure. Helix, Turn and Sheet
|
|
|
91 :param aa_seq: sequence to analyze
|
|
|
92 :return: dictionary with composition of each secondary structure
|
|
|
93 """
|
|
25
|
94 st_dic = {'Helix': 0, 'Turn': 0, 'Sheet': 0}
|
|
|
95 stu = ProteinAnalysis(aa_seq).secondary_structure_fraction()
|
|
|
96 st_dic['Helix'] = stu[0]
|
|
|
97 st_dic['Turn'] = stu[1]
|
|
|
98 st_dic['Sheet'] = stu[2]
|
|
|
99 return st_dic
|
|
|
100
|
|
|
101 self.feat={"SVM4311": ["DNA-A", "DNA-T", "DNA-G", "DNA-GC", "AA_Len", "G", "A", "S", "T", "N", "Turn", "Sheet",
|
|
|
102 "_PolarizabilityC1", "_PolarizabilityC3", "_SolventAccessibilityC1", "_SecondaryStrC1",
|
|
|
103 "_SecondaryStrC2", "_SecondaryStrC3", "_ChargeC2", "_ChargeC3", "_PolarityC1", "_NormalizedVDWVC1",
|
|
|
104 "_NormalizedVDWVC3", "_HydrophobicityC2", "_HydrophobicityC3", "_SecondaryStrT23",
|
|
|
105 "_NormalizedVDWVT13", "_PolarizabilityD1001", "_SolventAccessibilityD1001", "_SolventAccessibilityD2001",
|
|
|
106 "_SolventAccessibilityD3001", "_SecondaryStrD1025", "_ChargeD1075","_ChargeD2001", "_ChargeD2025",
|
|
|
107 "_ChargeD3025", "_ChargeD3050", "_PolarityD1075", "_PolarityD3025","_NormalizedVDWVD1001",
|
|
|
108 "_NormalizedVDWVD3050", "_HydrophobicityD2001", "DG", "DT", "GD"],
|
|
|
109 "ANN7185": ["DNA-GC", "AA_Len", "Aromaticity", "IsoelectricPoint", "G", "A", "L", "V", "I", "P", "F",
|
|
|
110 "S", "T", "C", "Y", "N", "Q", "D", "E", "R", "K", "H", "W", "M", "Turn", "Sheet", "_PolarizabilityC1",
|
|
|
111 "_PolarizabilityC2", "_PolarizabilityC3", "_SolventAccessibilityC1", "_SolventAccessibilityC2",
|
|
|
112 "_SecondaryStrC1", "_SecondaryStrC3", "_ChargeC1", "_ChargeC2", "_ChargeC3", "_PolarityC2",
|
|
|
113 "_NormalizedVDWVC2", "_NormalizedVDWVC3", "_HydrophobicityC1", "_HydrophobicityC2", "_SecondaryStrT13",
|
|
|
114 "_SecondaryStrT23", "_ChargeT12", "_ChargeT13", "_HydrophobicityT12", "_PolarizabilityD1001",
|
|
|
115 "_PolarizabilityD1025", "_PolarizabilityD1050", "_PolarizabilityD2001", "_PolarizabilityD3025",
|
|
|
116 "_PolarizabilityD3050", "_PolarizabilityD3075", "_SolventAccessibilityD1050", "_SolventAccessibilityD2001",
|
|
|
117 "_SolventAccessibilityD2025", "_SolventAccessibilityD2050", "_SolventAccessibilityD3025",
|
|
|
118 "_SolventAccessibilityD3050", "_SolventAccessibilityD3100", "_SecondaryStrD1025", "_SecondaryStrD1050",
|
|
|
119 "_SecondaryStrD1075", "_SecondaryStrD2001", "_SecondaryStrD2050", "_SecondaryStrD2075", "_ChargeD1050",
|
|
|
120 "_ChargeD1075", "_ChargeD1100", "_ChargeD2025", "_ChargeD3025", "_ChargeD3050", "_PolarityD2050",
|
|
|
121 "_PolarityD3050", "_NormalizedVDWVD1001", "_NormalizedVDWVD1050", "_NormalizedVDWVD2001", "_NormalizedVDWVD2025",
|
|
|
122 "_HydrophobicityD3001", "_HydrophobicityD3075", "AD", "AW", "AY", "RC", "RT", "NA", "NE",
|
|
|
123 "NG", "NP", "DE", "DQ", "DG", "DT", "DY", "CG", "CL", "CY", "CV", "EN", "QA", "QR", "QE",
|
|
|
124 "QI", "GA", "GR", "GD", "GQ", "GG", "GH", "GL", "GF", "GP", "GT", "GY", "HA", "HC", "HI",
|
|
|
125 "HK", "HP", "IC", "IG", "IS", "IT", "IW", "LA", "LR", "LH", "LI", "LK", "LP", "KQ", "KH",
|
|
|
126 "KS", "KT", "MQ", "MG", "MI", "FA", "FR", "FS", "FY", "PC", "PE", "PG", "PH", "PM", "PF",
|
|
|
127 "PT", "SA", "SD", "SC", "SQ", "SW", "TA", "TC", "TM", "WL", "WV", "YE", "YG", "YH", "YI",
|
|
|
128 "YL", "YK", "YM", "YS"]}
|
|
|
129
|
|
|
130 self.df_output = self.df.copy()
|
|
32
|
131 self.df_output.drop(['DNAseq', 'AAseq'], axis=1, inplace=True)
|
|
25
|
132 dna_feat = {}
|
|
|
133 aa_len = {}
|
|
|
134 aroma_dic = {}
|
|
|
135 iso_dic = {}
|
|
|
136 aa_content = {}
|
|
|
137 st_dic_master = {}
|
|
|
138 CTD_dic = {}
|
|
|
139 dp = {}
|
|
32
|
140 self.df1 = self.df[['ID']].copy()
|
|
|
141 self.df.drop(['ID'], axis=1, inplace=True)
|
|
25
|
142 for i in range(len(self.df)):
|
|
|
143 i_name = self.df.index[i]
|
|
32
|
144 dna_feat[i] = count_orf(self.df.iloc[i]['DNAseq'])
|
|
|
145 aa_len[i] = len(self.df.iloc[i]['AAseq'])
|
|
|
146 aroma_dic[i] = ProteinAnalysis(self.df.iloc[i]['AAseq']).aromaticity()
|
|
|
147 iso_dic[i] = ProteinAnalysis(self.df.iloc[i]['AAseq']).isoelectric_point()
|
|
|
148 aa_content[i] = count_aa(self.df.iloc[i]['AAseq'])
|
|
|
149 st_dic_master[i] = sec_st_fr(self.df.iloc[i]['AAseq'])
|
|
|
150 CTD_dic[i] = CalculateCTD(self.df.iloc[i]['AAseq'])
|
|
|
151 dp[i] = CalculateDipeptideComposition(self.df.iloc[i]['AAseq'])
|
|
25
|
152 for j in self.df.index:
|
|
|
153 self.df.loc[j, dna_feat[j].keys()] = dna_feat[j].values() #dic with multiple values
|
|
|
154 self.df.loc[j, 'AA_Len'] = int(aa_len[j]) #dic with one value
|
|
|
155 self.df.loc[j, 'Aromaticity'] = aroma_dic[j]
|
|
|
156 self.df.loc[j, 'IsoelectricPoint'] = iso_dic[j]
|
|
|
157 self.df.loc[j, aa_content[j].keys()] = aa_content[j].values()
|
|
|
158 self.df.loc[j, st_dic_master[j].keys()] = st_dic_master[j].values()
|
|
|
159 self.df.loc[j, CTD_dic[j].keys()] = CTD_dic[j].values()
|
|
|
160 self.df.loc[j, dp[j].keys()] = dp[j].values()
|
|
32
|
161 self.df.drop(['DNAseq', 'AAseq'], axis=1, inplace=True)
|
|
25
|
162
|
|
|
163 def Prediction(self):
|
|
32
|
164 """
|
|
|
165 Predicts the percentage of each CDS being depolymerase.
|
|
|
166 :return: None
|
|
|
167 """
|
|
|
168 ft_scaler = pd.DataFrame(self.scaler.transform(self.df.iloc[:, :]), index=self.df.index, columns=self.df.columns)
|
|
25
|
169 ft_scaler = ft_scaler.drop(columns=[col for col in self.df if col not in self.feat[self.name]], axis=1)
|
|
|
170 scores = self.model.predict_proba(ft_scaler)
|
|
|
171 pos_scores = np.empty((self.df.shape[0], 0), float)
|
|
|
172 for x in scores:
|
|
|
173 pos_scores = np.append(pos_scores, round(x[1]*100))
|
|
|
174 self.df_output.reset_index(inplace=True)
|
|
32
|
175 print(self.df_output.columns)
|
|
|
176 self.df_output.rename(columns={'index': 'CDS'}, inplace=True)
|
|
|
177 self.df_output['CDS'] += 1
|
|
|
178 self.df_output['{} DPO Prediction (%)'.format(self.name)] = pos_scores
|
|
|
179 print(self.df_output)
|
|
25
|
180 #self.df_output = self.df_output.sort_values(by='{} DPO Prediction (%)'.format(self.name), ascending=False)
|
|
|
181 self.df_output.to_html('output.html', index=False, justify='center')
|
|
|
182
|
|
32
|
183
|
|
25
|
184 if __name__ == '__main__':
|
|
|
185 __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
|
|
|
186
|
|
|
187 model = sys.argv[1]
|
|
|
188 fasta_file = sys.argv[2]
|
|
|
189
|
|
32
|
190 PDPO = PDPOPrediction(__location__, model, fasta_file)
|
|
25
|
191 PDPO.Datastructure()
|
|
|
192 PDPO.Prediction()
|
|
|
193
|
