def compute_statistics(out_dir, gold_prefix, source_prefix, same_samples,
                       use_sge, final_out_prefix):
    """Compute the statistics."""
    # Now, creating a temporary directory
    if not os.path.isdir(out_dir):
        try:
            os.mkdir(out_dir)
        except OSError:
            msg = "{}: file exists".format(out_dir)
            raise ProgramError(msg)

    # The out prefix
    out_prefix = os.path.join(out_dir, "tmp")

    # Subsetting the files
    logger.info("  - Subsetting the files")
    for k, (gold_sample, source_sample) in enumerate(same_samples):
        # Preparing the files
        try:
            filename = out_prefix + "_{}_gold.samples".format(k)
            with open(filename, 'w') as output_file:
                print >>output_file, "\t".join(gold_sample)

            filename = out_prefix + "_{}_source.samples".format(k)
            with open(filename, "w") as output_file:
                print >>output_file, "\t".join(source_sample)

        except IOError:
            msg = "can't write in dir {}".format(out_dir)
            raise ProgramError(msg)

    # Preparing the files
    nb = len(same_samples)
    keepSamples(
        [gold_prefix]*nb + [source_prefix]*nb,
        [out_prefix + "_{}_gold.samples".format(i) for i in range(nb)] +
        [out_prefix + "_{}_source.samples".format(i) for i in range(nb)],
        [out_prefix + "_{}_gold".format(i) for i in range(nb)] +
        [out_prefix + "_{}_source".format(i) for i in range(nb)],
        use_sge,
        transpose=True,
    )

    # Creating reports
    # The diff file
    diff_file = None
    try:
        diff_file = open(final_out_prefix + ".diff", 'w')
    except IOError:
        msg = "{}: can't write file".format(final_out_prefix + ".diff")
        raise ProgramError(msg)
    # The diff file header
    print >>diff_file, "\t".join(["name", "famID_source", "indID_source",
                                  "famID_gold", "indID_gold", "geno_source",
                                  "geno_gold"])

    # The concordance file
    concordance_file = None
    try:
        concordance_file = open(final_out_prefix + ".concordance", "w")
    except IOError:
        msg = "{}: can't write file".format(final_out_prefix + ".concordance")
        raise ProgramError(msg)
    # The concordance file header
    print >>concordance_file, "\t".join(["famID_source", "indID_source",
                                         "famID_gold", "indID_gold",
                                         "nb_geno", "nb_diff",
                                         "concordance"])

    for k in range(nb):
        # The samples
        gold_sample, source_sample = same_samples[k]

        # Reading the source freq file
        source_genotypes = {}
        filename = out_prefix + "_{}_source.tped".format(k)
        try:
            with open(filename, 'r') as input_file:
                for line in input_file:
                    row = line.rstrip("\r\n").split("\t")
                    marker_name = row[1]
                    geno = row[4]
                    if geno != "0 0":
                        geno = set(geno.split(" "))
                        source_genotypes[marker_name] = geno
        except IOError:
            msg = "{}: no such file".format(filename)
            raise ProgramError(msg)

        # Reading the gold freq file
        gold_genotypes = {}
        filename = out_prefix + "_{}_gold.tped".format(k)
        try:
            with open(filename, 'r') as input_file:
                for line in input_file:
                    row = line.rstrip("\r\n").split("\t")
                    marker_name = row[1]
                    geno = row[4]
                    if geno != "0 0":
                        geno = set(geno.split(" "))
                        gold_genotypes[marker_name] = geno
        except IOError:
            msg = "{}: no such file".format(filename)
            raise ProgramError(msg)

        # The genotyped markers in both, with their number
        genotyped_markers = source_genotypes.viewkeys() & \
            gold_genotypes.viewkeys()
        nb_genotyped = len(genotyped_markers)

        # Finding the number of differences, and creating the diff file
        nb_diff = 0
        for marker_name in genotyped_markers:
            # Getting the genotypes
            source_geno = source_genotypes[marker_name]
            gold_geno = gold_genotypes[marker_name]

            # Comparing the genotypes
            if source_geno != gold_geno:
                nb_diff += 1

                source_geno_print = list(source_geno - {"0"})
                source_geno_print.sort()
                if len(source_geno_print) == 1:
                    source_geno_print.append(source_geno_print[0])
                gold_geno_print = list(gold_geno - {"0"})
                gold_geno_print.sort()
                if len(gold_geno_print) == 1:
                    gold_geno_print.append(gold_geno_print[0])
                # We print in the diff file
                print >>diff_file, "\t".join([marker_name,
                                              "\t".join(source_sample),
                                              "\t".join(gold_sample),
                                              "/".join(source_geno_print),
                                              "/".join(gold_geno_print)])

        # Creating the concordance file
        concordance = 0.0
        if nb_genotyped != 0:
            concordance = (nb_genotyped - nb_diff) / float(nb_genotyped)
        print >>concordance_file, "\t".join(["\t".join(source_sample),
                                             "\t".join(gold_sample),
                                             str(nb_genotyped), str(nb_diff),
                                             str(concordance)])

    # Closing the output files
    diff_file.close()
    concordance_file.close()

    # Deleating the temporary directory
    try:
        shutil.rmtree(out_dir)
    except IOError:
        print >>sys.stderr, "   - Can't delete {}".format(out_dir)
def compute_statistics ( out_dir , gold_prefix , source_prefix , same_samples , use_sge , final_out_prefix ) : if not os . path . isdir ( out_dir ) : try : os . mkdir ( out_dir ) except OSError : msg = "{}: file exists" . format ( out_dir ) raise ProgramError ( msg ) out_prefix = os . path . join ( out_dir , "tmp" ) logger . info ( "  - Subsetting the files" ) for k , ( gold_sample , source_sample ) in enumerate ( same_samples ) : try : filename = out_prefix + "_{}_gold.samples" . format ( k ) with open ( filename , 'w' ) as output_file : print >> output_file , "\t" . join ( gold_sample ) filename = out_prefix + "_{}_source.samples" . format ( k ) with open ( filename , "w" ) as output_file : print >> output_file , "\t" . join ( source_sample ) except IOError : msg = "can't write in dir {}" . format ( out_dir ) raise ProgramError ( msg ) nb = len ( same_samples ) keepSamples ( [ gold_prefix ] * nb + [ source_prefix ] * nb , [ out_prefix + "_{}_gold.samples" . format ( i ) for i in range ( nb ) ] + [ out_prefix + "_{}_source.samples" . format ( i ) for i in range ( nb ) ] , [ out_prefix + "_{}_gold" . format ( i ) for i in range ( nb ) ] + [ out_prefix + "_{}_source" . format ( i ) for i in range ( nb ) ] , use_sge , transpose = True , ) diff_file = None try : diff_file = open ( final_out_prefix + ".diff" , 'w' ) except IOError : msg = "{}: can't write file" . format ( final_out_prefix + ".diff" ) raise ProgramError ( msg ) print >> diff_file , "\t" . join ( [ "name" , "famID_source" , "indID_source" , "famID_gold" , "indID_gold" , "geno_source" , "geno_gold" ] ) concordance_file = None try : concordance_file = open ( final_out_prefix + ".concordance" , "w" ) except IOError : msg = "{}: can't write file" . format ( final_out_prefix + ".concordance" ) raise ProgramError ( msg ) print >> concordance_file , "\t" . join ( [ "famID_source" , "indID_source" , "famID_gold" , "indID_gold" , "nb_geno" , "nb_diff" , "concordance" ] ) for k in range ( nb ) : gold_sample , source_sample = same_samples [ k ] source_genotypes = { } filename = out_prefix + "_{}_source.tped" . format ( k ) try : with open ( filename , 'r' ) as input_file : for line in input_file : row = line . rstrip ( "\r\n" ) . split ( "\t" ) marker_name = row [ 1 ] geno = row [ 4 ] if geno != "0 0" : geno = set ( geno . split ( " " ) ) source_genotypes [ marker_name ] = geno except IOError : msg = "{}: no such file" . format ( filename ) raise ProgramError ( msg ) gold_genotypes = { } filename = out_prefix + "_{}_gold.tped" . format ( k ) try : with open ( filename , 'r' ) as input_file : for line in input_file : row = line . rstrip ( "\r\n" ) . split ( "\t" ) marker_name = row [ 1 ] geno = row [ 4 ] if geno != "0 0" : geno = set ( geno . split ( " " ) ) gold_genotypes [ marker_name ] = geno except IOError : msg = "{}: no such file" . format ( filename ) raise ProgramError ( msg ) genotyped_markers = source_genotypes . viewkeys ( ) & gold_genotypes . viewkeys ( ) nb_genotyped = len ( genotyped_markers ) nb_diff = 0 for marker_name in genotyped_markers : source_geno = source_genotypes [ marker_name ] gold_geno = gold_genotypes [ marker_name ] if source_geno != gold_geno : nb_diff += 1 source_geno_print = list ( source_geno - { "0" } ) source_geno_print . sort ( ) if len ( source_geno_print ) == 1 : source_geno_print . append ( source_geno_print [ 0 ] ) gold_geno_print = list ( gold_geno - { "0" } ) gold_geno_print . sort ( ) if len ( gold_geno_print ) == 1 : gold_geno_print . append ( gold_geno_print [ 0 ] ) print >> diff_file , "\t" . join ( [ marker_name , "\t" . join ( source_sample ) , "\t" . join ( gold_sample ) , "/" . join ( source_geno_print ) , "/" . join ( gold_geno_print ) ] ) concordance = 0.0 if nb_genotyped != 0 : concordance = ( nb_genotyped - nb_diff ) / float ( nb_genotyped ) print >> concordance_file , "\t" . join ( [ "\t" . join ( source_sample ) , "\t" . join ( gold_sample ) , str ( nb_genotyped ) , str ( nb_diff ) , str ( concordance ) ] ) diff_file . close ( ) concordance_file . close ( ) try : shutil . rmtree ( out_dir ) except IOError : print >> sys . stderr , "   - Can't delete {}" . format ( out_dir )
Compute the statistics.