|
0
|
1 #!/usr/bin/env python
|
|
|
2 # title :fibronectin.py
|
|
|
3 # description :This script will analyze fasta files, look for restriction sites, cut the sequences around the restriction sites, translate the nucleic sequences into amino acids sequences.
|
|
|
4 # author :Fabienne Wong Jun Tai
|
|
|
5 # date :20121107
|
|
|
6 # version :1.0
|
|
|
7 # usage :python fibronectin.py -i file.fasta -o /output/dir/ -p pattern -5 seq_restric_5'-3 seq_restric_3'
|
|
|
8 # notes :
|
|
|
9 # python_version :3.7.11
|
|
|
10 # biopython_max_version :1.72
|
|
|
11 # ==============================================================================
|
|
|
12 import math
|
|
|
13 import matplotlib
|
|
|
14 import numpy
|
|
|
15 import re
|
|
|
16 import subprocess
|
|
|
17 import matplotlib.pyplot as plot
|
|
|
18 from args import *
|
|
|
19 from Bio import SeqIO, Seq
|
|
|
20 from Bio.SubsMat import MatrixInfo as matlist
|
|
|
21 from Bio import pairwise2
|
|
|
22 from Bio.pairwise2 import format_alignment
|
|
|
23 matplotlib.use('Agg')
|
|
|
24
|
|
|
25 args = Args()
|
|
|
26 print(sys.path[0])
|
|
|
27 # Variables initialization
|
|
|
28 fibronectin_script_dir = sys.path[0]
|
|
|
29 print(fibronectin_script_dir)
|
|
|
30 directory = args.output_dir
|
|
|
31 mcl_file = directory + "mcl.in"
|
|
|
32 mcl_output = directory + "mcl.out"
|
|
|
33 html_file = directory + "fibronectin_report.html"
|
|
|
34 graph_pic = directory + "distri.png"
|
|
|
35 input_file = os.path.basename(args.input)
|
|
|
36 site_res_5 = args.site_res_5
|
|
|
37 site_res_3 = args.site_res_3
|
|
|
38 tag = {'mut': [], 'ok_stop_ext': [], 'stop': [], 'no_restric': [], 'no_multiple': [], 'amber': []}
|
|
|
39 all_seq = []
|
|
|
40 all_seq_fasta = {} # dictionnary that will store information about all the sequences
|
|
|
41 good_seq = {} # dictionnary that will store information about the valid sequences
|
|
|
42 identical_clones = {}
|
|
|
43 var_seq_common = {} # dictionnary that will store the number of sequences that share the same variable parts
|
|
|
44 align_scores = []
|
|
|
45 nb_var_part = 0
|
|
|
46
|
|
|
47
|
|
|
48 def reverse_complement(seq):
|
|
|
49 # Generate the reverse complement
|
|
|
50 complement_nuc = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
|
|
|
51 rev_com = ""
|
|
|
52 for n in (seq[::-1]):
|
|
|
53 rev_com += complement_nuc[n]
|
|
|
54 return rev_com
|
|
|
55
|
|
|
56
|
|
|
57 def generate_aln(seq_dic, ids):
|
|
|
58 # Multiple Sequence Alignment via ClustalO
|
|
|
59 input = ''
|
|
|
60 for k in ids:
|
|
|
61 input += '>%s\n%s\n' % (k, re.sub("(.{80})", "\\1\n", seq_dic[k]['prot'], re.DOTALL))
|
|
|
62 p = subprocess.Popen("clustalo -i - --outfmt clu", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, stdin=subprocess.PIPE, universal_newlines=True)
|
|
|
63 aln_out, aln_err = p.communicate(input=input)
|
|
|
64 print(type(aln_out))
|
|
|
65 return aln_out
|
|
|
66
|
|
|
67
|
|
|
68 def report_html(html_file, tag, all_seq, good_seq, all_seq_fasta, identical_clones, nb_var_part, var_seq_common, align_scores, args):
|
|
|
69 # Generate the html file for the report
|
|
|
70 all_seq.sort()
|
|
|
71 no_restric = tag['no_restric']
|
|
|
72 no_restric.sort()
|
|
|
73 no_multiple = tag['no_multiple']
|
|
|
74 no_multiple.sort()
|
|
|
75 stop = tag['stop']
|
|
|
76 stop.sort()
|
|
|
77 amber = tag['amber']
|
|
|
78 amber.sort()
|
|
|
79 mut = tag['mut']
|
|
|
80 mut.sort()
|
|
|
81 # good_ids = good_seq.keys()
|
|
|
82
|
|
|
83 good_seq = dict(sorted(good_seq.items()))
|
|
|
84 good_ids = good_seq.keys()
|
|
|
85
|
|
|
86 # good_ids.sort()
|
|
|
87
|
|
|
88 w = open(html_file, 'w')
|
|
|
89 w.write(
|
|
|
90 '<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN""http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8" /><title>Fibronectin Report</title><link href="http://twitter.github.com/bootstrap/assets/css/bootstrap.css" rel="stylesheet" /><style type="text/css">body {padding-top: 40px;}.subhead {padding: 40px 0;}.subhead h1 {font-size: 60px;}.fasta { font-family: Monaco, Menlo, Consolas, "Courier New", monospace; font-size: 12px;}code.grey{color: #636D71;}</style></head><body><a id="top"></a><div class="navbar navbar-fixed-top"><div class="navbar-inner"><div class="container"><a class="brand" href="#top">Fibronectin Report</a><div class="nav-collapse collapse"><ul class="nav"><li><a href="#input">Input data</a></li><li><a href="#analysis">Sequences analysis</a></li><li><a href="#variable">Variable regions analysis</a></li><li><a href="#cluster">Clustering</a></li><li><a href="#stat">Statistics</a></li><li><a href="#annex">Annex</a></li></ul></div></div></div></div><div class="container-fluid"><header class="subhead"><h1>Fibronectin Report</h1></header><div class="page-header"><a id="input"></a><h2>Input data</h2></div>')
|
|
|
91 # Input data
|
|
|
92 w.write(
|
|
|
93 '<p>Input file:<br/><code class="grey">%s</code></p><p>Number of sequences in input file:<br/><code class="grey">%d</code></p><p>Pattern of the sequence bank:<br/><code class="grey">%s</code></p><p>5\' restriction site:<br/><code class="grey">%s</code></p><p>3\' restriction site:<br/><code class="grey">%s</code></p>' % (
|
|
|
94 input_file, len(all_seq), args.pattern, args.site_res_5, args.site_res_3))
|
|
|
95 # Sequence analysis
|
|
|
96 w.write(
|
|
|
97 '<div class="page-header"><a id="analysis"></a><h2>Sequences analysis</h2></div><p>Caption:</p><ul><li class="text-success">Valid sequences that will be part of the next analysis </li><li class="text-warning">Good sequences but will not be part of the next analysis</li><li class="text-error">Rejected sequences</li></ul><table class="table table-striped table-bordered"><tr><th class="text-error">Absence of restriction sites</th><th class="text-error">Incorrect number of nucleotides between the restriction sites</th><th class="text-error">Stop codon <u>inside</u> the area of interest</th><th class="text-warning">Mutation in the conserved regions</th><th class="text-success">Valid sequences</th><th>Amber codon in the sequence (<u>inside</u> the area of interest)</th></tr>')
|
|
|
98 w.write(
|
|
|
99 '<tr><td class="text-error">%d sequence(s) (%.2f%%)</td><td class="text-error">%d sequence(s) (%.2f%%)</td><td class="text-error">%d sequence(s) (%.2f%%)</td><td class="text-warning">%d sequence(s) (%.2f%%)</td><td class="text-success">%d sequence(s) (%.2f%%)</td><td>%d sequence(s)</td></tr>' % (
|
|
|
100 len(no_restric), float(len(no_restric)) / float(len(all_seq)) * 100, len(no_multiple), float(len(no_multiple)) / float(len(all_seq)) * 100, len(stop),
|
|
|
101 float(len(stop)) / float(len(all_seq)) * 100, len(mut), float(len(mut)) / float(len(all_seq)) * 100, len(good_ids), float(len(good_ids)) / float(len(all_seq)) * 100,
|
|
|
102 len(amber)))
|
|
|
103 w.write(
|
|
|
104 '<tr><td class="text-error">%s</td><td class="text-error">%s</td><td class="text-error">%s</td><td class="text-warning">%s</td><td class="text-success">%s</td><td>%s</td></tr></table>' % (
|
|
|
105 '<br/>'.join(no_restric), '<br/>'.join(no_multiple), '<br/>'.join(stop), '<br/>'.join(mut), '<br/>'.join(good_ids), '<br/>'.join(amber)))
|
|
|
106 # Variable regions analysis
|
|
|
107 w.write(
|
|
|
108 '<div class="page-header"><a id="variable"></a><h2>Variable regions analysis</h2></div><p>The following group of sequences are identical clones on the variable regions:</p>')
|
|
|
109 identical_clones_seq = identical_clones.keys()
|
|
|
110 if identical_clones_seq:
|
|
|
111 for seq in identical_clones_seq:
|
|
|
112 ids = list(set(identical_clones[seq])) # return only one occurrence of each item in the list
|
|
|
113 w.write('<div class="row-fluid"><div class="span5"><pre>%d sequences (%.2f%% of valid sequences)<br/>%s</pre></div>' % (
|
|
|
114 len(ids), float(len(ids)) / float(len(good_ids)) * 100, '<br/>'.join(ids)))
|
|
|
115 w.write('<div class="span3"><table class="table table-striped table-bordered"><thead><tr><th>Variable region</th><th>Repeated sequence</th></tr></thead><tbody>')
|
|
|
116 for z in range(len(good_seq[ids[0]]['var'])):
|
|
|
117 w.write('<td>%d</td><td>%s</td></tr>' % (z + 1, good_seq[ids[0]]['var'][z]))
|
|
|
118 w.write('</tbody></table></div></div>')
|
|
|
119 else:
|
|
|
120 w.write('<p>No clone was found.</p>')
|
|
|
121
|
|
|
122 first = True
|
|
|
123 for i in range(nb_var_part):
|
|
|
124 keys = []
|
|
|
125 for k in (var_seq_common[str(i + 1)].keys()):
|
|
|
126 nb = var_seq_common[str(i + 1)][k]
|
|
|
127 if nb > 1:
|
|
|
128 if first:
|
|
|
129 w.write(
|
|
|
130 '<p>Here\'s the distribution of the repeated sequences in variable regions:</p><table class="table table-striped table-bordered"><thead><tr><th>Variable region</th><th>Repeated sequence</th><th>Number of occurrences (percentage of valid sequences)</th></tr></thead><tbody>')
|
|
|
131 first = False
|
|
|
132 keys.append(k)
|
|
|
133 else:
|
|
|
134 keys.append(k)
|
|
|
135 nb = len(keys)
|
|
|
136 if nb != 0:
|
|
|
137 w.write('<tr>')
|
|
|
138 for z in range(nb):
|
|
|
139 if z == 0:
|
|
|
140 w.write('<td rowspan="%d">%d</td>' % (nb, i + 1))
|
|
|
141 w.write('<td>%s</td><td>%d (%.2f%%)</td></tr>' % (
|
|
|
142 keys[z], var_seq_common[str(i + 1)][keys[z]], float(var_seq_common[str(i + 1)][keys[z]]) / float(len(good_ids)) * 100))
|
|
|
143 w.write('</tbody></table>')
|
|
|
144 # Clustering
|
|
|
145 w.write('<div class="page-header"><a id="cluster"></a><h2>Clustering</h2></div><p>The following clusters were generated by MCL:</p>')
|
|
|
146 for line in open(mcl_output, 'r'):
|
|
|
147 w.write('<div class="row-fluid"><div class="span6"><pre>%d sequences (%.2f%% of valid sequences)<br/>%s</pre></div></div>' % (
|
|
|
148 len(line.split("\t")), float(len(line.split("\t"))) / float(len(good_ids)) * 100, '<br/>'.join(line.split("\t"))))
|
|
|
149 # Statistics
|
|
|
150 w.write('<div class="page-header"><a id="stat"></a><h2>Statistics</h2></div>')
|
|
|
151 w.write('<p>Here\'s some statistics about the valid sequences:</p><p>Mean for the pairwise alignement scores: %.2f<br/>Standard deviation: %.2f</p>' % (
|
|
|
152 numpy.mean(align_scores), numpy.std(align_scores)))
|
|
|
153 w.write('<div class="row-fluid"><div class="span6"><img src="%s" alt="Distribution of the pairwise alignment score"></div>' % os.path.basename(graph_pic))
|
|
|
154 w.write('<div class="span6"><table class="table table-striped table-bordered"><thead><tr><th>Pairwise Alignment Score</th><th>Number of occurrences</th></tr></thead><tbody>')
|
|
|
155 uniq_scores = sorted(list(set(align_scores)))
|
|
|
156 scores_dic = {}
|
|
|
157 for score in uniq_scores:
|
|
|
158 scores_dic[score] = align_scores.count(score)
|
|
|
159
|
|
|
160 scores_dic = dict(sorted(scores_dic.items()))
|
|
|
161 scores = scores_dic.items()
|
|
|
162 # scores.sort()
|
|
|
163
|
|
|
164 for el in scores:
|
|
|
165 w.write('<tr><td>%.2f</td><td>%d</td></tr>' % (el[0], el[1]))
|
|
|
166 w.write('</tbody></table></div></div>')
|
|
|
167 # Annex
|
|
|
168 w.write('<div class="page-header"><a id="annex"></a><h2>Annex</h2></div>')
|
|
|
169 w.write('<p><strong>Valid protein sequences</strong> in FASTA format:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">')
|
|
|
170 for _id in good_ids:
|
|
|
171 w.write('>%s\n%s\n' % (_id, re.sub("(.{80})", "\\1\n", good_seq[_id]['prot'], re.DOTALL)))
|
|
|
172 w.write('</textarea>')
|
|
|
173 aln_out = generate_aln(good_seq, good_ids)
|
|
|
174 print(str(aln_out))
|
|
|
175 w.write(
|
|
|
176 '<p>Multiple sequence alignment of the <strong>valid sequences</strong> generated by Clustal Omega:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">%s</textarea>' % str(
|
|
|
177 aln_out))
|
|
|
178
|
|
|
179 if no_multiple:
|
|
|
180 w.write(
|
|
|
181 '<p><strong>Protein sequences with an incorrect number of nucleotides between the restriction sites</strong> in FASTA format:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">')
|
|
|
182 for _id in no_multiple:
|
|
|
183 w.write('>%s\n%s\n' % (_id, re.sub("(.{80})", "\\1\n", all_seq_fasta[_id]['prot'], re.DOTALL)))
|
|
|
184 w.write('</textarea>')
|
|
|
185
|
|
|
186 if mut:
|
|
|
187 w.write('<p><strong>Mutated protein sequences</strong> in FASTA format:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">')
|
|
|
188 for _id in mut:
|
|
|
189 w.write('>%s\n%s\n' % (_id, re.sub("(.{80})", "\\1\n", all_seq_fasta[_id]['prot'], re.DOTALL)))
|
|
|
190 w.write('</textarea>')
|
|
|
191 aln_out = generate_aln(all_seq_fasta, mut)
|
|
|
192
|
|
|
193 w.write(
|
|
|
194 '<p>Multiple sequence alignment of the <strong>mutated sequences</strong> generated by Clustal Omega:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">%s</textarea>' % str(
|
|
|
195 aln_out))
|
|
|
196
|
|
|
197 if stop:
|
|
|
198 w.write('<p><strong>Protein sequences with a stop codon</strong> in FASTA format:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">')
|
|
|
199 for _id in stop:
|
|
|
200 w.write('>%s\n%s\n' % (_id, re.sub("(.{80})", "\\1\n", all_seq_fasta[_id]['prot'], re.DOTALL)))
|
|
|
201 w.write('</textarea>')
|
|
|
202
|
|
|
203 if amber:
|
|
|
204 w.write('<p><strong>Protein sequences with an amber codon</strong> in FASTA format:</p><textarea class="span8 fasta" type="text" rows="20" readonly="readonly">')
|
|
|
205 for _id in amber:
|
|
|
206 w.write('>%s\n%s\n' % (_id, re.sub("(.{80})", "\\1\n", all_seq_fasta[_id]['prot'], re.DOTALL)))
|
|
|
207 w.write('</textarea>')
|
|
|
208
|
|
|
209 w.write('</div></body></html>')
|
|
|
210 w.close()
|
|
|
211
|
|
|
212
|
|
|
213 nb_seq = len(list(SeqIO.parse(args.input, "fasta")))
|
|
|
214
|
|
|
215 for seq_record in SeqIO.parse(args.input, "fasta"):
|
|
|
216 seq_id = seq_record.id
|
|
|
217 seq = str(seq_record.seq)
|
|
|
218 seq = seq.upper()
|
|
|
219 all_seq.append(seq_id)
|
|
|
220 # Checking if both restriction sites are present in the sequence
|
|
|
221 if site_res_5 in seq and site_res_3 in seq:
|
|
|
222 valid = True
|
|
|
223 else:
|
|
|
224 valid = False
|
|
|
225 tag['no_restric'].append(seq_id)
|
|
|
226 # If sequence has both restriction sites, checking if it is necessary to take the reverse complement strand
|
|
|
227 if valid:
|
|
|
228 site_res_5_pos = seq.index(site_res_5)
|
|
|
229 site_res_3_pos = seq.index(site_res_3)
|
|
|
230 # If site_res_5_pos > site_res_3_pos, reverse complement strand has to be calculated
|
|
|
231 if site_res_5_pos > site_res_3_pos:
|
|
|
232 # Checking if the number of nucleic acids between the restriction sites is a multiple of 3
|
|
|
233 length = math.fabs((site_res_5_pos + len(site_res_5)) - site_res_3_pos)
|
|
|
234 valid = length % 3 == 0
|
|
|
235 cut_seq = seq[:site_res_5_pos + len(site_res_5)]
|
|
|
236 cut_seq = reverse_complement(cut_seq)
|
|
|
237 # Else if site_res_5_pos < site_res_3_pos, use the sequence as it is
|
|
|
238 else:
|
|
|
239 # Checking if the number of nucleic acids between the restriction sites is a multiple of 3
|
|
|
240 length = math.fabs((site_res_3_pos + len(site_res_3)) - site_res_5_pos)
|
|
|
241 valid = length % 3 == 0
|
|
|
242 cut_seq = seq[site_res_5_pos:]
|
|
|
243 # If the number of nucleic acids between the restriction sites isn't a multiple of 3, put the sequence away
|
|
|
244 if not valid:
|
|
|
245 tag['no_multiple'].append(seq_id)
|
|
|
246 prot_seq = Seq.translate(cut_seq)
|
|
|
247 all_seq_fasta[seq_id] = {}
|
|
|
248 all_seq_fasta[seq_id]['prot'] = prot_seq
|
|
|
249 else:
|
|
|
250 # Translate nucleic sequence into amino acid sequence
|
|
|
251 prot_seq = Seq.translate(cut_seq)
|
|
|
252 all_seq_fasta[seq_id] = {}
|
|
|
253 all_seq_fasta[seq_id]['prot'] = prot_seq
|
|
|
254
|
|
|
255 # Looking for stop codon in the sequence and getting their position in the sequence
|
|
|
256 if '*' in prot_seq:
|
|
|
257 pos_stop = [m.start() for m in re.finditer("\*", prot_seq)]
|
|
|
258 stop = False
|
|
|
259 # Checking if stop codon is between the restriction sites, also checking if it is an amber codon. if stop codon other than amber codon -> tag stop
|
|
|
260 for i in range(len(pos_stop)):
|
|
|
261 if pos_stop[i] < length / 3:
|
|
|
262 stop_codon_nuc = cut_seq[pos_stop[i] * 3:pos_stop[i] * 3 + 3]
|
|
|
263 if stop_codon_nuc != "TAG":
|
|
|
264 tag['stop'].append(seq_id)
|
|
|
265 stop = True
|
|
|
266 break
|
|
|
267 else:
|
|
|
268 if seq_id not in tag['amber']:
|
|
|
269 tag['amber'].append(seq_id)
|
|
|
270 # If stop codon wasn't found between the restriction sites
|
|
|
271 if not stop:
|
|
|
272 """
|
|
|
273 # Checking if there is a stop codon outside the restriction sites. If yes -> tag ok_stop_ext
|
|
|
274 for i in range(len(pos_stop)):
|
|
|
275 if (pos_stop[i] > length/3):
|
|
|
276 stop_codon_nuc = cut_seq[pos_stop[i]*3:pos_stop[i]*3+3]
|
|
|
277 if stop_codon_nuc != "TAG":
|
|
|
278 tag['ok_stop_ext'].append(seq_id)
|
|
|
279 stop = True
|
|
|
280 break
|
|
|
281 else:
|
|
|
282 if (seq_id not in tag['amber']):
|
|
|
283 tag['amber'].append(seq_id)
|
|
|
284 """
|
|
|
285 # Checking if there was a mutation in the fix part, if yes -> tag mut else retrieve variable parts
|
|
|
286 mut = False
|
|
|
287 pattern_part = args.pattern.split(":")
|
|
|
288 tmp_prot_seq = prot_seq
|
|
|
289 var_parts = []
|
|
|
290 for i in range(len(pattern_part) - 1): # not checking the latest fix part
|
|
|
291 part = pattern_part[i]
|
|
|
292 # If part is fix
|
|
|
293 if not part[0].isdigit():
|
|
|
294 # If part not in prot_seq -> mutation, flag then break
|
|
|
295 if part not in tmp_prot_seq:
|
|
|
296 mut = True
|
|
|
297 tag['mut'].append(seq_id)
|
|
|
298 break
|
|
|
299 # Else, store the variable part if exist then remove the fix part + variable part (tmp_prot_seq starts at the end of part)
|
|
|
300 else:
|
|
|
301 pos_fix = tmp_prot_seq.index(part)
|
|
|
302 if pos_fix != 0:
|
|
|
303 var_parts.append(tmp_prot_seq[0:pos_fix])
|
|
|
304 tmp_prot_seq = tmp_prot_seq[pos_fix + len(part):]
|
|
|
305 # Else part is variable
|
|
|
306 else:
|
|
|
307 nb_var_part += 1
|
|
|
308 # Treating latest fix part if no mutation before
|
|
|
309 if not mut:
|
|
|
310 last_part = pattern_part[-1]
|
|
|
311 last_var = pattern_part[-2]
|
|
|
312 if '-' in last_var:
|
|
|
313 var_max = int(last_var.split('-')[1])
|
|
|
314 else:
|
|
|
315 var_max = int(last_var)
|
|
|
316 last_part = last_part[0:var_max + 1]
|
|
|
317 if last_part not in tmp_prot_seq:
|
|
|
318 mut = True
|
|
|
319 tag['mut'].append(seq_id)
|
|
|
320 else:
|
|
|
321 pos_fix = tmp_prot_seq.index(last_part)
|
|
|
322 if pos_fix != 0:
|
|
|
323 var_parts.append(tmp_prot_seq[0:pos_fix])
|
|
|
324 # If no mutation the sequence is validated and all the info are stored
|
|
|
325 if not mut:
|
|
|
326 good_seq[seq_id] = {}
|
|
|
327 good_seq[seq_id]['dna'] = cut_seq
|
|
|
328 good_seq[seq_id]['prot'] = prot_seq
|
|
|
329 good_seq[seq_id]['var'] = var_parts
|
|
|
330
|
|
|
331 # If all sequences are invalid, the program will exit as there is no data to continue
|
|
|
332 if not good_seq:
|
|
|
333 print("All sequences are invalid. At least 2 valid sequences are necessary to proceed to the next step. The program will now exit.")
|
|
|
334 sys.exit()
|
|
|
335 elif len(good_seq.keys()) == 1:
|
|
|
336 print("There is only one valid sequence among the input data. At least 2 valid sequences are necessary to proceed to the next step. The program will now exit")
|
|
|
337 sys.exit()
|
|
|
338
|
|
|
339 # Initialization of dict var_seq_common
|
|
|
340 for n in range(nb_var_part):
|
|
|
341 var_seq_common[str(n + 1)] = {}
|
|
|
342
|
|
|
343 # Opening the file where the mcl input will be written
|
|
|
344 mcl = open(mcl_file, 'w')
|
|
|
345
|
|
|
346 id = good_seq.keys()
|
|
|
347 for i in range(len(id)):
|
|
|
348 var_1 = good_seq[list(id)[i]]['var']
|
|
|
349
|
|
|
350 # Classifying variable sequences
|
|
|
351 for k in range(len(var_1)):
|
|
|
352 try:
|
|
|
353 var_seq_common[str(k + 1)][var_1[k]] += 1
|
|
|
354 except KeyError:
|
|
|
355 var_seq_common[str(k + 1)][var_1[k]] = 1
|
|
|
356
|
|
|
357 for j in range(i + 1, len(id)):
|
|
|
358 var_2 = good_seq[list(id)[j]]['var']
|
|
|
359 # Comparing the sequences' variable parts to find identical clones
|
|
|
360 if var_1 == var_2:
|
|
|
361 try:
|
|
|
362 s = "".join(var_1)
|
|
|
363 identical_clones[s].extend([id[i], id[j]])
|
|
|
364 except KeyError:
|
|
|
365 identical_clones[s] = [id[i], id[j]]
|
|
|
366
|
|
|
367 # Align the 2 sequences using NWalign_PAM30
|
|
|
368 seq_1 = ''.join(var_1)
|
|
|
369 seq_2 = ''.join(var_2)
|
|
|
370 print(seq_1)
|
|
|
371 print(seq_2)
|
|
|
372 matrix = matlist.pam30
|
|
|
373 cpt = 0
|
|
|
374 if len(seq_2) > len(seq_1):
|
|
|
375 print(pairwise2.align.globalds(seq_1, seq_2, matrix, -11, -1))
|
|
|
376 for a in pairwise2.align.globalds(seq_1, seq_2, matrix, -11, -1):
|
|
|
377 for k in range(a[4]):
|
|
|
378 if a[0][k] == a[1][k]:
|
|
|
379 cpt += 1
|
|
|
380 print(format_alignment(*a, full_sequences=True))
|
|
|
381 else:
|
|
|
382 print(pairwise2.align.globalds(seq_2, seq_1, matrix, -11, -1))
|
|
|
383 for a in pairwise2.align.globalds(seq_2, seq_1, matrix, -11, -1):
|
|
|
384 for k in range(a[4]):
|
|
|
385 if a[0][k] == a[1][k]:
|
|
|
386 cpt += 1
|
|
|
387 print(format_alignment(*a, full_sequences=True))
|
|
|
388 print("######################################@")
|
|
|
389 print(cpt)
|
|
|
390
|
|
|
391 if len(seq_2) > len(seq_1):
|
|
|
392 p = subprocess.Popen(fibronectin_script_dir + "/NWalign_PAM30 %s %s 3" % (seq_1, seq_2), shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
|
|
|
393 else:
|
|
|
394 p = subprocess.Popen(fibronectin_script_dir + "/NWalign_PAM30 %s %s 3" % (seq_2, seq_1), shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
|
|
|
395
|
|
|
396 out, err = p.communicate()
|
|
|
397
|
|
|
398 print(out)
|
|
|
399 print("######################################@")
|
|
|
400 lines = out.split(bytes("\n", encoding='utf8'))
|
|
|
401 print(lines[5].split(bytes(' ', encoding='utf8'))[5])
|
|
|
402 score = float(lines[5].split(bytes(' ', encoding='utf8'))[5]) * 100
|
|
|
403 align_scores.append(score)
|
|
|
404 mcl.write('%s\t%s\t%0.2f\n' % (list(id)[i], list(id)[j], score))
|
|
|
405 mcl.close()
|
|
|
406
|
|
|
407 # Clusters formation
|
|
|
408 subprocess.call("mcl %s --abc -I 6.0 -o %s" % (mcl_file, mcl_output), shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
|
|
|
409
|
|
|
410 # Producing distribution graph
|
|
|
411 plot.hist(align_scores, bins=numpy.arange(0, 101, 2))
|
|
|
412 plot.xlabel('Pairwise Alignment Score')
|
|
|
413 plot.ylabel('Number of occurrences')
|
|
|
414 plot.title('Distribution of the pairwise alignment score')
|
|
|
415 plot.grid(True)
|
|
|
416 plot.savefig(graph_pic)
|
|
|
417
|
|
|
418 # Generating html report
|
|
|
419 report_html(html_file, tag, all_seq, good_seq, all_seq_fasta, identical_clones, nb_var_part, var_seq_common, align_scores, args)
|
|
|
420
|
|
|
421 # Removing intermediate files
|
|
|
422 subprocess.call("rm %s %s " % (mcl_file, mcl_output), shell=True)
|
|
|
423
|
|
|
424 print("HTML report has been generated in the output directory. The program will now exit.")
|
