@transform( importExpressionTracks, suffix("_data.import"), ".ttest.expdiff" )
def buildExpressionDifferencesTTest( infile, outfile ):
    '''compute stats on expression differences between each set
    and the appropriate unstimulated set.

    A FDR Q-Value is computed as well for each pairwise comparison.
    '''

    nskipped = 0

    track = infile[:-len("_data.import")]

    outs = open( outfile, "w" )
    outs.write( "cluster_id\tmean1\tvar1\tmean2\tvar2\tdf\tpvalue\tqvalue\tdiff\tdiff95lower\tdiff95upper\tfold\tfold95lower\tfold95upper\n" )

    try:
        control, cluster_ids, treatments, controls = PExpression.getExpressionMeasurements( track )
    except P.PipelineError, msg:
        E.warn("skipped: %s" % msg )
        outs.close()
        return

    E.info( "track=%s, control=%s: computing ttest for %i probesets, %i samples, %i control" % (track, control,
                                                                                                len(cluster_ids),
                                                                                                len(treatments),
                                                                                                len(controls)))

    result, nskipped = Expression.WelchsTTest()( cluster_ids, treatments, controls )

    for s in result:
        # not that values are logged (log2), hence the fold-change
        # is exp( mean1 - mean2 )
        fold = math.pow(2,s.mMean1-s.mMean2)
        outs.write ("%s\t%f\t%f\t%f\t%f\t%5.2f\t%e\t%e\t%f\t%f\t%f\t%f\t%f\t%f\n" % \
                    (s.mProbeset, 
                     s.mMean1,
                     s.mSampleVariance1,
                     s.mMean2, 
                     s.mSampleVariance2,
                     s.mDegreesFreedom,
                     s.mPValue, 
                     s.mQValue, 
                     s.mDifference,
                     s.mDifferenceLower,
                     s.mDifferenceUpper,
                     fold,
                     math.pow(2,s.mDifferenceLower),
                     math.pow(2,s.mDifferenceUpper),
                     ) )
    outs.close()

############################################################
############################################################
############################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD.py/buildExpressionDifferencesTTest/original