resfile_customquantitation.py

Example program illustrating basic usage of ms_customquantitation.

#!/usr/bin/python
 
 
import msparser
import sys
import re
 
def main():
    
# 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.
 
    protein_peptide_mapping = {'ALBU_HUMAN': [ [1, 1], [2, 1], [3, 1] ], 
                           'ALBU_BOVIN': [ [2, 2], [5, 1], [7, 1] ]}
 
    peptide_intensities = {'q1p1': 1000, 'q2p1': 500, 'q3p1': 750, 'q2p2': 750, 'q5p1': 2000, 'q7p1': 1000}
 
    qmethod = msparser.ms_quant_method()
    average = msparser.ms_quant_average()
 
    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)
 
    p = msparser.ms_quant_protocol()
    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.
 
    rationame = "amount"
    h = {'normalisation': 'none', 'outliers': 'none', 'protein_ratio_type': 'average', 'min_num_peptides': 2}
    set_params(qmethod, h)
 
    customquant = msparser.ms_customquantitation(qmethod)
    
    if not customquant.isValid():
        print(customquant.getLastErrorString())
        sys.exit(1)
        
    dump_warnings(customquant.getErrorHandler())
    customquant.clearAllErrors()
    
    dump_quant_method(qmethod, rationame)
    
# Create the protein-peptide mapping
 
    for accession in protein_peptide_mapping:
        peptides = protein_peptide_mapping[accession]
        for peptide in peptides:
            customquant.addPeptideQuantKey(accession, 0, msparser.ms_peptide_quant_key(peptide[0], peptide[1]))
            
# Add peptide intensities
 
    for keys in peptide_intensities:
        match = re.match(r"[a-z](\d+)[a-z](\d+)", keys)
        q = int(match.group(1))
        p = int(match.group(2))
        ratioValue = peptide_intensities[keys]
        quantKey = msparser.ms_peptide_quant_key(q, p)
        ratio = msparser.ms_peptide_quant_ratio(quantKey, rationame, ratioValue) 
        customquant.addPeptideRatio(ratio)
        
    print("\n")
    
# Print it all out immediately:
        
    for accession in protein_peptide_mapping:
        print("%s:" %(accession))
        ratio = customquant.getProteinRatio(accession, 0, rationame)
        
        if ratio.isMissing():
            print("%10s: ---" %(rationame))
        else:
            print("%10s = %10g, stderr = %10g, N = %d, normal-p = %6g, hypo-p = %6g" %(rationame, ratio.getValue(), ratio.getStandardError_only(), ratio.getSampleSize(), ratio.getNormalityPvalue_only(), ratio.getHypothesisPvalue_only()))
            
        print(" %10s    %s" %('Peptide', rationame))
        
        peptides = protein_peptide_mapping[accession]
        for i in ms_range(1, len(peptides)-1):
            peptide = peptides[i]
            for i in range(len(peptide)):
                q = peptide[0]
                p = peptide[1]
                        
            values = ''
            
            ratio = customquant.getPeptideRatio(msparser.ms_peptide_quant_key(q, p), rationame)
            
            if ratio.isMissing():
                values = "---"
            elif ratio.isInfinite():
                values = "###"
            else:
                values = ratio.getValue()
                
            print("    %-9s  %-10s" %("q%d_p%d" %(q, p), str(values)))
            
        print("\n")
 
# ms_range function has been added to make the range supplied consistent with other languages such as Perl and C#
def ms_range(start, stop, step=1):
    i = start
    while i <= stop:
        yield i
        i += step
        
def set_params(qmethod, h):
 
    if qmethod.haveQuality():
        quality = qmethod.getQuality()
    else:
        quality = msparser.ms_quant_quality()
    
    if qmethod.haveNormalisation():
        normalisation = qmethod.getNormalisation()
    else:
        normalisation = msparser.ms_quant_normalisation()
    
    if qmethod.haveOutliers():
        outliers = qmethod.getOutliers()
    else:
        outliers = msparser.ms_quant_outliers()
    
    normalisation.setMethod(h['normalisation'])
    outliers.setMethod(h['outliers'])
    qmethod.setMinNumPeptides(h['min_num_peptides'])
    qmethod.setProteinRatioType(h['protein_ratio_type'])
    quality.setMinPrecursorCharge(1)
    quality.setPepThresholdType('maximum expect')
    quality.setPepThresholdValue('0.05')
 
    qmethod.setQuality(quality)
    qmethod.setNormalisation(normalisation)
    qmethod.setOutliers(outliers)
    
def dump_quant_method(qmethod, rationames):
 
    comps = []
    
    for i in ms_range(1, qmethod.getNumberOfComponents()-1):
        comp = qmethod.getComponentByNumber(i)
        comps.append(comp.getName())
        
    print("Components: %s" %(comps))
    print("Report ratios: %s" %(rationames))
    
    print("Protein ratio type = %s" %(qmethod.getProteinRatioType()))
    print("Min. number of peptides: %d" %(qmethod.getMinNumPeptides()))
    
    if qmethod.haveQuality():
        q = qmethod.getQuality()
        print("Quality: min. precursor charge = %s" %(q.getMinPrecursorCharge()))
        print("QualityL pep. threshold type = %s" %(q.getPepThresholdType()))
        print("Quality: pep. threshold value = %s" %(q.getPepThresholdValue()))
    else:
        print("Quality: no restrictions")
        
    if qmethod.haveNormalisation():
        q = qmethod.getNormalisation()
        print("Normalisation = %s" %(q.getMethod()))
    else:
        print("Normalisation: none")
        
    if qmethod.haveOutliers():
        q = qmethod.getOutliers()
        print("Outliers = %s" %(q.getMethod()))
    else: 
        print("Outliers: none")
        
        
def dump_warnings(err):
    for i in range(1, err.getNumberOfErrors()):
        print("Warning: %s" %(err.getErrorString(i)))
    
 
if __name__ == "__main__":
    main()
    
'''
resfile_customquantitation.py
Will give the following output:
 
Components: []
Report ratios: amount
Protein ratio type = average
Min. number of peptides: 2
Quality: min. precursor charge = 1
QualityL pep. threshold type = maximum expect
Quality: pep. threshold value = 0.05
Normalisation = none
Outliers = none
 
 
ALBU_HUMAN:
    amount =    721.125, stderr =    1.22269, N = 3, normal-p = 0.813228, hypo-p = 0.000932101
    Peptide    amount
    q1_p1      1000.0
    q2_p1      500.0
    q3_p1      750.0
 
 
ALBU_BOVIN:
    amount =    1144.71, stderr =    1.33789, N = 3, normal-p = 0.552542, hypo-p = 0.00170392
    Peptide    amount
    q2_p2      750.0
    q5_p1      2000.0
    q7_p1      1000.0
    
'''