resfile_customquantitation.pl

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_3_0_0-2024-09-24-0-g93ebaeb4f4 $ 
##############################################################################
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";
    }
}