Example program illustrating basic usage of ms_customquantitation.
#!/usr/local/bin/perl ############################################################################## # file: resfile_customquantitation.pl # # Mascot Parser toolkit example code # ############################################################################## # COPYRIGHT NOTICE # # Copyright 1998-2013 Matrix Science Limited All Rights Reserved. # # # ############################################################################## # $Source: parser/examples/test_perl/resfile_customquantitation.pl $ # $Author: villek@matrixscience.com $ # $Date: 2018-07-30 16:23:53 +0100 $ # $Revision: 1b450440f9c97e1e41d0fc6016a27d68951d4532 | MSPARSER_REL_2_8_1-0-gea32989045 $ ############################################################################## use strict; use warnings; ############################################################################## use msparser; # This example uses fictitious data to illustrate how ms_customquantitation can # be used with the Average protocol. In a real situation peptide intensities # and the protein-peptide mapping would come from some external source. my @protein_peptide_mapping = ( { accession => 'ALBU_HUMAN', peptides => [ [1, 1], [2, 1], [3, 1] ] }, { accession => 'ALBU_BOVIN', peptides => [ [2, 2], [5, 1], [7, 1] ] }, ); my %peptide_intensities = ( 'q1p1' => 1_000, 'q2p1' => 500, 'q3p1' => 750, 'q2p2' => 750, 'q5p1' => 2_000, 'q7p1' => 1_000, ); my $qmethod = msparser::ms_quant_method->new(); do { my $average = msparser::ms_quant_average->new; $average->setReferenceAccession('ALBU_HUMAN'); $average->setReferenceAmount('1.0'); # We set the number of top N peptides to 2, so we expect the following # peptide intensities to be used: # ALBU_HUMAN = 1000, 750 # ALBU_BOVIN = 2000, 1000 $average->setNumPeptides(2); my $p = msparser::ms_quant_protocol->new; $p->setAverage($average); $qmethod->setProtocol($p); }; # Since ms_customquantitation does not use the ratio definitions in the # quantitation method, we need not define them at all. Ratios are meaningless # in the Average protocol anyway, so we can use any name we want to denote # the intensity values. my $rationame = "amount"; set_params($qmethod, normalisation => 'none', outliers => 'none', protein_ratio_type => 'average', min_num_peptides => 1, ); my $customquant = msparser::ms_customquantitation->new($qmethod); if (not $customquant->isValid) { print STDERR $customquant->getLastErrorString(), "\n"; exit 1; } dump_warnings($customquant->getErrorHandler()); $customquant->clearAllErrors(); dump_quant_method($qmethod, [$rationame]); # Create the protein-peptide mapping: for (@protein_peptide_mapping) { my ($acc, $peptides) = @$_{qw(accession peptides)}; for (@$peptides) { # Note that the q,p arguments to ms_peptide_quant_key must be integers. # Otherwise the wrong constructor is chosen. Adding to 0 converts the # values to integers. $customquant->addPeptideQuantKey($acc, 0, msparser::ms_peptide_quant_key->new(0+$$_[0], 0+$$_[1])); } } # Add peptide intensities: for (keys %peptide_intensities) { my ($q, $p) = $_ =~ /q(\d+)p(\d+)/; my $ratio = msparser::ms_peptide_quant_ratio->new(msparser::ms_peptide_quant_key->new(0+$q, 0+$p), $rationame, $peptide_intensities{$_}); $customquant->addPeptideRatio($ratio); } print "\n"; # Print it all out immediately: for my $mapping (@protein_peptide_mapping) { printf "%s:\n", $$mapping{'accession'}; my $ratio = $customquant->getProteinRatio($$mapping{'accession'}, 0, $rationame); if ($ratio->isMissing()) { printf "%10s: ---\n", $rationame; } else { printf "%10s = %10g, stderr = %10g, N = %d, normal-p = %6g, hypo-p = %6g\n", $rationame, $ratio->getValue(), $ratio->getStandardError(), $ratio->getSampleSize(), $ratio->getNormalityPvalue(), $ratio->getHypothesisPvalue(); } print "\n"; print join(' ', map { sprintf("%10s", $_) } 'Peptide', $rationame), "\n"; for (@{ $$mapping{'peptides'} }) { my ($q, $p) = @$_; my @values = (); my $ratio = $customquant->getPeptideRatio(msparser::ms_peptide_quant_key->new(0+$q, 0+$p), $rationame); if ($ratio->isMissing()) { push @values, sprintf("%10s", '---'); } elsif ($ratio->isInfinite()) { push @values, sprintf("%10s", '###'); } else { push @values, sprintf("%10g", $ratio->getValue()); } printf "%10s %s\n", sprintf("q%d_p%d", $q, $p), join(' ', @values); } print "\n"; } sub set_params { my ($qmethod, %h) = @_; my $quality; if ($qmethod->haveQuality) { $quality = $qmethod->getQuality(); } else { $quality = msparser::ms_quant_quality->new(); } my $normalisation; if ($qmethod->haveNormalisation) { $normalisation = $qmethod->getNormalisation(); } else { $normalisation = msparser::ms_quant_normalisation->new(); } my $outliers; if ($qmethod->haveOutliers) { $outliers = $qmethod->getOutliers(); } else { $outliers = msparser::ms_quant_outliers->new(); } $normalisation->setMethod($h{'normalisation'} || 'none'); $outliers->setMethod($h{'outliers'} || 'none'); $qmethod->setMinNumPeptides($h{'min_num_peptides'} || 2); $qmethod->setProteinRatioType($h{'protein_ratio_type'} || 'average'); $quality->setMinPrecursorCharge($h{'min_precursor_charge'} || 1); $quality->setPepThresholdType($h{'pep_threshold_type'} || 'maximum expect'); $quality->setPepThresholdValue($h{'pep_threshold_value'} || '0.05'); $qmethod->setQuality($quality); $qmethod->setNormalisation($normalisation); $qmethod->setOutliers($outliers); } sub dump_quant_method { my ($qmethod, $rationames) = @_; do { my @comps = (); for my $i (0 .. $qmethod->getNumberOfComponents()-1) { my $comp = $qmethod->getComponentByNumber($i); $comps[$i] = $comp->getName(); } printf "Components: [%s]\n", join(', ', @comps); }; printf "Report ratios: [%s]\n", join(', ', @$rationames); printf "Protein ratio type = %s\n", $qmethod->getProteinRatioType(); printf "Min. number of peptides = %d\n", $qmethod->getMinNumPeptides(); if ($qmethod->haveQuality()) { my $q = $qmethod->getQuality(); printf "Quality: min. precursor charge = %s\n", $q->getMinPrecursorCharge(); printf "Quality: pep. threshold type = %s\n", $q->getPepThresholdType(); printf "Quality: pep. threshold value = %s\n", $q->getPepThresholdValue(); } else { printf "Quality: no restrictions\n"; } if ($qmethod->haveNormalisation()) { my $q = $qmethod->getNormalisation(); printf "Normalisation = %s\n", $q->getMethod(); } else { printf "Normalisation: none\n"; } if ($qmethod->haveOutliers()) { my $q = $qmethod->getOutliers(); printf "Outliers = %s\n", $q->getMethod(); } else { printf "Outliers: none\n"; } } sub dump_warnings { my ($err) = @_; for my $i (1 .. $err->getNumberOfErrors()) { print 'Warning: ', $err->getErrorString($i), "\n"; } }
Copyright © 2022 Matrix Science Ltd. All Rights Reserved. Generated on Thu Mar 31 2022 01:12:29 |