Python Benchmark.add_method方法代码示例

# 需要导入模块: from benchmark import Benchmark [as 别名]
# 或者: from benchmark.Benchmark import add_method [as 别名]

#.........这里部分代码省略.........
    proteins_in_testset = dataset_properties.query('dataset_id in @benchmark_dataset_id')['protein'].values

    print("Start processing alignment files...")

    ########## Iterate over proteins
    counter=0
    it = -1
    while counter < n_proteins:
        it += 1
        proteins_subset = proteins_in_testset[(it * n_proteins):( (it+1)*n_proteins)]
        np.random.shuffle(proteins_subset)
        for protein in proteins_subset:

            if counter >=  n_proteins:
                break

            protein = protein.strip()

            alignment_file  = alignment_dir + "/" + protein + ".filt.psc"
            psipred_file    = psipred_dir + "/" + protein + ".filt.withss.a3m.ss2"
            netsurfp_file   = netsurfp_dir + "/" + protein + ".filt.netsurfp"
            mi_file = mi_dir + "/" + protein + ".filt.mi.pc.mat"
            omes_file = omes_dir + "/" + protein + "filt.omes.fodoraldrich.mat"
            eval_file = evaluation_dir + "/" + protein.strip() + "." + method_name
            pll_braw_file = None
            cd_braw_file = None
            pcd_braw_file = None
            bayposterior_mat_file = None
            bayfactor_mat_file = None


            if os.path.exists(eval_file) and not update:
                print("Evaluation file {0} already exists. Skip this protein. ".format(eval_file))
                continue

            if not os.path.exists(alignment_file):
                print("Alignment file {0} does not exist. Skip this protein. ".format(alignment_file))
                continue

            if not os.path.exists(psipred_file):
                print("    Psipred file {0} does not exist. Skip protein {1}!".format(psipred_file, protein))
                continue

            if not os.path.exists(netsurfp_file):
                print("    NetsurfP file {0} does not exist. Skip protein {1}!".format(netsurfp_file, protein))
                continue


            if pll_braw_dir is not None:
                pll_braw_file = pll_braw_dir + "/" + protein.strip() + ".filt.braw.gz"
                if not os.path.exists(pll_braw_file):
                    print("    pLL braw file {0} does not exist. Skip protein {1}!".format(pll_braw_file, protein))
                    continue

            if cd_braw_dir is not None:
                cd_braw_file = cd_braw_dir + "/" + protein.strip() + ".filt.braw.gz"
                if not os.path.exists(cd_braw_file):
                    print("    CD braw file {0} does not exist. Skip protein {1}!".format(cd_braw_file, protein))
                    continue

            if pcd_braw_dir is not None:
                pcd_braw_file = pcd_braw_dir + "/" + protein.strip() + ".filt.braw.gz"
                if not os.path.exists(pcd_braw_file):
                    print("    PCD braw file {0} does not exist. Skip protein {1}!".format(pcd_braw_file, protein))
                    continue

            if bayposterior_mat_dir is not None:
                bayposterior_mat_file = bayposterior_mat_dir + "/" + protein.strip() + ".bayesian_3comp_pLL.mat"
                if not os.path.exists(bayposterior_mat_file):
                    print("    bayesian posterior mat file {0} does not exist. Skip protein {1}!".format(bayposterior_mat_file, protein))
                    continue

            if bayesfactor_mat_dir is not None:
                bayfactor_mat_file = bayesfactor_mat_dir + "/" + protein.strip() + ".bayesian_3comp_pLL.mat"
                if not os.path.exists(bayfactor_mat_file):
                    print("    bayes factor mat file {0} does not exist. Skip protein {1}!".format(bayfactor_mat_file, protein))
                    continue


            BCP = BayesianContactPredictor(alignment_file)
            BCP.set_contact_prior_model(rf_clf, rf_meta)

            BCP.set_n_threads(n_threads)
            BCP.set_sequence_separation(sequence_separation)
            BCP.set_contact_threshold(contact_threshold)

            BCP.contact_prior(
                psipred_file, netsurfp_file, mi_file, omes_file,
                pll_braw_file, cd_braw_file, pcd_braw_file, bayposterior_mat_file, bayfactor_mat_file
            )


            contact_prior_mat = BCP.get_contact_prior(contact=1)
            meta = {
                'opt_code' : 1,
                'rf' : rf_meta
            }

            b.add_method(protein, method_name, contact_prior_mat, meta, apc=False, update=update)
            counter += 1

# 需要导入模块: from benchmark import Benchmark [as 别名]
# 或者: from benchmark.Benchmark import add_method [as 别名]

#.........这里部分代码省略.........


    #debugging
    # evaluation_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/evaluation/"
    # property_files_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/dataset/dataset_properties/"
    # alignment_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/psicov/"
    # psipred_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/psipred/hhfilter_results_n5e01/"
    # netsurfp_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/netsurfp/"
    # mi_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/contact_prediction/local_methods/mi_pc/"
    # omes_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/contact_prediction/local_methods/omes_fodoraldrich/"
    #
    # method_name = "pLL_3comp_reg100prec01mu_100k_ncthr8"
    # braw_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/contact_prediction/ccmpred-pll-centerv/braw/"
    # qij_dir = "/home/vorberg/work/data/benchmarkset_cathV4.1/contact_prediction/ccmpred-pll-centerv/qij/"
    # contact_prior_model_file = "/home/vorberg/work/data/bayesian_framework/contact_prior/random_forest/new_pipeline_5folds/random_forest/classweightNone_noncontactthr8/100000contacts_500000noncontacts_5window_8noncontactthreshold_maxfeatures030/random_forest_nestimators1000_maxfeatures0.3_maxdepth100_minsamplesleaf10_75features.pkl"
    # coupling_prior_parameters_file = "/home/vorberg/work/data//bayesian_framework/mle_for_couplingPrior_cath4.1/ccmpred-pll-centerv/3/reg_prec100_mu01/diagonal_100000_nrcomponents3_noncontactthr8/parameters"
    # sequence_separation = 8
    # contact_threshold = 8
    # n_threads = 8
    # n_proteins = 50


    ###########  Setup dataset_id
    dataset_properties = pd.DataFrame()
    for id, property_file in enumerate(sorted(glob.glob(property_files_dir + "/*"))):
        properties = pd.read_table(property_file)
        properties['id'] = id + 1
        properties.columns = ['protein', 'resol', 'CATH-topology', 'domlength', 'alilength', 'dataset_id']
        dataset_properties = dataset_properties.append(properties, ignore_index=True)


    ########## Setup Benchmark framework
    b = Benchmark(evaluation_dir)

    ######### Load contact prior model
    rf_clf, rf_meta = BayesianContactPredictor.load_contact_prior_model(contact_prior_model_file)

    ######### Load coupling prior parameters
    coupling_prior_parameters = BayesianContactPredictor.load_coupling_prior_hyperparameters(coupling_prior_parameters_file)

    #get all existing braw files
    benchmark_dataset_id = [6, 7, 8]
    proteins_in_testset = dataset_properties.query('dataset_id in @benchmark_dataset_id')['protein'].values
    braw_files = [braw_dir+"/" + protein.strip() + ".filt.braw.gz" for protein in proteins_in_testset]
    braw_files_existing = [braw_file for braw_file in braw_files if os.path.exists(braw_file)]
    braw_files_shuff = random.sample(braw_files_existing[:n_proteins], len(braw_files_existing[:n_proteins]))

    print("Start processing {0} braw files...".format(len(braw_files_shuff)))

    ########## Iterate over proteins
    for braw_file in braw_files_shuff:
        #braw_file ='/home/vorberg/work/data/benchmarkset_cathV4.1/contact_prediction/ccmpred-pll-centerv/braw//2j5yA00.filt.braw.gz'
        protein = braw_file.split("/")[-1].split(".")[0]

        alignment_file = alignment_dir + "/" + protein + ".filt.psc"
        psipred_file = psipred_dir + "/" + protein + ".filt.withss.a3m.ss2"
        netsurfp_file = netsurfp_dir + "/" + protein + ".filt.netsurfp"
        mi_file = mi_dir + "/" + protein + ".filt.mi.pc.mat"
        omes_file = omes_dir + "/" + protein + ".filt.omes.fodoraldrich.mat"
        qij_file = qij_dir + "/" + protein + ".filt.bqij.gz"

        if not os.path.exists(alignment_file):
            print("Alignment file {0} does not exist. Skip this protein. ".format(alignment_file))
            continue

        if not os.path.exists(braw_file):
            print("binary Raw file {0} does not exist. Skip this protein. ".format(braw_file))
            continue

        if not os.path.exists(qij_file):
            print("Qij file {0} does not exist. Skip this protein. ".format(qij_file))
            continue


        print("Compute posterior probabilities for contact with Bayesian model for protein {0}".format(protein))
        BCP = BayesianContactPredictor(alignment_file)
        BCP.set_contact_prior_model(rf_clf, rf_meta)
        BCP.set_coupling_prior_parameters(coupling_prior_parameters)

        BCP.set_n_threads(n_threads)
        BCP.set_sequence_separation(sequence_separation)
        BCP.set_contact_threshold(contact_threshold)

        BCP.contact_prior(psipred_file, netsurfp_file, mi_file, omes_file)
        BCP.contact_likelihood(braw_file, qij_file)
        BCP.contact_posterior()

        contact_prior_mat       = BCP.get_contact_prior(contact=1)
        contact_posterior_mat   = BCP.get_contact_posterior(contact=1)
        posterior_meta = BCP.get_meta()

        b.add_method(protein.strip(), method_name, contact_posterior_mat, posterior_meta, apc=False, update=True)
        b.add_method(protein.strip(), "rf_contact_prior", contact_prior_mat, posterior_meta, apc=False, update=True)

        if evaluate_bayes_factor:
            contact_likelihood_mat = BCP.get_contact_likelihood(contact=1, normalized=False, bayes_factor=True)
            b.add_method(protein.strip(), method_name+"_logbf", contact_likelihood_mat, posterior_meta, apc=False, update=True)
        if evaluate_likelihood:
            contact_likelihood_mat = BCP.get_contact_likelihood(contact=1, normalized=True, bayes_factor=False)
            b.add_method(protein.strip(), method_name+"_llik", contact_likelihood_mat, posterior_meta, apc=False, update=True)