annotate DPOGALAXY.py @ 41:ddd9e3b58525 draft default tip

Uploaded
author jose_duarte
date Fri, 02 Feb 2024 16:28:15 +0000
parents 1a450881f456
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
35
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
1 import pickle
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
2 from Bio import SeqIO
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
3 import os
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
4 import pandas as pd
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
5 import numpy as np
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
6 from local_ctd import CalculateCTD
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
7 from local_AAComposition import CalculateDipeptideComposition
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
8 import sys
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
9 from Bio.SeqUtils.ProtParam import ProteinAnalysis
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
10
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
11
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
12 class PDPOPrediction:
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
13
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
14 def __init__(self, folder='location', mdl='', seq_file='fasta_file.fasta', ttable=11):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
15 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
16 Initialize PhageDPO prediction.
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
17 :param folder: data path
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
18 :param mdl: ml model, in this case ANN or SVM
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
19 :param seq_file: fasta file
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
20 :param ttable: Translational table. By default, The Bacterial, Archaeal and Plant Plastid Code Table 11
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
21 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
22 self.records = []
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
23 self.data = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
24 self.df_output = None
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
25 self.seqfile = seq_file
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
26 self.__location__ = os.path.realpath(os.path.join(os.getcwd(), folder))
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
27
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
28 with open(os.path.join(self.__location__, mdl), 'rb') as m:
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
29 self.model0 = pickle.load(m)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
30 self.model = self.model0.named_steps['clf']
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
31 self.scaler = self.model0.named_steps['scl']
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
32 self.selectk = self.model0.named_steps['selector']
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
33 self.name = 'model'
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
34
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
35 for seq in SeqIO.parse(os.path.join(self.__location__, self.seqfile), 'fasta'):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
36 record = []
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
37 DNA_seq = seq.seq
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
38 AA_seq = DNA_seq.translate(table=ttable)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
39 descr_seq = seq.description.replace(' ', '')
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
40 self.data[descr_seq] = [DNA_seq._data, AA_seq._data]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
41 record.append(seq.description)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
42 record.append(DNA_seq._data)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
43 record.append(AA_seq._data)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
44 self.records.append(record)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
45
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
46 columns = ['ID', 'DNAseq', 'AAseq']
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
47 self.df = pd.DataFrame(self.records, columns=columns)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
48 #self.df = self.df.set_index('ID')
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
49 self.df.update(self.df.DNAseq[self.df.DNAseq.apply(type) == list].str[0])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
50 self.df.update(self.df.AAseq[self.df.AAseq.apply(type) == list].str[0])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
51
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
52 def Datastructure(self):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
53 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
54 Create dataset with all features
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
55 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
56 def count_orf(orf_seq):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
57 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
58 Function to count open reading frames
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
59 :param orf_seq: sequence to analyze
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
60 :return: dictionary with open reading frames
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
61 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
62 dic = {'DNA-A': 0, 'DNA-C': 0, 'DNA-T': 0, 'DNA-G': 0, 'DNA-GC': 0}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
63 for letter in range(len(orf_seq)):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
64 for k in range(0, 4):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
65 if str(orf_seq[letter]) in list(dic.keys())[k][-1]:
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
66 dic[list(dic.keys())[k]] += 1
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
67 dic['DNA-GC'] = ((dic['DNA-C'] + dic['DNA-G']) / (
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
68 dic['DNA-A'] + dic['DNA-C'] + dic['DNA-T'] + dic['DNA-G'])) * 100
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
69 return dic
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
70
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
71 def count_aa(aa_seq):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
72 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
73 Function to count amino acids
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
74 :param aa_seq: sequence to analyze
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
75 :return: dictionary with amino acid composition
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
76 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
77 dic = {'G': 0, 'A': 0, 'L': 0, 'V': 0, 'I': 0, 'P': 0, 'F': 0, 'S': 0, 'T': 0, 'C': 0,
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
78 'Y': 0, 'N': 0, 'Q': 0, 'D': 0, 'E': 0, 'R': 0, 'K': 0, 'H': 0, 'W': 0, 'M': 0}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
79 for letter in range(len(aa_seq)):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
80 if aa_seq[letter] in dic.keys():
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
81 dic[aa_seq[letter]] += 1
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
82 return dic
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
83
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
84 def sec_st_fr(aa_seq):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
85 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
86 Function to analyze secondary structure. Helix, Turn and Sheet
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
87 :param aa_seq: sequence to analyze
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
88 :return: dictionary with composition of each secondary structure
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
89 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
90 st_dic = {'Helix': 0, 'Turn': 0, 'Sheet': 0}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
91 stu = ProteinAnalysis(aa_seq).secondary_structure_fraction()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
92 st_dic['Helix'] = stu[0]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
93 st_dic['Turn'] = stu[1]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
94 st_dic['Sheet'] = stu[2]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
95 return st_dic
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
96
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
97
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
98 self.df_output = self.df.copy()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
99 self.df_output.drop(['DNAseq', 'AAseq'], axis=1, inplace=True)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
100 dna_feat = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
101 aa_len = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
102 aroma_dic = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
103 iso_dic = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
104 aa_content = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
105 st_dic_master = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
106 CTD_dic = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
107 dp = {}
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
108 self.df1 = self.df[['ID']].copy()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
109 self.df.drop(['ID'], axis=1, inplace=True)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
110 for i in range(len(self.df)):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
111 i_name = self.df.index[i]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
112 dna_feat[i] = count_orf(self.df.iloc[i]['DNAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
113 aa_len[i] = len(self.df.iloc[i]['AAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
114 aroma_dic[i] = ProteinAnalysis(self.df.iloc[i]['AAseq']).aromaticity()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
115 iso_dic[i] = ProteinAnalysis(self.df.iloc[i]['AAseq']).isoelectric_point()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
116 aa_content[i] = count_aa(self.df.iloc[i]['AAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
117 st_dic_master[i] = sec_st_fr(self.df.iloc[i]['AAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
118 CTD_dic[i] = CalculateCTD(self.df.iloc[i]['AAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
119 dp[i] = CalculateDipeptideComposition(self.df.iloc[i]['AAseq'])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
120 for j in self.df.index:
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
121 self.df.loc[j, dna_feat[j].keys()] = dna_feat[j].values() #dic with multiple values
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
122 self.df.loc[j, 'AA_Len'] = int(aa_len[j]) #dic with one value
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
123 self.df.loc[j, 'Aromaticity'] = aroma_dic[j]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
124 self.df.loc[j, 'IsoelectricPoint'] = iso_dic[j]
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
125 self.df.loc[j, aa_content[j].keys()] = aa_content[j].values()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
126 self.df.loc[j, st_dic_master[j].keys()] = st_dic_master[j].values()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
127 self.df.loc[j, CTD_dic[j].keys()] = CTD_dic[j].values()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
128 self.df.loc[j, dp[j].keys()] = dp[j].values()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
129 self.df.drop(['DNAseq', 'AAseq'], axis=1, inplace=True)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
130
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
131 def Prediction(self):
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
132 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
133 Predicts the percentage of each CDS being depolymerase.
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
134 :return: model prediction
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
135 """
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
136 scores = self.model0.predict_proba(self.df.iloc[:, :])
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
137 pos_scores = np.empty((self.df.shape[0], 0), float)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
138 for x in scores:
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
139 pos_scores = np.append(pos_scores, round(x[1]*100))
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
140 self.df_output.reset_index(inplace=True)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
141 self.df_output.rename(columns={'index': 'CDS'}, inplace=True)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
142 self.df_output['CDS'] += 1
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
143 self.df_output['{} DPO Prediction (%)'.format(self.name)] = pos_scores
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
144 self.df_output.to_html('output.html', index=False, justify='center')
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
145
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
146
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
147 if __name__ == '__main__':
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
148 __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
149
40
1a450881f456 Uploaded
jose_duarte
parents: 38
diff changeset
150 model = 'svm'
38
9558da071ec9 Uploaded
jose_duarte
parents: 35
diff changeset
151 fasta_file = sys.argv[1]
35
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
152
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
153 PDPO = PDPOPrediction(__location__, model, fasta_file)
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
154 PDPO.Datastructure()
a662eb3f87c2 Uploaded
jose_duarte
parents:
diff changeset
155 PDPO.Prediction()