msparser_full.py

Large example / test program that illustrates most functionality.

#!/usr/bin/python
 
 
import msparser
import sys
import re
import time
 
def main() :
    if len(sys.argv) < 2 :
        print("Must specify results filename as parameter")
        return 1
 
    resfile = msparser.ms_mascotresfilebase.createResfile(sys.argv[1])
 
    if not resfile.isValid() :
        print("Cannot process file '%s':" % sys.argv[1])
        print("Error number: %s" % resfile.getLastError())
        print("Error string: %s" % resfile.getLastErrorString())
        return 1
 
    checkErrorHandler(resfile)
    searchInformation(resfile)
    searchParameters(resfile)
    inputData(resfile)
 
    print("-------------------------------------------------------------")
    print("---   Full Protein summary report                         ---")
    print("-------------------------------------------------------------")
 
    show_results(
        resfile, 
        0, 
        msparser.ms_mascotresults.MSRES_NOFLAG,
        1.0 / 20.0,      # minProbability
        5                # Number of hits
        )
 
    print("-------------------------------------------------------------")
    print("---   Concise Protein summary report                      ---")
    print("-------------------------------------------------------------")
 
    show_results(
        resfile,
        0, 
        msparser.ms_mascotresults.MSRES_GROUP_PROTEINS 
            | msparser.ms_mascotresults.MSRES_SHOW_SUBSETS,
        0,               # minProbability
        3                # Number of hits
        )
 
    if resfile.isMSMS() :
        print("-------------------------------------------------------------")
        print( "---   Peptide summary report                              ---")
        print("-------------------------------------------------------------")
 
        show_results(
            resfile,
            1, # peptide summary
            msparser.ms_mascotresults.MSRES_GROUP_PROTEINS 
                | msparser.ms_mascotresults.MSRES_SHOW_SUBSETS,
            0,               # minProbability
            8,               # number of hits
            0,               # ignoreIonsScoreBelow
            0                # minimum peptide length
            )
    else :
        print("Not an MS-MS results file - cannot show peptide summary report")
 
 
def show_results(resfile, usePeptideSummary, flags, minProteinProb, maxHits, minIonsScore = 0, minPepLenInPepSummary = 0) :
    """
    - file is ms_mascotresfilebase                           
    - usePeptideSummary is true for peptide summary, false for protein
    - flags is the flags for display                      
    - minProteinProb is the minimum protein probability to display 
    - maxHits is the maximum number of hits to display      
    - minIonsScore is the minimum ions score to use              
    - minPepLenInPepSummary is the minPepLenInPepSummary to use           
    """
 
    if usePeptideSummary :
        results = msparser.ms_peptidesummary(
            resfile, flags, minProteinProb, maxHits, "", minIonsScore, minPepLenInPepSummary
            )
    else :
        results = msparser.ms_proteinsummary(
            resfile, flags, minProteinProb, maxHits
            )
 
    if resfile.getLastError() :
        print("Error : %s" % resfile.getLastErrorString())
        return
    
    hit  = 1
    prot = results.getHit(hit)
 
    while prot :
        accession   = prot.getAccession()
        description = results.getProteinDescription(accession)
        mass        = results.getProteinMass(accession)
 
        print("Protein Hit %d" % hit)
        print("===================")
        print("Accession   : %s" % accession)
        print("Description : %s" % description)
        print("Score       : %s" % prot.getScore())
        print("Mass        : %s" % mass)
        print("Frame       : %s" % prot.getFrame())
        print("Coverage    : %s" % prot.getCoverage())
        print("RMS error   : %s" % prot.getRMSDeltas(results))
        print("Peptides    : %s" % prot.getNumDisplayPeptides())
 
        # Each protein has a number of peptides that matched - list them:
        num_peps = prot.getNumPeptides()
 
        for i in range(1, 1+ num_peps) :
            query = prot.getPeptideQuery(i)
            p     = prot.getPeptideP(i)
 
            isDupSameQuery = prot.getPeptideDuplicate(i) != msparser.ms_protein.DUPE_DuplicateSameQuery 
            if p != -1 and query != -1 and not isDupSameQuery :
                pep = results.getPeptide(query, p)
                if not pep:
                    continue
 
                displayPeptideInfo(
                    0, pep, results, 
                    prot.getPeptideDuplicate(i) == msparser.ms_protein.DUPE_Duplicate,
                    prot.getPeptideIsBold(i),
                    prot.getPeptideShowCheckbox(i)
                    )
 
        # Now display list of all proteins that contained subsets or and
        # identical list of the matching peptides
 
        if flags & msparser.ms_mascotresults.MSRES_GROUP_PROTEINS :
            print("Proteins matching the same set of peptides:")
 
            i = 1
            similar_prot = results.getNextSimilarProtein(hit, 1)
            while similar_prot :
                accession = prot.getAccession()
                print(accession, "Total score:", similar_prot.getScore(),)
                print(" Peptides matched:", similar_prot.getNumDisplayPeptides())
                i += 1
                similar_prot = results.getNextSimilarProtein(hit, i)
            
            if flags & msparser.ms_mascotresults.MSRES_SHOW_SUBSETS :
                print("Proteins matching a subset of these peptides:")
 
                i = 1
                subset_prot = results.getNextSubsetProtein(hit, 1)
                while subset_prot :
                    accession = subset_prot.getAccession()
                    print(accession,  "Total score:", subset_prot.getScore(),)
                    print(" Peptides matched:", subset_prot.getNumDisplayPeptides(), "")
                    i += 1
                    subset_prot = results.getNextSubsetProtein(hit, i)
                
        print(" ")
 
        hit += 1
        prot = results.getHit(hit)
 
 
    results.createUnassignedList(msparser.ms_mascotresults.SCORE)
 
    if results.getNumberOfUnassigned() :
        print("Unassigned list")
        print("---------------")
 
        for u in range(1, 1 + results.getNumberOfUnassigned()) :
            pep = results.getUnassigned(u)
            displayPeptideInfo(0, pep, results, 0, 1, 1)
        
 
    if usePeptideSummary :
        print(" ")
        displayYellowPopupInfo(results, 97)
 
 
def displayYellowPopupInfo(results, q) :
    """
    Shows the equivalent of the yellow popup box for given query
    - results is the results object
    - q is the query number
    """
 
    fmt = "%5s %5s %9s %7s %7s"
    print(fmt % ("Score", "Delta", "Hit", "Protein", "Peptide"))
 
    for p in range(1, 11) :
        pep = results.getPeptide(q, p)
        if not pep: continue
 
        seq = pep.getPeptideStr()
        if not seq: continue
 
        tmp = results.getProteinsWithThisPepMatch(q, p)
 
        (hit, protein) = ('', '')
 
        if tmp :
            match = re.search('(\d+):([^ ]*)[ ]*(.*)', tmp)
            hit, protein, morethan = match.group(1), match.group(2), match.group(3)
 
            if morethan :
                hit += "+"
 
        print(fmt % (pep.getIonsScore(), pep.getDelta(), hit, protein, seq))
 
    p = 1
    print("Accessions that matched query %s rank %s :- %s" % (q, p, results.getProteinsWithThisPepMatch(q, p)))
 
 
def displayPeptideInfo(showFullDetails, p, results, isDuplicate, isBold, showCheckBox) :
    q = p.getQuery()
 
    if not showFullDetails :
        fmt = "%2s %4s %4d %11f %4d(%4d) %-20s %s%3.2f%s %3s"
 
        cb, bold = "--", "dim"
        if showCheckBox : cb = "CB" 
        if isBold : bold = "BOLD" 
        
        paren1, paren2 = "", ""
        if isDuplicate : paren1 = "(" ; paren2 = ")" 
        
        print(fmt % (
            cb,
            bold,
            q,
            p.getObserved(),
            p.getRank(),
            p.getPrettyRank(),
            p.getPeptideStr(),
            paren1,
            p.getIonsScore(),
            paren2,
            results.getReadableVarMods(q, p.getRank())
            ))
 
        return
    
 
    print("Peptide hit")
 
    if p.getAnyMatch() :
        fmt = "    %-12s: %s"
        print(fmt % ('Query'       , q))
        print(fmt % ('Rank'        , p.getRank()))
        print(fmt % ('Matched'     , p.getAnyMatch()))
        print(fmt % ('missedCleave', p.getMissedCleavages()))
        print(fmt % ('mrCalc'      , p.getMrCalc()))
        print(fmt % ('delta'       , p.getDelta()))
        print(fmt % ('observed'    , p.getObserved()))
        print(fmt % ('charge'      , p.getCharge()))
        print(fmt % ('mrExp'       , p.getMrExperimental()))
        print(fmt % ('ionsMatched' , p.getNumIonsMatched()))
        print(fmt % ('peptideStr'  , p.getPeptideStr()))
        print(fmt % ('peaksUsed1'  , p.getPeaksUsedFromIons1()))
        print(fmt % ('varModsStr'  , p.getVarModsStr()))
        print(fmt % ('readable mod', results.getReadableVarMods(q, p.getRank)))
        print(fmt % ('ionsScore'   , p.getIonsScore()))
        print(fmt % ('seriesUsedS' , p.getSeriesUsedStr()))
        print(fmt % ('peaksUsed2'  , p.getPeaksUsedFromIons2()))
        print(fmt % ('peaksUsed3'  , p.getPeaksUsedFromIons3()))
        print(fmt % ('idth, hth, p', ', '.join(
            results.getPeptideIdentityThreshold(q, 20),
            results.getHomologyThreshold(q, 20),
            results.getProbOfPepBeingRandomMatch(p.getIonsScore(), q)
            )))
        print(" ")
    else :
        print("    No match")
 
 
def checkErrorHandler(resfile) :
    """
    Calls a couple of functions with invalid arguments, then prints
    their error messages.
    """
 
    print("Testing the error handling...            ")
    print("=========================================")
 
    numQueries = resfile.getNumQueries()
    resfile.getObservedCharge(numQueries + 40);  # Should fail
 
    print("Error number: %s" % resfile.getLastError())
    print("Error string: %s" % resfile.getLastErrorString())
 
    resfile.clearAllErrors()
    print("Cleared all errors - should have no errors left: %s errors left" %resfile.getNumberOfErrors())
    
    for x in range(1, 21) :
        resfile.getObservedCharge(numQueries + x);  # Should fail
    
    # Now, the best way, print out all errors.
    print("More errors added - there are now %s errors" % resfile.getNumberOfErrors())
 
    for i in range(1, 1 + resfile.getNumberOfErrors()) :
        print("Error number: %s : %s" % (resfile.getErrorNumber(i), resfile.getErrorString(i)))
           
 
    print(" ")
    resfile.clearAllErrors()
 
 
def searchInformation(resfile) :
    """
    Display parameters from the ms_mascotresfilebase object. The functions
    anyPMF, anyMSMS, anySQ should normally be used in preference to isPMF etc
    because some people submit MSMS though the sequence query window etc.
    """
 
    fmt = "%-20s: %s"
 
    print("Search information from ms_mascotresfilebase")
    print("========================================")
    print(fmt % ("Number of queries"   , resfile.getNumQueries()))
    print(fmt % ("Number of sequences" , resfile.getNumSeqs()))
    print(fmt % ("Sequences after tax" , resfile.getNumSeqsAfterTax()))
    print(fmt % ("Number of residues"  , resfile.getNumResidues()))
    print(fmt % ("Execution time"      , resfile.getExecTime()))
    print(fmt % ("Date (seconds)"      , resfile.getDate()))
 
    date = time.localtime(resfile.getDate())
    Wdays = "Mon Tue Wed Thu Fri Sat Sun".split(" ")
    Mons = "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec".split(" ")
 
    print(fmt %( "Date", "%s %s %d %02d:%02d:%02d %d" % (
        Wdays[date.tm_wday],
        Mons[date.tm_mon - 1],
        date.tm_mday,
        date.tm_hour,
        date.tm_min,
        date.tm_sec,
        date.tm_year
        )))
 
    print(fmt % ("Mascot version"       , resfile.getMascotVer()))
    print(fmt % ("Fasta version"        , resfile.getFastaVer()))
    print(fmt % ("Is PMF?"              , resfile.isPMF()))
    print(fmt % ("Is MSMS?"             , resfile.isMSMS()))
    print(fmt % ("Is SQ?"               , resfile.isSQ()))
    print(fmt % ("Is Error tolerant"    , resfile.isErrorTolerant()))
    print(fmt % ("Any PMF?"             , resfile.anyPMF()))
    print(fmt % ("Any MSMS?"            , resfile.anyMSMS()))
    print(fmt % ("Any SQ?"              , resfile.anySQ()))
    print(fmt % ("Any peptide matches"  , resfile.anyFastaMatches()))
 
    print(" ")
 
 
 
def searchParameters(resfile) :
    """
    Display parameters from the ms_searchparams object.
    The values come from the parameters and the masses sections of the file.
    """
 
    params = resfile.params()
    fmt = "%-20s: %s"
 
    print("Search parameters from ms_searchparams")
    print("=========================================")
    print(fmt % ("License"             , params.getLICENSE()))
    print(fmt % ("Search title"        , params.getCOM()))
    print(fmt % ("SEG mass"            , params.getSEG()))
    print(fmt % ("Peptide tol"         , params.getTOL()))
    print(fmt % ("Peptide tol units"   , params.getTOLU()))
    print(fmt % ("Fragment tol"        , params.getITOL()))
    print(fmt % ("Fragment tol units"  , params.getITOLU()))
    print(fmt % ("Missed cleavages"    , params.getPFA()))
    print(fmt % ("Database"            , params.getDB()))
    print(fmt % ("Static mods"         , params.getMODS()))
    print(fmt % ("Average/monoisotopic", params.getMASS()))
    print(fmt % ("Enzyme"              , params.getCLE()))
    print(fmt % ("Raw data file name"  , params.getFILENAME()))
    print(fmt % ("Input data"          , params.getQUE()))
    print(fmt % ("Type of search"      , params.getSEARCH()))
    print(fmt % ("User name"           , params.getUSERNAME()))
    print(fmt % ("User email"          , params.getUSEREMAIL()))
    print(fmt % ("Charge state"        , params.getCHARGE()))
    print(fmt % ("Repeat search file"  , params.getINTERMEDIATE()))
    print(fmt % ("Num hits to display" , params.getREPORT()))
    print(fmt % ("Show overview"       , params.getOVERVIEW()))
    print(fmt % ("Data file format"    , params.getFORMAT()))
    print(fmt % ("Form version"        , params.getFORMVER()))
    print(fmt % ("Variable mods"       , params.getIT_MODS()))
 
    for i in range(12) :
        print(fmt % ( "User%02d" % i, params.getUSERField(i)))
    
    print(fmt % ("Precursor mass"      , params.getPRECURSOR()))
    print(fmt % ("Taxonomy filter"     , params.getTAXONOMY()))
    print(fmt % ("Type of report"      , params.getREPTYPE()))
    print(fmt % ("Accessions to search", params.getACCESSION()))
    print(fmt % ("Subcluster used"     , params.getSUBCLUSTER()))
    print(fmt % ("ICAT search?"        , params.getICAT()))
    print(fmt % ("Instrument type"     , params.getINSTRUMENT()))
    print(fmt % ("Error tolerant?"     , params.getERRORTOLERANT()))
    print(fmt % ("Rules (ions series)" , params.getRULES()))
 
    for ch in range(ord('A'), 1 + ord('Z')) :
        print(fmt % ("Residue " + chr(ch), params.getResidueMass(chr(ch))))
    
    print(fmt % ("C terminus mass" , params.getCTermMass()))
    print(fmt % ("N terminus mass" , params.getNTermMass()))
    print(fmt % ("Mass of hydrogen", params.getHydrogenMass()))
    print(fmt % ("Mass of oxygen"  , params.getOxygenMass()))
    print(fmt % ("Mass of carbon"  , params.getCarbonMass()))
    print(fmt % ("Mass of nitrogen", params.getNitrogenMass()))
    print(fmt % ("Mass of electron", params.getElectronMass()))
 
    i = 1
    while params.getVarModsName(i) :
        print(fmt % ("Variable mod name"   , params.getVarModsName(i)))
        print(fmt % ("Variable mod delta"  , params.getVarModsDelta(i)))
        print(fmt % ("Variable mod neutral", params.getVarModsNeutralLoss(i)))
        i += 1
 
    print(" ")
 
 
def inputData(resfile) :
    # display input data via inputquery get functions
    for i in range(1, 2) :
        # replace range with (1, 1 + resfile.getNumQueries()) to get all input
        # data.
 
        print("Input data for query number %d" % i)
        print("=========================================")
 
        q = msparser.ms_inputquery(resfile, i)
        fmt = "    %-16s: %s"
 
        print(fmt % ("title"       , q.getStringTitle(True)))
        print(fmt % ("mass_min"    , q.getMassMin()))
        print(fmt % ("mass_max"    , q.getMassMax()))
        print(fmt % ("int_min"     , q.getIntMin()))
        print(fmt % ("int_max"     , q.getIntMax()))
        print(fmt % ("num_vals"    , q.getNumVals()))
        print(fmt % ("num_used1"   , q.getNumUsed()))
        print(fmt % ("ions1"       , q.getStringIons1()))
        print(fmt % ("ions2"       , q.getStringIons2()))
        print(fmt % ("ions3"       , q.getStringIons3()))
        print(fmt % ("peptol"      , q.getPepTol()))
        print(fmt % ("peptol units", q.getPepTolUnits()))
        print(fmt % ("peptol str"  , q.getPepTolString()))
        print(fmt % ("repeat srch" , resfile.getRepeatSearchString(i)))
 
        num_peaks = q.getNumberOfPeaks(1)
        for j in range (1, 1+ num_peaks) :
            print("%f, %f" % (q.getPeakMass(1, j), q.getPeakIntensity(1, j)))
                
    print(" ")
 
 
if __name__ == "__main__" :
    sys.exit(main())