view release on metacpan or  search on metacpan

scripts/peptSpectra/getIonIntensities.pl  view on Meta::CPAN

      $massIndex = $itemIndex;
    }
    elsif (($_{type} eq 'intensity') || ($_{type} eq 'height')){
      $intensityIndex = $itemIndex;
    }
    $itemIndex++;
  }

  undef($curChar);

} # StartTag


sub EndTag
{
  my($p, $el)= @_;
  if ($el eq 'idi:sequence'){
    $peptide = $curChar;
  }
  elsif ($el eq 'idi:modif'){
    $modif = $curChar;
  }
  elsif ($el eq 'idi:charge'){
    $charge = $curChar;
  }
  elsif ($el eq 'ple:peaks'){
    $peaks = $curChar;
  }
  elsif (($el eq 'idi:OneIdentification') && (!defined($imposedCharge) || ($charge == $imposedCharge))){

    # Builds the peak list
    my @lines = split(/\n/, $peaks);
    my (@peaks, $maxIntensity);
    foreach my $line (@lines){
      $line =~ s/\r//g;
      my @part = split(/\s+/, $line);
      if ($part[$massIndex] > 0.0){
	push(@peaks, [$part[$massIndex], $part[$intensityIndex]]);
	$maxIntensity = $part[$intensityIndex] if ($part[$intensityIndex] > $maxIntensity);
      }
    }

    if (defined($maxIntensity)){
      # Matches with theoretical masses
      my %spectrum;
      if ($matchSel eq 'closest'){
	matchSpectrumClosest(pept=>$peptide, modif=>$modif, spectrum=>\%spectrum, expSpectrum=>\@peaks, fragTypes=>\@fragTypes);
      }
      elsif ($matchSel eq 'greedy'){
	matchSpectrumGreedy(pept=>$peptide, modif=>$modif, spectrum=>\%spectrum, expSpectrum=>\@peaks, fragTypes=>\@fragTypes, order=>\@fragTypes, tol=>$tol);
      }
      elsif ($matchSel eq 'mostintense'){
	matchSpectrumGreedy(pept=>$peptide, modif=>$modif, spectrum=>\%spectrum, expSpectrum=>\@peaks, fragTypes=>\@fragTypes, tol=>$tol);
      }
      else{
	CORE::die("Unknown match type [$matchSel]");
      }

      if (defined($withPlots)){
	my $msms = new InSilicoSpectro::InSilico::MSMSOutput(spectrum=>\%spectrum, prec=>2, modifLvl=>1, expSpectrum=>\@peaks, intSel=>'order', tol=>$tol, minTol=>$minTol);
	$msms->plotSpectrumMatch(fname=>"$peptide-$$-$nPlot", format=>'png', fontChoice=>'default:Large', changeColModifAA=>1, legend=>'right', plotIntern=>1);
	$nPlot++;
      }

      # Normalizes intensities
      my %normInt;
      normalizeIntensities($intSel, \@peaks, \%normInt);

      # Extracts statistics
      my $len = length($peptide);
      foreach my $frag (keys(%{$spectrum{mass}{term}})){
	for (my $i = 0; $i < @{$spectrum{ionType}{$frag}}; $i++){
	  for (my $j = $i*$len; $j < ($i+1)*$len; $j++){
	    if (defined($spectrum{mass}{term}{$frag}[$j]) && defined($spectrum{match}{term}{$frag}[$j])){
	      my $theoMass = $spectrum{mass}{term}{$frag}[$j];
	      my $expMass = $spectrum{match}{term}{$frag}[$j][0];
	      if ((abs($theoMass-$expMass)/($theoMass+$expMass)*2.0e6 <= $tol) || (abs($theoMass-$expMass) <= $minTol)){
		push(@{$intensStat{$spectrum{ionType}{$frag}[$i]}}, $normInt{$spectrum{match}{term}{$frag}[$j][0]});
	      }
	    }
	  }
	}
      }
      foreach my $frag (keys(%{$spectrum{match}{intern}})){
	foreach my $aa (keys(%{$spectrum{match}{intern}{$frag}})){
	  if (defined($spectrum{match}{intern}{$frag}{$aa})){
	    my $theoMass = $spectrum{mass}{intern}{$frag}{$aa};
	    my $expMass = $spectrum{match}{intern}{$frag}{$aa}[0];
	    if ((abs($theoMass-$expMass)/($theoMass+$expMass)*2.0e6 <= $tol) || (abs($theoMass-$expMass) <= $minTol)){
	      push(@{$intensStat{"$frag($aa)"}}, $normInt{$spectrum{match}{intern}{$frag}{$aa}[0]});
	    }
	  }
	}
      }
    }
  }

} # EndTag