display_quant_time_align.py

Example program illustrating usage of ms_ms1quant_time_align.

#!/usr/bin/python
 
 
from __future__ import print_function
 
import sys
from datetime import datetime
import msparser
 
 
def checkErrorHandler(errs) :
   for i in range(1, 1 + errs.getNumberOfErrors()) :
        print("Error number: %s : %s" % (errs.getErrorNumber(i), errs.getErrorString(i)))
 
 
# If no arguments supplied, print help and exit
 
if len(sys.argv) < 3:
    print("Usage: %s <time_alignment.xml> <schema_directory>" % sys.argv[0])
    sys.exit(1)
 
timeAlignmentXmlFile = sys.argv[1]
schemaDirectory = sys.argv[2]
start_time = datetime.now()
 
time_align = msparser.ms_ms1quant_time_align()
if time_align.loadXmlFile(timeAlignmentXmlFile, schemaDirectory) :
    elapsed = datetime.now() - start_time
    print("File loaded: %s in %d seconds" % (timeAlignmentXmlFile, elapsed.seconds))
    print("")
    
    data = time_align.getFinalResults()
    num_sub_projects = len(data)
    fraction_numbers = time_align.getSubProjectToFractionMap()
    
    # Print project details
    print("Prj\tFractn\tmeanErrsRaw\tmeanErrsAlig\tstdevErrsRaw\tstdevErrsAlig\tpearsonRaw\tpearsonAlign")
    for prj in range(num_sub_projects):
        fraction = fraction_numbers[prj]
        (success, meanErrorsRaw, meanErrorsAligned, stdevErrorsRaw, stdevErrorsAligned, pearsonCoefficientRaw,  pearsonCoefficientAligned) = time_align.getEvaluation(1 + prj)
        print("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t" % (1+prj, fraction, meanErrorsRaw, meanErrorsAligned, stdevErrorsRaw, stdevErrorsAligned, pearsonCoefficientRaw,  pearsonCoefficientAligned))
    print("")
    
    # Print fraction limits
    print("Fract\tminimum m/z\tmaximum mz\tminimum rt\tmaximum rt\tminbin\tmaxbin\tminrt\tmaxrt\tnum q\t#q>thr\tprojects")
    fraction_to_prj = time_align.getFractionToSubProjectMap()
    for fraction,proj in fraction_to_prj.items():
        sub_proj_str = ""
        for x in proj:
            if sub_proj_str:
                sub_proj_str += ", "
            sub_proj_str = sub_proj_str + str(1+x)
        limits = time_align.getCombinedLimits(fraction)
        print("%d\t%f\t%f\t%f\t%f\t%d\t%d\t%d\t%d\t%d\t%d\t%s" % (fraction, limits.getMinMZ(), limits.getMaxMZ(), limits.getMinRT(), limits.getMaxRT(),limits.getCombinedMinBinNum(), limits.getCombinedMaxBinNum(), limits.getCombinedMinRT(), limits.getCombinedMaxRT(), limits.getNumQueries(), limits.getNumQueriesAboveThreshold(), sub_proj_str))
    print("")
 
 
    showAllNonZero = 0
    if (showAllNonZero) :
        for prj in range(num_sub_projects):
            fraction = fraction_numbers[prj]
            limits = time_align.getCombinedLimits(fraction)
            print("Sub project: %d, fraction: %d" % (prj, fraction))
            num_bins = len(data[prj])
            print("  Number of bins: %d" % (num_bins))
            for bin_num in range(num_bins) :
                print("    Bin: %d" % (bin_num))
                num_rts = len(data[prj][bin_num])
                print("      Num rts: %d" % (num_rts))
                for rt in range(num_rts) :
    #               print("data[%d][%d][%d]: %f" % (prj, bin_num, rt, data[prj][bin_num][rt]))
                    if (data[prj][bin_num][rt] != 0) :
                        print("SubProj: %d, bin: %d, rt: %d, offset: %f" % (1+prj, bin_num + limits.getCombinedMinBinNum(), rt + limits.getCombinedMinRT(), data[prj][bin_num][rt]))
    else :
        mid_project = int(num_sub_projects /2) -1
        mid_bin = int(len(data[mid_project]) /2)
        num_rts = len(data[mid_project][mid_bin])
        fraction = fraction_numbers[mid_project]
        limits = time_align.getCombinedLimits(fraction)
        numNonZero = 0
        elapsed = datetime.now() - start_time
        print("Calculated mid points for: %s, elapsed time: %d seconds" % (timeAlignmentXmlFile, elapsed.seconds))
        print("Showing first 50 non zero deltas from sub project: %d, bin number: %d, Num rts: %d, Min rt: %d" % (1 + mid_project, mid_bin + limits.getCombinedMinBinNum(),num_rts, limits.getCombinedMinRT()))
        for rt in range(num_rts) :
            if (data[mid_project][mid_bin][rt] != 0) :
                if (numNonZero < 50) :
                    numNonZero = 1 + numNonZero
                    print("SubProj: %d, bin: %d, rt: %d, offset: %f" % (1+mid_project, mid_bin + limits.getCombinedMinBinNum(), rt + limits.getCombinedMinRT(), data[mid_project][mid_bin][rt]))
 
        print("")
        print("Showing all the feature widths for sub project: %d, bin number: %d" % (1 + mid_project, mid_bin + limits.getCombinedMinBinNum()))
        feature_widths = time_align.getFeatureWidths()
        features = feature_widths[mid_project][mid_bin]
        for rt, width in features.items() :
            print("%f\t%f" % (rt, width))
 
    print("Time: %s" % (str(datetime.now())))
    
checkErrorHandler(time_align.getErrorHandler())