lib/SGN/Controller/InSilicoPcr.pm

   1
   2 package SGN::Controller::InSilicoPcr;
   3
   4 use Moose;
   5 use File::Temp qw | tempfile |;
   6 use File::Basename qw | basename |;
   7 use Storable qw | nstore retrieve |;
   8 use Bio::Graphics::Gel;
   9 use File::NFSLock qw | uncache |;
  10 use File::Copy qw | copy |;
  11 use Data::Dumper;
  12
  13 use Config::Any;
  14 use Try::Tiny;
  15 use Tie::UrlEncoder; our %urlencode;
  16 use File::Spec qw | catfile |;
  17 use File::Slurp qw | read_file write_file |;
  18 use CXGN::Tools::Run;
  19 use Bio::Seq;
  20 use CXGN::Page::FormattingHelpers qw/ html_break_string/;
  21
  22
  23 BEGIN { extends 'Catalyst::Controller'; }
  24
  25 sub AUTO {
  26     my $self = shift;
  27     my $c = shift;
  28     SGN::Schema::BlastDb->dbpath($c->config->{blast_db_path});
  29 }
  30
  31 sub index :Path('/tools/in_silico_pcr/') :Args(0) {
  32   my $self = shift;
  33   my $c = shift;
  34
  35   my $db_id = $c->req->param('db_id');
  36
  37   my $seq = $c->req->param('seq');
  38   my $schema = $c->dbic_schema("SGN::Schema");
  39
  40   my $group_rs = $schema->resultset("BlastDbGroup")->search( name => "In Silico PCR datasets");
  41   # my $group_rs = $schema->resultset("BlastDbGroup")->search( undef, { order_by=>'ordinal' });
  42
  43   my $databases = {};
  44   my $dataset_groups = [];
  45
  46   my $preselected_database = $c->config->{preselected_blastdb};
  47   my $preselected_category = '';
  48
  49   if ($db_id) {
  50     my $rs = $schema->resultset("BlastDb")->search( { blast_db_id => $db_id }, { join => 'blast_db_group' });
  51
  52     if ($rs == 0) {
  53       $c->throw( is_error => 0, message => "The blast database with id $db_id could not be found.");
  54     }
  55
  56     $preselected_database = $rs->first()->blast_db_id(); # first database of the category
  57     $preselected_category = $rs->first()->blast_db_group_id();
  58   }
  59
  60   foreach my $g ($group_rs->all()) {
  61     my @blast_dbs = $g->blast_dbs();
  62     push @$dataset_groups, [ $g->blast_db_group_id, $g->name() ];
  63
  64     my @dbs_AoA;
  65
  66     foreach my $db (@blast_dbs) {
  67       push @dbs_AoA, [ $db->blast_db_id(), $db->title(), $db->type() ];
  68     }
  69
  70     my @arr = sort {$a->[1] cmp $b->[1]} @dbs_AoA;
  71     $databases->{ $g->blast_db_group_id } = \@arr;
  72   }
  73
  74   $c->stash->{preselected_database} = $preselected_database;
  75   $c->stash->{preselected_category} = $preselected_category;
  76   $c->stash->{preload_id} = $c->req->param('preload_id');
  77   $c->stash->{preload_type} = $c->req->param('preload_type');
  78
  79   $c->stash->{databases} = $databases;
  80   $c->stash->{dataset_groups} = $dataset_groups;
  81
  82   $c->stash->{template} = '/tools/in_silico_pcr/in_silico_pcr.mas';
  83 }
  84
  85
  86 sub _reverse_complement{
  87         my $seq = shift;
  88         my $rev_seq = reverse $seq;
  89   $rev_seq =~ tr/ACGTacgt/TGCAtgca/;
  90
  91         return $rev_seq;
  92 }
  93
  94
  95 sub run_pcr_blast :Path('/tools/pcr_results') :Args(0) {
  96   my ($self, $c) = @_;
  97
  98   my @errors; #to store erros as they happen
  99
 100   my $fprimer = $c->req->param("fprimer");
 101   my $rprimer = $c->req->param("rprimer");
 102   my $productLength = $c->req->param("productLength");
 103   my $allowedMismatches = $c->req->param('allowedMismatches');
 104   my $frevcom = $c->req->param('frevcom'); #forward primer reverse complement
 105   my $rrevcom = $c->req->param('rrevcom'); #reverse primer reverse complement
 106
 107   my $params = $c->req->params();
 108
 109   #getting the length of the primers
 110   my $flength = length($fprimer);
 111   my $rlength = length($rprimer);
 112   $params->{flength} = $flength;
 113   $params->{rlength} = $rlength;
 114
 115   #reverse complement if checked
 116   if ($frevcom){
 117         $fprimer = _reverse_complement($fprimer);
 118   }
 119   if ($rrevcom){
 120         $rprimer = _reverse_complement($rprimer);
 121   }
 122
 123   # print STDERR "fprimer: $fprimer, rprimer: $rprimer\n";
 124   # print STDERR "DB id: $blast_db_id\n";
 125
 126   #validating productLength
 127   if ($productLength <= 0 or $productLength !~ /^[\d]*$/g) {
 128     push (@errors , "Max Product Length should be a positive digit\n");
 129   }
 130
 131   #validating AllowedMismatches
 132   if ($allowedMismatches < 0 or $allowedMismatches !~ /^[\d]*$/g) {
 133     push (@errors , "Allowed mismatches should be a positive digit\n");
 134   }
 135
 136   # return errors
 137   if (scalar (@errors) > 0){
 138     $c->stash->{errors} = join("<BR>" , @errors);
 139     $c->stash->{template} = '/tools/in_silico_pcr/insilicopcr_output.mas';
 140     return;
 141   }
 142
 143   # giving the primers a fasta format
 144   $fprimer = ">FORWARD-PRIMER\n$fprimer";
 145   $rprimer = ">REVERSE-PRIMER\n$rprimer";
 146
 147   my $sequence = "$fprimer\n$rprimer\n";
 148   $sequence =~ s/^\s+|\s+$|\n\s*\n//g; #< trim out leading and trailing whitespace and blank lines
 149
 150   # create blast input and output file
 151   my $blast_tmp_output = $c->config->{cluster_shared_tempdir}."/blast";
 152   mkdir $blast_tmp_output if ! -d $blast_tmp_output;
 153
 154   my ($seq_fh, $seqfile) = tempfile(
 155         "blast_XXXXXX",
 156         DIR=> $blast_tmp_output,
 157         );
 158
 159   print STDERR "Opening file for sequence ($seqfile)... ";
 160   open(my $FH, ">", $seqfile) || die "Can't open file for query ($seqfile)\n";
 161   print $FH $sequence if $sequence;
 162   close($FH);
 163
 164   my $jobid = basename($seqfile);
 165   print STDERR "JOB ID CREATED: $jobid\n";
 166
 167   print STDERR "Done.\n";
 168
 169
 170
 171   # get matrix
 172   my $m = $params->{matrix};
 173   $m =~ /^(BLOSUM|PAM)\d+$/ or $c->throw( is_error => 0, message => "invalid matrix '$m'" );
 174
 175   # get evalue
 176   my $evalue = $params->{evalue} ? $params->{evalue} : 1;
 177
 178   # get filter
 179   my $filter = $params->{filterq} ? 'T' : 'F';
 180
 181   # get blast database
 182   my $schema = $c->dbic_schema("SGN::Schema");
 183   my $bdb = $schema->resultset("BlastDb")->find($params->{database} ) or die "could not find bdb with file_base '$params->{database}'";
 184   my $basename = File::Spec->catfile($c->config->{blast_db_path},$bdb->file_base());
 185
 186
 187   my $result_file =  $seqfile.".out";
 188   # my $result_file = File::Spec->catfile($c->config->{basepath}, $c->tempfiles_subdir('blast'), $jobid.".out");
 189
 190   # create blast command
 191   my $blast_cmd = "blastall -p blastn -i $seqfile -d $basename -m 8 -M $m -F $filter -e $evalue -o $result_file";
 192
 193   print STDERR "COMMAND: $blast_cmd\n";
 194   my $blast_error = system($blast_cmd);
 195
 196   if ($blast_error) {
 197     print STDERR "An error occurred! $blast_error\n";
 198     $c->stash->{rest} = { error => $blast_error };
 199   }
 200
 201   my $pcr_img_path = $c->tempfiles_subdir('temp_images');
 202   my $gel_img_tempdir = $c->path_to( $c->tempfiles_subdir('temp_images') );
 203
 204   my ($pcr_seq,$gel_url) = _blast_to_pcr($result_file,$params,$gel_img_tempdir,$basename);
 205
 206
 207   # print STDERR "result_file: ".read_file($result_file)."\n";
 208   print STDERR "gel_img_tempdir: $gel_img_tempdir\n";
 209   print STDERR "gel_url: ".$pcr_img_path."/".$gel_url."\n";
 210
 211   $c->stash->{blast_res} = read_file($result_file);
 212   $c->stash->{pcr_seq} = $pcr_seq;
 213   $c->stash->{gel_url} = $pcr_img_path."/".$gel_url;
 214   $c->stash->{template} = '/tools/in_silico_pcr/insilicopcr_output.mas';
 215 }
 216
 217 sub _blast_to_pcr {
 218   my $blast_file = shift;
 219   my $params = shift;
 220   my $gel_img_tempdir = shift;
 221   my $basename = shift;
 222
 223   #Parsing the blast m8 result file
 224
 225   #Report parsing method was taken from /util/blast_result_handle.pl by Aure
 226
 227   my (%query_id, %subject_id, %identity, %align_length, %mismatches, %gap_openings, %q_start, %q_end, %s_start, %s_end,
 228         %e_value, %hit_score, %orientation);
 229
 230   my (@fprimer_ids ,@rprimer_ids);
 231
 232   open(my $res_fh, "<", $blast_file) or die "$! opening $blast_file for reading";
 233
 234   my $line=0;
 235
 236   while (<$res_fh>) {
 237       $line++;
 238
 239       my @data=split(/\t/, $_);
 240
 241       #separating forward primers from reverse primers using 2 arrays
 242       push (@fprimer_ids , $line) if ($data[0] eq 'FORWARD-PRIMER');
 243       push (@rprimer_ids , $line) if ($data[0] eq 'REVERSE-PRIMER');
 244
 245       # print STDERR "$data[1]\n";
 246
 247       $query_id{$line}=$data[0];
 248       $subject_id{$line}=$data[1];
 249       $identity{$line}=$data[2];
 250       $align_length{$line}=$data[3];
 251       $mismatches{$line}=$data[4];
 252       $gap_openings{$line}=$data[5];
 253       $q_start{$line}=$data[6];
 254       $q_end{$line}=$data[7];
 255       $s_start{$line}=$data[8];
 256       $s_end{$line}=$data[9];
 257       $e_value{$line}=$data[10];
 258       $hit_score{$line}=$data[11];
 259
 260       #finding the orientation of the strand "+" is  5'->3' and "-" 3'->5'
 261       $orientation{$line}= ($s_end{$line}-$s_start{$line} > 0)? '+' : '-';
 262   }
 263
 264   close $res_fh;
 265
 266
 267   #Finding Results
 268   my @pcr_results; #is an array of array references [forward Primer # , reverse primer #, + or - for parent orientation]
 269
 270   foreach my $forward (@fprimer_ids){
 271
 272     # print STDERR "fwd: ".$subject_id{$forward}."\t".$s_start{$forward}."\t".$align_length{$forward}."\t".$mismatches {$forward}."\n";
 273     # print STDERR "params: product length: ".$params->{productLength}."\tmm: ".$params->{allowedMismatches}."\tfwd length: ".$params->{flength}."\trev length: ".$params->{rlength}."\n";
 274
 275     foreach my $reverse (@rprimer_ids){
 276       # print STDERR "rev: ".$subject_id{$reverse}."\t".$s_start{$reverse}."\t".$align_length{$reverse}."\t".$mismatches {$reverse}."\n";
 277
 278                 if ($subject_id{$forward} eq $subject_id{$reverse}    #both on the same subject seq
 279         and  $s_start{$reverse}- $s_start{$forward}<= $params->{productLength} #product Length is within user's choice
 280         and  $mismatches {$forward} <= $params->{allowedMismatches}  #Allowed mismatches by user
 281         and  $mismatches {$reverse} <= $params->{allowedMismatches}
 282         and  $align_length{$forward} == $params->{flength}  #primers match exact length
 283         and  $align_length {$reverse} == $params->{rlength}
 284                 ) {
 285         # print STDERR "\ninside the if\n";
 286         # print STDERR "orientation fwd: ".$orientation{$forward}." s_end fwd: ".$s_end{$forward}."\n";
 287         # print STDERR "orientation rev: ".$orientation{$reverse}." s_end rev: ".$s_end{$reverse}."\n";
 288
 289           #if the product is in the + starnd of parent seq add a + sign in the array
 290
 291           if ( $orientation{$forward} eq '+'     #forward is on the + strand
 292           and  $orientation{$reverse} eq '-'     #reverse is on the - strand b/c its a complement
 293           and  $s_end{$forward} < $s_end{$reverse}  #end of forward located upstream of beginning of reverse
 294           ){
 295             # print STDERR "first if";
 296                 push (@pcr_results , [$forward,$reverse, '+']) ;
 297            }
 298
 299           #if the product is in the - strand of the parent seq add a - sign in the array
 300           elsif ( $orientation{$forward} eq '-'     #forward is on the - strand (complemet here)
 301              and  $orientation{$reverse} eq '+'     #reverse is on the + strand
 302              and  $s_end{$forward} > $s_end{$reverse}  #end of forward located downstream of beginning of reverse
 303             )
 304              {
 305                # print STDERR "scond if\n";
 306                 push (@pcr_results , [$forward,$reverse, '-']);
 307              }
 308              else {
 309                # print STDERR "\nin else\n";
 310              }
 311        }
 312     }#end of the 4each loop
 313
 314   }
 315
 316   ##############################################################################################################################
 317
 318   my $fs = Bio::BLAST::Database->open(full_file_basename => "$basename",);
 319
 320   my $find_seq;
 321   my $find_subseq;
 322   my $find_id;
 323   my $pcr_seq;
 324
 325   my @product_sizes; #for the gel
 326
 327
 328   if (scalar(@pcr_results) ==  0 ){
 329     print STDERR "No PCR Product Found\n";
 330
 331     $pcr_seq = "No PCR Product Found";
 332     my $pcr_img_name = "insilicopcr_no_result.png";
 333
 334     system("cp static/img/$pcr_img_name $gel_img_tempdir");
 335     return ($pcr_seq,$pcr_img_name);
 336   }
 337   else{
 338
 339     foreach my $result (@pcr_results){
 340
 341       #finding parent sequence
 342       $find_id = $subject_id{$result->[0]};
 343       $find_seq = $fs->get_sequence($find_id);
 344       # $find_seq = $bdb->get_sequence($find_id);
 345       # print STDERR "id: $find_id\n";
 346       # print STDERR "seq: $find_seq\n";
 347
 348       # print STDERR "s_sstart1: ".$s_start{$result->[1]}." s_start0: ".$s_start{$result->[0]}." res2: ".$result->[2]."\n";
 349
 350         #finding the pcr result sequence
 351         $find_subseq = $find_seq->subseq($s_start{$result->[0]},$s_start{$result->[1]}) if $result->[2] eq '+';
 352       $find_subseq = $find_seq->subseq($s_start{$result->[1]},$s_start{$result->[0]}) if $result->[2] eq '-';
 353
 354         ######################################################################################
 355
 356         #generating sequence object for the result to be able to find its molecular weight
 357         my $seq_obj = Bio::Seq->new(-seq       => $find_subseq ,
 358                                   -alphabet  => 'dna'
 359       );
 360
 361       my $seq_size = $seq_obj->length;
 362       push (@product_sizes , $seq_size);
 363
 364       #finding the ID of the sequence and adding + sign if it is on the plus strand and - if its on minus strand and some coordinates
 365       $find_id = $find_seq->id();
 366       $find_id .= $result->[2] eq '+' ? ' strand = plus, ' : ' strand = minus, ';
 367       $find_id .= " start = $s_start{$result->[0]}, end = $s_start{$result->[1]}, size = $seq_size bp";
 368       #######################################################################################
 369
 370       #reverse complementing $find_subseq if the orientation is '-'
 371       $pcr_seq = $seq_obj->revcom->seq if $result->[2] eq '-';
 372
 373       $pcr_seq = ">$find_id\n$pcr_seq"
 374     }
 375
 376     ##############################################################################################################################
 377     #Generating a gel of the results
 378
 379     my $gel = Bio::Graphics::Gel->new('pcr' => \@product_sizes,
 380                       -lane_length => 400,
 381                       -lane_width => 100,
 382                       -band_thickness => 2,
 383                     # -min_frag => 1000,
 384                       -font_size => 16,
 385                       -gelcolor => [240,240,240],
 386                       -bandcolor => [170,170,170]
 387                     );
 388
 389     my $gel_img = $gel->img->png;
 390
 391     # create the gel img temp name
 392     my ($fh ,$temp_file) = tempfile( DIR => $gel_img_tempdir, TEMPLATE=>"gel_XXXXXX", SUFFIX => ".png");
 393     print $fh $gel_img;
 394
 395     my $pcr_img_name = basename ($temp_file);
 396
 397     return ($pcr_seq,$pcr_img_name);
 398   }
 399
 400 }
 401
 402 1;