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

    nskipped = 0

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

    target_path = os.path.join( os.path.abspath( PARAMS["exportdir"] ), "SAM" )
    if not os.path.exists( target_path): 
        try:
            os.makedirs( target_path )
        except OSError: 
            pass

    outs = open( outfile, "w" )
    outs.write( "cluster_id\tmean1\tstd1\tmean2\tstd2\tpvalue\tqvalue\tdiff\tfold\tcalled\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 SAM for %i probesets, %i samples, %i control" % (track, control,
                                                                                              len(cluster_ids),
                                                                                              len(treatments),
                                                                                              len(controls)))

    genes, summary, cutoffs = Expression.SAM()( cluster_ids, treatments, controls,
                                                pattern = os.path.join(target_path, outfile + "%s"),
                                                fdr = float(PARAMS["expression_sam_fdr"]),
                                                ngenes = float(PARAMS["expression_sam_ngenes"]),
                                                ndelta = float(PARAMS["expression_sam_ndelta"]),
                                                npermutations = PARAMS["expression_sam_permutations"],
                                                method = PARAMS["expression_sam_method"],
                                                use_excel_sam = False )

    if summary == None:
        E.warn( "no cutoff when running sam for %s" % infile )

    if cutoffs:
        logs = open( outfile + ".cutoffs", "w" )
        logs.write("\t".join( cutoffs[0]._fields) + "\n" )
        for x in cutoffs:
            logs.write( "\t".join(map(str,x) ) + "\n" )
        logs.close()
    for s in genes:
        outs.write ("%s\t%f\t%f\t%f\t%f\t%e\t%f\t%f\t%f\t%i\n" % \
                        (s.probeset, 
                         s.mean1,
                         s.stddev1,
                         s.mean2,
                         s.stddev2,
                         s.pvalue, 
                         s.qvalue, 
                         s.difference,
                         s.fold,
                         s.called,
                         ) )
    outs.close()
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD.py/buildExpressionDifferencesSAM/original