· 8 years ago · Jun 03, 2018, 04:54 PM
1package Greenphyl::Load::Family;
2
3use strict;
4use warnings;
5no warnings 'once';
6use Carp;
7use Bio::Seq;
8use Bio::SearchIO;
9use Bio::SeqIO;
10use lib "../../lib";
11use Greenphyl::Path;
12use Greenphyl::SequenceFilter 'filter_text_id';
13
14umask 0000; # set the appropriate permissions for created files/folder to anydody as deamon gives only read permission to others
15
16my ( $mysql, $dbh );
17
18sub new {
19 my ( $class, %p ) = @_;
20 my $self = {};
21
22 $mysql = $p{mysql};
23 $dbh = $p{dbh};
24
25 bless( $self, $class );
26 return $self;
27}
28
29#calculate the number of sequence in the fasta input file
30## inefficient. could be done with a grep in 1 line.
31## qx/ grep '>' $file | wc -l /;
32
33sub nbSequence {
34 my ( $file ) = @_;
35 my $total_size = 0;
36 open( IN, $file );
37
38 while ( <IN> ) {
39 $total_size++ if /^>/;
40 }
41
42 close( IN );
43 return $total_size;
44}
45
46# $liste_family is a text file containing family name that were modified and created
47# this list will be useful for the maintenance when running MEME and the phylo pipeline
48
49# premier hit ignore car sequence contre elle-meme = meilleur hit
50# second hit : -> soit match contre une famille -> si longueur bonne integration dans la famille
51# -> sinon ORPHAN
52# -> match avec ORPHAN -> longueur bonne = nouvelle famille et on regarde les hit suivant -> si longueur bonne integration dans la nouvelle famille
53# -> autre arret de l'algo (longueur mauvaise ou plus de hit)
54# -> longueur mauvaise = ORPHAN
55# -> math avec rien -> ORPHAN
56
57sub loadFamily {
58 my ( $self, $file, $count, $evalue, $database, $out_number, $file_out, $file_temp, $species_name, $restart, %family_rem ) = @_;
59
60 my $pass = new Greenphyl::Path( $species_name );
61 my $cpt = 0;
62 my $cpt_new_family = 0;
63 my $cpt_modif_family = 0;
64 my $cpt_progress = 0;
65 my $nb_sequence = nbSequence( $file );
66
67 my %liste_pourcentage;
68
69 open( LOG, ">> $Greenphyl::Path::data" )
70 or die "cannot create the file $Greenphyl::Path::data $!";
71 open( FAMILY, ">> $pass->{file_result}$Greenphyl::Path::liste_family " )
72 || die( "error not create the file $pass->{file_result}$Greenphyl::Path::liste_family\n" );
73 open( NFAM, ">> $pass->{file_result}$Greenphyl::Path::new_family" )
74 || die( "error not create the file $pass->{file_result}$Greenphyl::Path::new_family\n" );
75 open( AFF, ">> $Greenphyl::Path::data_path/$Greenphyl::Path::log" )
76 || die( "error not create the file $Greenphyl::Path::data_path/$Greenphyl::Path::log\n" );
77 my $in = Bio::SeqIO->newFh( -file => $file, -format => 'fasta' );
78
79 while ( <$in> ) {
80
81 # extract a single query sequence from the fasta file with the newly added sequences
82
83 my $id = filter_text_id( $_->display_id ); #query
84 my $seq = $_->seq;
85 my $fasta = $_;
86 my $cpt_hit = 0; #counter for the number of hit with blast
87 my $hit_while = 0; #to last the loop
88
89 $cpt_progress++;
90
91 if ( $cpt_progress > $restart ) {
92 print "\n----------------------------------------------\n";
93 print AFF "Gene family modified: " . $cpt_progress . " / " . $nb_sequence . "\n";
94 print $cpt_progress. " / " . $nb_sequence . " => $id \n";
95
96 $cpt++;
97 ##check if query is in the database
98 my $sql = "SELECT seq_id FROM sequences WHERE seq_textid = '$id'";
99 my @rep = $mysql->select( $sql );
100
101 unless ( $rep[0] ) {
102 print "ERROR : sequence is not in the db!<br />";
103 last;
104 }
105
106 my $seq = $_->seq;
107
108 #++++ should be in run::blast
109 ##creation input fasta file as input for blast software
110 $file_temp = clipboard( $seq, $file_temp );
111 system( "chmod 777 $file_temp" );
112 $count++;
113
114 # run blast to compare it to blast database containing other sequences across all species
115 my $cmd = "$Greenphyl::Config::BLASTALL -i $file_temp -e $evalue -p blastp -F none -d $database -K $out_number -v $out_number -b $out_number -o $file_out";
116 system( $cmd );
117
118 my $found = 0;
119 if ( -e $file_out ) # if blast output exist
120 {
121
122 # load blast results
123 my $in = new Bio::SearchIO(
124 -format => 'blast',
125 -file => "$file_out"
126 );
127 system( "chmod 666 $file_out" );
128
129 #++++
130
131 # collect all sequences that are similar to this one from the blast result
132 # and see if the query sequence can be added to a family or combined with another sequence to form a new family
133 # making sure that it is not accidentally combined with itself
134 while ( my $result = $in->next_result and $hit_while == 0 ) {
135 my $level1 = 0;
136 my $cpt = 0;
137
138 while ( my $hit = $result->next_hit
139 and $hit_while == 0 ) #sortir ici
140 {
141 my $seq_textid = filter_text_id( $hit->name );
142 $cpt_hit++;
143
144 if ( $seq_textid ne $id ) # ignore first hit
145 {
146 $found++;
147
148 $seq_textid = ucfirst( lc( $seq_textid ) );
149 my $sql = "SELECT S.seq_textid, F.family_name, F.family_id, FO.class_id
150 FROM sequences S, seq_is_in SI, family F, found_in FO
151 WHERE S.seq_textid =\"$seq_textid\" and S.seq_id = SI.seq_id and SI.family_id = F.family_id and F.family_id= FO.family_id";
152 my @rep7 = $mysql->select( $sql );
153 my $count_verif; #pour verifier la taille de la sequence
154
155 # hit found in an existing family
156 if ( @rep7 ) {
157 print "Find family for $seq_textid \t";
158 my $family_id = $rep7[0][2]; #retrieve family id at level 1
159 $count_verif = &count( $family_id, $id ); #verifier si ce sont les bons parametres
160
161 if ( $count_verif == 1 ) # length ok
162 {
163 ## retrieve family_id at level 1 down to 4 (if exists)
164 for my $i ( 0 .. $#rep7 ) {
165 my $family_id = $rep7[$i][2];
166
167 if ( $found == 1 ) {
168
169 #retrieve seq_id for the sequence to be inserted in seq_is_in
170 my $sql = "SELECT seq_id FROM sequences WHERE seq_textid = '" . $id . "'";
171 my @rep = $mysql->select( $sql );
172 my $seq_id = $rep[0];
173
174 ## check if the sequence is not already classified
175 my $sql1 = "SELECT SI.family_id
176 FROM seq_is_in SI, found_in F
177 WHERE SI.seq_id = '"
178 . $seq_id . "' and F.family_id = SI.family_id and F.class_id = '" . $rep7[$i][3] . "'";
179 my @rep2 = $mysql->select( $sql1 );
180
181 ##if classified
182 if ( $rep2[0] ) {
183
184 #Deja classifiee au niveau $rep7[$i][3] !!<br />\n";
185 last;
186 }
187
188 #classify second hit
189 else {
190 print FAMILY "$family_id\n";
191 print LOG "Insertion of the sequence $id in the $family_id family\n";
192 %liste_pourcentage = &dataFamily( $family_id, %liste_pourcentage );
193 my $sql2 = "INSERT INTO seq_is_in (seq_id, family_id) VALUES(\"$seq_id\",\"$family_id\")";
194 $mysql->insert( $sql2 );
195 $cpt_modif_family++;
196 }
197 }
198 }
199 }
200 else #length not ok : query = orphan
201 {
202 $hit_while = 1; #last the loop
203 next;
204 }
205 }
206
207 # hit found est un orphan
208 else {
209 ##recup seq_id de query
210 my $sql = "SELECT seq_id FROM sequences WHERE seq_textid = '$id'";
211 my @rep = $mysql->select( $sql );
212 my $seq_id = $rep[0];
213
214 ## verifie que query deja pas classifié dans la base avec 1er orphan de la sortie du blast
215 my $sql2 = "SELECT family_id FROM seq_is_in WHERE seq_id = '$seq_id'";
216 my @rep2 = $mysql->select( $sql2 );
217 my $family_query = $rep2[0];
218
219 #pour premier orphan trouve cree une famille, voir else
220
221 ##deja classe:va voir si autre sorties de blast sont des orphans et si oui les classe dans meme famille
222 if ( $family_query ) {
223 ##recup id des autres sortie du blast
224 my $sql = "SELECT seq_id FROM sequences WHERE seq_textid = '$seq_textid'";
225 my @rep = $mysql->select( $sql );
226 my $seq_id = $rep[0];
227
228 if ( $seq_id ) {
229 ## verifie les autres sortie blast
230 my $sql = "SELECT family_id FROM seq_is_in WHERE seq_id = '$seq_id'";
231 my @rep = $mysql->select( $sql );
232 my $family_id2 = $rep[0];
233
234 if ( $family_id2 ) {
235 $hit_while = 1; #on va pas regarder les autres hits, on stop la boucle ici.
236 }
237 else #sortie est un orphan
238 {
239 $count_verif = &count( $family_query, $id ); #test la taille de la sequence
240
241 if ( $count_verif == 1 ) #le test de taille est concluant
242 {
243 %liste_pourcentage = &dataFamily( $family_query, %liste_pourcentage );
244 my $sql2 = "INSERT INTO seq_is_in (seq_id, family_id) VALUES(\"$seq_id\",\"$family_query\")";
245 $mysql->insert( $sql2 );
246 $dbh->do( "TRUNCATE go_family_cache;" );
247 print FAMILY "$family_query\n";
248 $cpt_new_family++;
249 }
250 }
251 }
252 }
253
254 #pas de famille cree : premier passage dans la boucle : creation d'une nouvelle famille avec l'orphan
255 #match avec un orphan
256 else {
257 print "new_family with $seq_textid\t";
258
259 ###compare sequence length between the query and the hit
260 my $sql00 = "SELECT seq_length FROM sequences WHERE seq_textid = \"$id\""; #query
261 my $sql01 = "SELECT seq_length FROM sequences WHERE seq_textid = \"$seq_textid\""; #hit
262 my @rep00 = $mysql->select( $sql00 );
263 my @rep01 = $mysql->select( $sql01 );
264
265 # should be set at the begining of the script
266 my $soixante = ( $rep01[0] * 60 ) / 100;
267 my $centquarante = ( $rep01[0] * 140 ) / 100;
268
269 # TODO: this size comparison is pretty broken and will need to be replaced once cleaning is done ( minimum to value is wider range than max to value )
270 if ( ( $rep00[0] < $centquarante )
271 and ( $rep00[0] > $soixante ) ) #if length control is ok
272 {
273
274 # va chercher max family_number pour creer nouvelle famille
275 my $sql1 = 'SELECT MAX(family_number) FROM family';
276 my @rep = $mysql->select( $sql1 );
277 my $old_number = $rep[0] || 0;
278 my $family_number = $old_number + 1;
279
280 ## create a new family in the db table family
281 my $new = "New created family by $species_name";
282
283 my $sql = "INSERT INTO family (family_name, family_number, inference) VALUES(\"$new\", \"$family_number\", '')";
284
285 $mysql->insert( $sql );
286 $dbh->do( "TRUNCATE go_family_cache;" );
287
288 ## recup le family_id creer
289
290 my $sql3 = "SELECT family_id FROM family WHERE family_number LIKE \"$family_number\"";
291
292 my @rep1 = $mysql->select( $sql3 );
293 my $family_id = $rep1[0];
294
295 ## new family is classify at level 1
296 my $nb = 1;
297 my $sql4 = "INSERT INTO found_in (family_id, class_id) values(\"$family_id\", \"$nb\")";
298
299 $mysql->insert( $sql4 );
300 $dbh->do( "TRUNCATE go_family_cache;" );
301
302 print "New family level 1: Family_id: $family_id\n";
303
304 ## 2 sequences in the new family
305
306 ##Query
307 my $sql7 = "SELECT seq_id FROM sequences WHERE seq_textid = '" . $seq_textid . "'";
308 my @rep3 = $mysql->select( $sql7 );
309 my $seq_id = $rep3[0];
310 my $sql2 = "INSERT INTO seq_is_in (seq_id, family_id) values(\"$seq_id\",\"$family_id\")";
311 $mysql->insert( $sql2 );
312 $dbh->do( "TRUNCATE go_family_cache;" );
313
314 ##Ancien orphan
315 my $sql6 = "SELECT seq_id FROM sequences WHERE seq_textid = '" . $id . "'";
316 my @rep2 = $mysql->select( $sql6 );
317 my $seq_id2 = $rep2[0];
318 my $sql5 = "INSERT INTO seq_is_in (seq_id, family_id) values(\"$seq_id2\",\"$family_id\")";
319 $mysql->insert( $sql5 );
320 $dbh->do( "TRUNCATE go_family_cache;" );
321
322 %liste_pourcentage = &dataFamily( $family_id, %liste_pourcentage );
323 print FAMILY "$family_id\n";
324 print NFAM "$family_id\n";
325
326 }
327 else { #taille est incorrecte = orphan
328 print " too small ORPHAN \n";
329 $hit_while = 1;
330 }
331 }
332 }
333 } #end control if hit diff query
334 } # end boucle sur list HIT
335 }
336 }
337 }
338 }
339
340 $dbh->do( "call countFamilySeqBySpecies()" ) if $ENV{TESTING}; # TODO: this needs to be removed and count_family_seq updated properly!
341
342 print LOG "number of families amended : $cpt_modif_family \n";
343 print LOG "number of families created : $cpt_new_family \n";
344
345 close FAMILY;
346 close AFF;
347
348 &pourcentage( $species_name, \%liste_pourcentage, \%family_rem );
349 return $cpt_new_family, $cpt_modif_family;
350}
351
352#track the number of sequence added to the family
353#by tracking the number at the beginning and added, we can calculate the number of sequence removed
354#allow to calculate the percentage of modification
355sub dataFamily {
356 my ( $family_id, %liste_pourcentage ) = @_;
357
358 if ( keys( %liste_pourcentage ) ) {
359
360 if ( exists $liste_pourcentage{$family_id} ) #family was modified at least once
361 {
362 $liste_pourcentage{$family_id} = $liste_pourcentage{$family_id} + 1;
363 }
364 else {
365 my $sql = "SELECT seq_id FROM seq_is_in WHERE family_id = $family_id";
366 my @rep = $mysql->select( $sql );
367 my $old_size = scalar( @rep );
368 $liste_pourcentage{$family_id} = 1;
369 }
370 }
371
372 else {
373 my $sql = "SELECT seq_id FROM seq_is_in WHERE family_id = $family_id";
374 my @rep = $mysql->select( $sql );
375 $liste_pourcentage{$family_id} = 1;
376 }
377
378 return %liste_pourcentage;
379}
380
381#track the number of changes in a family (number of new and/or removed sequences)
382#calculate the percentage of modification for a family
383# %liste_family contains id of each famiily with new sequences members : key = family_id value = number of sequences added
384# %family contains id of each famiily where sequences were removed
385
386sub pourcentage (\%) #BAD way to specify paramter type
387{
388 my ( $species_name, $ref_liste_pourcentage, $ref_family_rem ) = @_;
389 my %liste_pourcentage = %$ref_liste_pourcentage;
390 my %family_rem = %$ref_family_rem;
391
392 my $pass = new Greenphyl::Path( $species_name );
393
394 open( POURCENTAGE, "> $pass->{file_result}$Greenphyl::Path::pourcentage " )
395 || die( "error not create the file $pass->{file_result}$Greenphyl::Path::liste_family\n" );
396 for my $key ( keys %liste_pourcentage ) #regarde toutes les familles modifiees avec des ADD
397 {
398 my $sql = "SELECT seq_id FROM seq_is_in WHERE family_id = \"$key\"";
399 my @rep = $mysql->select( $sql );
400 my $new_size = scalar( @rep ); #family size
401
402 if ( exists $family_rem{$key} ) #famille a ete modifiee en + et -, on a une valeur de cle commune a %liste_family et %family_rem
403 {
404 my $old_size = $new_size + $liste_pourcentage{$key} - $family_rem{$key}; #ancienne valeur de la taille
405 my $pourcentage = ( ( $liste_pourcentage{$key} + $family_rem{$key} ) / $old_size ) * 100;
406 print POURCENTAGE "Family_id => " . $key . "\ttaille origine =>" . $old_size . "\t taille actuelle => " . $new_size . "\t Number of add= > " . $liste_pourcentage{$key} . "\tnumber of remove => " . $family_rem{$key} . "\t pourcentage of modification => " . $pourcentage . "\n";
407 delete( $family_rem{$key} ); #efface de la table des familles - cette entree car deja traitee
408 }
409
410 else #famille a ete modifiee en + seulement
411 {
412 my $old_size = $new_size - $liste_pourcentage{$key}; #ancienne taille = taille actuelle - ce qui a ete ajoute
413 my $pourcentage = ( $liste_pourcentage{$key} / $old_size ) * 100;
414 print POURCENTAGE "Family_id => " . $key . "\ttaille origine =>" . $old_size . "\t taille actuelle = > " . $new_size . "\t Number of add => " . $liste_pourcentage{$key} . "\tnumber of remove = /\t pourcentage of modification => " . $pourcentage . "\n";
415 }
416 }
417
418 for my $key ( keys %family_rem ) # que ce qui a ete modifie en - (car on a enleve ce qui etait commun avec %liste_family
419 {
420 my $sql = "SELECT seq_id FROM seq_is_in WHERE family_id = $key";
421 my @rep = $mysql->select( $sql );
422 my $new_size = scalar( @rep ); #taille de la famille
423 my $old_size = $new_size + $family_rem{$key};
424 my $pourcentage = ( $family_rem{$key} / $old_size ) * 100;
425 print POURCENTAGE "Family_id => " . $key . "\ttaille origine => " . $old_size . "\t taille actuelle => " . $new_size . "\t Number of add => /\tnumber of remove => " . $family_rem{$key} . "\t pourcentage of modification => " . $pourcentage . "\n";
426 }
427 close POURCENTAGE;
428}
429
430# if new_size = 0 obsolete family was probably removed from the db!
431
432# function duplicated in many scripts, i don't know how to get rid of that function; TODO
433sub clipboard {
434 my ( $clipboard, $file_temp ) = @_;
435 my @donnees = ( split( /\n/, $clipboard ) );
436 my ( @sequence, $name, $seq );
437 $name ||= '';
438 open( OUT, ">$file_temp" ) or die "cannot open file : $!";
439
440 for my $i ( 0 .. $#donnees ) {
441
442 if ( substr( $donnees[$i], 0, 1 ) eq ">" ) { $name = "$donnees[$i]"; }
443 else {
444 $seq .= $donnees[$i];
445 }
446
447 print OUT $name . $seq . "\n";
448 $name = "";
449 $seq = "";
450 }
451 close OUT;
452 return $file_temp;
453
454}
455
456=pod
457
458=head2 count
459
460B<Description>: perform control length on the sequences. sequences are inserted in the classification if the pass the test.
461
462B<ArgsCount>: 2
463
464=over 4
465
466=item $seq_textid: (string)
467
468sequence name of the query
469
470=item $family_id: (int)
471
472family identifier in which the query sequence may be classified
473
474=back
475
476=cut
477
478sub count {
479 my ( $family_id, $seq_textid ) = @_;
480
481 my $cpt = 0;
482 my $somme_length = 0;
483
484 my %result;
485
486 my $sql = "SELECT seq_length FROM sequences WHERE seq_textid = \"$seq_textid\"";
487 my @rep = $mysql->select( $sql );
488
489 if ( @rep ) {
490 my $seq_length = $rep[0];
491
492 #select les sequences qui appartiennent a la meme famille que la seq que l'on veut etudier (pour le niveau 1 uniquement)
493 my $sql = "SELECT S.seq_textid, S.seq_length FROM sequences S, seq_is_in SI WHERE S.seq_id = SI.seq_id and SI.family_id = $family_id";
494 my @rep = $mysql->select( $sql );
495
496 if ( @rep ) #si trouve des sequences (si pas trouve alors que toutes les sequences de ALL sont dans la base = probleme !!!!!!)
497 {
498 for my $i ( 0 .. $#rep ) {
499 $cpt++;
500 $result{$cpt}{'length'} = $rep[$i][1]; #recupere pour chaque sequence de la famille dans une hashatable, la taille de chaque sequence
501 }
502
503 for my $test ( keys %result ) {
504 $somme_length += $result{$test}{'length'}; #fait la somme des longueurs des sequences pour la famille
505 }
506
507 my $moy_length = $somme_length / $cpt; #fait la moyenne
508
509 # should be more visible for configuration and should work with median
510 my $soixante = ( $moy_length * 60 ) / 100;
511 my $centquarante = ( $moy_length * 140 ) / 100;
512
513 if ( $seq_length < $soixante ) {
514 print "Too small sequence $seq_textid !!! $seq_length < $soixante \n";
515 return 0; # orphan, too small, sequence not integrated
516 }
517 if ( $seq_length > $centquarante ) {
518 print "Too large sequence $seq_textid !!! $seq_length > $centquarante \n";
519 return 0; #mauvais c'est un orphan, il est trop grand, on integre pas cette sequence.
520 }
521 else {
522 print "OK min $soixante max $centquarante for $seq_length\n";
523 return 1; #bon, on integre cette sequence.
524 }
525 }
526 else {
527 print "ERROR: the sequence is not in the Database ! \n";
528 exit;
529 }
530 }
531 else {
532 print "cannot find length for $seq_textid , maybe the sequence is not in the database\n";
533 }
534}
535
5361;