resfile_ms2quant_to_custom.py

Example program converting ms_ms2quantitation to ms_customquantitation.

#!/usr/bin/python
 
 
import msparser
import sys
 
# If no arguments supplied, print help and exit
 
if len(sys.argv) != 5:
    print("Usage: %s <quantitation schema 1 filepath> <quantitation schema 2 filepath> <unimod schema filepath> <results file>" % sys.argv[0])
    print("The location of schema files should be e.g.")
    print("  ../html/xmlns/schema/quantitation_1/quantitation_1.xsd")
    print("  ../html/xmlns/schema/quantitation_2/quantitation_2.xsd")
    print("  ../html/xmlns/schema/unimod_2/unimod_2.xsd")
    print("if running from the Mascot cgi directory.")
    sys.exit(1)
    
quant_schema1_filepath = sys.argv[1]
quant_schema2_filepath = sys.argv[2]
unimod_schema_filepath = sys.argv[3]
 
resfile = msparser.ms_mascotresfilebase.createResfile(sys.argv[4], 1)
 
if not resfile.isValid():
    print(resfile.getLastErrorString())
    sys.exit(1)
 
# 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 dump_warnings(err):
    
    for i in ms_range(1, err.getNumberOfErrors()):
        print('Warning: %s' %(err.getErrorString(i)))
 
# And print out related warnings, if any (we know these are only warnings,
# because isValid() returned true earlier):
 
dump_warnings(resfile.getErrorHandler())
resfile.clearAllErrors()
 
resfile.setXMLschemaFilePath(msparser.ms_mascotresfilebase.XML_SCHEMA_QUANTITATION, "http://www.matrixscience.com/xmlns/schema/quantitation_2 " + quant_schema2_filepath + " http://www.matrixscience.com/xmlns/schema/quantitation_1 " + quant_schema1_filepath)
resfile.setXMLschemaFilePath(msparser.ms_mascotresfilebase.XML_SCHEMA_UNIMOD, "http://www.unimod.org/xmlns/schema/unimod_2 " + unimod_schema_filepath)
 
qmethod = msparser.ms_quant_method()
 
if not resfile.getQuantitationMethod(qmethod):
    print(resfile.getLastErrorString())
    sys.exit(1)
 
datfile = msparser.ms_datfile()
undef, flags, minprob, maxhits, iisb, minpeplen, use_pepsum, flags2 = resfile.get_ms_mascotresults_params(datfile.getMascotOptions())
 
if not use_pepsum:
    print("Results File cannot be opened as a peptide summary")
    sys.exit(1)
   
pepsum = msparser.ms_peptidesummary(resfile, flags, minprob, maxhits, '', iisb, minpeplen, '', flags2)
 
if not resfile.isValid():
    print(resfile.getLastErrorString)
    sys.exit(1)
 
dump_warnings(resfile.getErrorHandler())
resfile.clearAllErrors()
 
# Quantitate the file using ms_ms2quantitation:
 
ms2quant = msparser.ms_ms2quantitation(pepsum, qmethod)   
 
if not ms2quant.isValid():
    print(ms2quant.getLastErrorString())
    sys.exit(1)
 
if ms2quant.hasIntensityNormalisation():
    ms2quant.normaliseIntensities()
elif ms2quant.hasPeptideRatioNormalisation():
    ms2quant.normalisePeptideRatios()
 
dump_warnings(ms2quant.getErrorHandler())
ms2quant.clearAllErrors()
 
# Then import the data into ms_customquantitation:
customquant = msparser.ms_ms1quantitation(pepsum, qmethod)
 
if not customquant.isValid():
    print(customquant.getLastErrorString())
    sys.exit(1)
 
err = customquant.getErrorHandler()
# And print out related warnings, if any (we know these are only warnings,
# because isValid() returned true earlier):
 
for i in range(1, err.getNumberOfErrors()):
    print("Warning: %s" %(err.getErrorString(i)))
    
customquant.clearAllErrors()
 
# Refresh the quantitation method object, in case customquant altered
# some of the parameters during conversion:
 
qmethod = customquant.getQuantitationMethod()
 
# Now we can inspect peptide and protein data in customquant.
# But first, dump quantitation method parameters to see what we have:
 
comps = []
 
for i in ms_range(0, qmethod.getNumberOfComponents()-1):
    comp = qmethod.getComponentByNumber(i)
    comps.append(comp.getName())
    
print("Components: %s" %(comps))
 
rationames = []
 
for i in ms_range(0, qmethod.getNumberOfReportRatios()-1):
    ratioByNumber = qmethod.getReportRatioByNumber(i)
    rationames.append(ratioByNumber.getName())
    
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("Quality: 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()))
    if q.haveProteins():
        p = q.getProteins()
        print("getProteins")
        for i in ms_range(0, p.getNumberOfProteins()-1):
            proteins = p.getProteins(i)
            protein = proteins.getAccession()
            print(" - Accessions = %s" %(protein))
 
    if q.havePeptides():
        p = q.getPeptides()
        peptides = p.getPeptides()
        for i in ms_range(0, p.getNumberOfPeptides()-1):
            peptides = p.getPeptides(i)
            peptide = peptides.getSequence()
            print(" - Sequences = %s" %(peptide))
            
    component_bases = []
    for i in ms_range(0, qmethod.getNumberOfComponents()-1):
        gcbn = qmethod.getComponentByNumber(i)
        base = ms2quant.getIntensityNormalisationBase(gcbn.getName())
        component_bases.append(base)
        
    print("Intensity normalisation constants: %s" %(component_bases))
    print("Ration normalisation constants: %s" %(rationames))
    
else:
    print("Normalisation: none")
    
if qmethod.haveOutliers():
    q = qmethod.getOutliers()
    print("Outliers = %s" %(q.getMethod()))
else:
    print("Outliers: none")
    
print("\n")
 
# Collect proteins we're interested in:
 
proteins = []
 
for i in ms_range(1, pepsum.getNumberOfHits()):
    hit = pepsum.getHit(i)
    proteins.append(hit)
    j = 0
    protein = pepsum.getNextFamilyProtein(i, j)
    while protein !=None:
        proteins.append(protein)
        j += 1
        
for protein in proteins:
    print("%d::%s:" %(protein.getDB(), protein.getAccession()))
    
    for rationame in rationames:
        ratio = customquant.getProteinRatio(protein.getAccession(), protein.getDB(), rationame)
        
        if ratio.isMissing():
            print("%10s: ---" %(rationame))
        else:
            print("%10s = %10g, stderr = %10g, N = %d, normal-p = %6g, hyp-p = %6g" %(rationame, ratio.getValue(), ratio.getStandardError_only(), ratio.getSampleSize(), ratio.getNormalityPvalue_only(), ratio.getHypothesisPvalue_only()))
            
    print("  %10s %s" %('Peptide:', rationame))
    
    for i in ms_range(1, protein.getNumPeptides()):
        if protein.getPeptideDuplicate(i) == msparser.ms_protein.DUPE_DuplicateSameQuery:
            continue
        q = protein.getPeptideQuery(i)
        p = protein.getPeptideP(i)
        
        peptide = pepsum.getPeptide(q, p)
        
        if not peptide:
            continue
        
        values = ''
        for rationame in rationames:
            # Note that the q,p arguments to ms_peptide_quant_key must be
            # integers. Otherwise the wrong constructor is chosen. Here they
            # are integers, as they are return values from Parser functions
            # that return ints.
            quantKey = msparser.ms_peptide_quant_key(q, p)
            ratio = ms2quant.getPeptideRatio(quantKey, rationame)
            
            if ratio.isMissing():
                values = "---"
                continue
                
            if ratio.isInfinite():
                values = "###"
            else:
                values = ratio.getValue()
                
        print("    %-9s  %-10s" %("q%d_p%d" %(q, p), str(values)))
 
    print("\n")