<jupyter_start><jupyter_text>Read Data<jupyter_code>pd.set_option('display.max_colwidth', None)
pd.set_option('display.max_rows', None)
path = Path('..\Data')
figs_output = Path('Paper\Figures')
train_genomes_path = path / 'train_genomes.fasta'
train_genomes_metedata_path = path / 'train_genomes_info.csv'

test_genomes_path = path / 'test_genomes.fasta'
test_genomes_metedata_path = path / 'test_genomes_verified.csv'

train_dataset = GenomesData(train_genomes_path, train_genomes_metedata_path)
test_dataset = GenomesData(test_genomes_path, test_genomes_metedata_path)<jupyter_output><empty_output><jupyter_text>Check split train dataset to train and validation according to insertion date<jupyter_code>print('train dataset', len(train_dataset.data))

split = feature_selection.split_by_insertion_date(train_dataset, proportion=0.2)

train_idx, validation_idx = next(split.split())

X_train_raw, X_valid_raw = train_dataset.data[train_idx], train_dataset.data[validation_idx]
y_train, y_valid = train_dataset.y[train_idx], train_dataset.y[validation_idx]

print('train', len(X_train_raw), 'validation', len(X_valid_raw))
y_valid.value_counts()
DATE_INSERTED = 'Date Inserted'

min_validation_date = min(train_dataset.metadata.iloc[validation_idx, :][DATE_INSERTED])
min_validation_date
max_train_date = max(train_dataset.metadata.iloc[train_idx, :][DATE_INSERTED])
max_train_date
assert min_validation_date > max_train_date<jupyter_output><empty_output><jupyter_text>Feature Selection - best k features according to chi2  20% validation<jupyter_code>min_val = 50
max_val = 600
inc = 50

k_range = range(min_val, max_val + 1, inc)

split = feature_selection.split_by_insertion_date(train_dataset, proportion=0.2)
grid_search = feature_selection.perform_fs_k_best(train_dataset, k_range, split=split)
print(grid_search.best_params_)

param_name = 'k_best__k'
kvals_df = feature_selection.grid_search_results_to_df(grid_search, param_name)
kvals_df
print('Validation results per k value')

fig_path = figs_output / 'feat_selection_k.png'

fs_dict = dict(zip(list(kvals_df.index), list(kvals_df['mean_test_roc_auc'])))
create_barplots(fs_dict, width=30, x_label='Number features selected in the first feature selection step (k)', y_label=AUROC, out_path=fig_path, round_by=3, bar_label_font_size=13, save=True)<jupyter_output>Validation results per k value<jupyter_text>Feature Selection - select best feature from each cluster<jupyter_code>min_val = 0
max_val = 0.85
inc = 0.06
t_range = np.arange(min_val, max_val, inc)

# created corr matrix using feature_selection.create_corr_matrix(X_train_raw, y_train) 
corr_matrix_path = path / 'X_train_corr_mat_0.8_k450.csv'
corr_matrix_train = pd.read_csv(corr_matrix_path, index_col=0)
dist_matrix_train = feature_selection.feature_corr_to_dist_matrix(corr_matrix_train)

split = feature_selection.split_by_insertion_date(train_dataset, proportion=0.2)
grid_search_cluster = feature_selection.perform_fs_clusters(train_dataset, dist_matrix_train, t_range, split=split)
print(grid_search_cluster.best_params_)

param_name = 'cluster__threshold'
thresholds_df = feature_selection.grid_search_results_to_df(grid_search_cluster, param_name)
thresholds_df
print('Validation results per threshold value')

fig_path = figs_output / 'feat_selection_t.png'
fs_dict = dict(zip(list(thresholds_df.index), list(thresholds_df['mean_test_roc_auc'])))

create_barplots(fs_dict, width=0.04, x_label='Clustering threshold (t)', y_label=AUROC, out_path=fig_path, round_by=3, bar_label_font_size=13, y_lim=[0.7, 0.905], save=True)<jupyter_output>Validation results per threshold value<jupyter_text>Validation Results<jupyter_code>best_k = 450
best_t = 0.18<jupyter_output><empty_output><jupyter_text>Correlation 0.8 of train<jupyter_code>fs_pipeline = feature_selection.get_fs_pipeline(dist_matrix_train, k=best_k, threshold=best_t)
fs_pipeline.fit(X_train_raw, y_train)

all_pgfams_train = pd.array(fs_pipeline['vectorize'].get_feature_names())
k_best_pgfams_train = all_pgfams_train[fs_pipeline['k_best'].get_support()]
cluster_pgfams_train = k_best_pgfams_train[fs_pipeline['cluster'].get_support()]
corr_matrix_train_clust = corr_matrix_train.loc[cluster_pgfams_train, cluster_pgfams_train]
corr_clust_matrix = corr_matrix_train_clust
corr_all_matrix = corr_matrix_train

NO_CORR_REDUCTION = f'Before Correlation Reduction ({len(corr_all_matrix)} features)'
CORR_REDUCTION = f'After Correlation Reduction ({len(corr_clust_matrix)} features)'

corr_clust = [i for i in corr_matrix_train_clust.values.flatten() if i != 1]
corr_all = [i for i in corr_matrix_train.values.flatten() if i != 1]

dataset_names = [CORR_REDUCTION, NO_CORR_REDUCTION]

fig_path = figs_output / 'feature_correlation_train_without_valid.png'

create_correlation_plot(corr_clust, corr_all, dataset_names, out_path=fig_path, fig_size=(9, 4), save=True)<jupyter_output><empty_output><jupyter_text>no clustering<jupyter_code>fs_pipeline = feature_selection.get_fs_pipeline(dist_matrix_train, k=best_k, threshold=0)

fs_pipeline.fit(X_train_raw, y_train)

results_valid_no_clust = predict_and_print_results(X_valid_raw, y_valid, fs_pipeline)
fs_pipeline = feature_selection.get_fs_pipeline(dist_matrix_train, k=best_k, threshold=best_t)

fs_pipeline.fit(X_train_raw, y_train)

results_valid = predict_and_print_results(X_valid_raw, y_valid, fs_pipeline)
bar_groups = [BACC, SEN, SPE, AUPR, AUROC]
bars = [CORR_REDUCTION, NO_CORR_REDUCTION]

results_for_plot = {NO_CORR_REDUCTION: results_valid_no_clust, CORR_REDUCTION: results_valid}

fig_path = figs_output / 'validation_results.png'
create_grouped_barplots(results_for_plot, bar_groups, bars, out_path=fig_path, fig_size=(9, 6), bar_width=0.4, legend_loc='upper left', x_ticks_font_size=18, y_lim=[0.4, 1.09], save=True)<jupyter_output><empty_output><jupyter_text>Final model  Perform feature selection process on the entire train dataset, using best parameters<jupyter_code>best_k = 450
best_t = 0.18
X_train_dataset_raw, X_test_raw = train_dataset.data, test_dataset.data
y_train_dataset, y_test = train_dataset.y, test_dataset.y

# created corr matrix using feature_selection.create_corr_matrix(X_train_dataset_raw, y_train, k=450) 
corr_matrix_path = path / 'X_train_corr_mat_k450.csv'
X_train_corr_matrix_k450 = pd.read_csv(corr_matrix_path, index_col=0)
X_train_dist_mat = feature_selection.feature_corr_to_dist_matrix(X_train_corr_matrix_k450)

fs_pipeline = feature_selection.get_fs_pipeline(X_train_dist_mat, k=best_k, threshold=best_t)

fs_pipeline.fit(X_train_dataset_raw, y_train_dataset)

results = predict_and_print_results(X_test_raw, y_test, fs_pipeline)
y_train_dataset.value_counts()<jupyter_output><empty_output><jupyter_text>Extract final features<jupyter_code>all_pgfams = pd.array(fs_pipeline['vectorize'].get_feature_names())
k_best_pgfams = all_pgfams[fs_pipeline['k_best'].get_support()]
cluster_pgfams = k_best_pgfams[fs_pipeline['cluster'].get_support()]<jupyter_output><empty_output><jupyter_text>Examine features<jupyter_code>len(all_pgfams)
vectorizer = CountVectorizer(lowercase=False, binary=True, vocabulary=cluster_pgfams)
X_train_fs_clust = train_dataset.vectorize_data(vectorizer)
assert all(cluster_pgfams == vectorizer.get_feature_names())
assert all(cluster_pgfams == X_train_fs_clust.columns)
vectorizer = CountVectorizer(lowercase=False, binary=True, vocabulary=k_best_pgfams)
X_train_fs_k_best = train_dataset.vectorize_data(vectorizer)
chi2_vals, __ = sklearn.feature_selection.chi2(X_train_fs_k_best, train_dataset.y)
min(chi2_vals)
X_train_fs_k_best_hp = X_train_fs_k_best.loc[:, set(pgfam_to_count_hp.keys()).intersection(X_train_fs_k_best.columns)]
X_train_fs_k_best_nhp = X_train_fs_k_best.loc[:, set(pgfam_to_count_nhp.keys()).intersection(X_train_fs_k_best.columns)]

chi2_vals_fs_k_best_hp, __ = sklearn.feature_selection.chi2(X_train_fs_k_best_hp, train_dataset.y)
chi2_vals_fs_k_best_nhp, __ = sklearn.feature_selection.chi2(X_train_fs_k_best_nhp, train_dataset.y)
min(chi2_vals_fs_k_best_hp)
min(chi2_vals_fs_k_best_nhp)
def get_hp_to_nhp_count(X, y, pgfam):
    return feature_importance.count_hp_vs_nhp(X[X[pgfam]==1].index, y)


pgfam_to_count_hp = {}
pgfam_to_count_nhp = {}

for pgfam in k_best_pgfams:
    hp, nhp = get_hp_to_nhp_count(X_train_fs_k_best, train_dataset.y, pgfam)
    if hp > nhp:
        pgfam_to_count_hp[pgfam] = hp
    else:
        pgfam_to_count_nhp[pgfam] = nhp


min(pgfam_to_count_hp.values()), min(pgfam_to_count_nhp.values())
vectorizer = CountVectorizer(lowercase=False, binary=True, vocabulary=cluster_pgfams)
X_train_final_features = train_dataset.vectorize_data(vectorizer)

X_train_final_features_hp = X_train_final_features.loc[:, set(pgfam_to_count_hp.keys()).intersection(X_train_final_features.columns)]
X_train_final_features_nhp = X_train_final_features.loc[:, set(pgfam_to_count_nhp.keys()).intersection(X_train_final_features.columns)]

chi2_vals_final_hp, __ = sklearn.feature_selection.chi2(X_train_final_features_hp, train_dataset.y)
chi2_vals_final_nhp, __ = sklearn.feature_selection.chi2(X_train_final_features_nhp, train_dataset.y)
chi2_final_feats_hp = pd.Series(chi2_vals_final_hp, index=X_train_final_features_hp.columns)
chi2_final_feats_hp = chi2_final_feats_hp.sort_values(ascending=False)
chi2_final_feats_hp.head()
chi2_final_feats_nhp = pd.Series(chi2_vals_final_nhp, index=X_train_final_features_nhp.columns)
chi2_final_feats_nhp = chi2_final_feats_nhp.sort_values(ascending=False)
chi2_final_feats_nhp.head()
pgfams = {}

with open(path / 'PATRIC_pgfams_wspc.txt') as f:
    header = f.readline()

    for line in f:
        pgfam, product = line.strip().split('\t')
        pgfams[pgfam] = product
for pgfam in chi2_final_feats_hp.index[:15]:
    print(pgfams[pgfam])
for pgfam in chi2_final_feats_nhp.index[:15]:
    print(pgfams[pgfam])<jupyter_output>Cytidylate kinase (EC 2.7.4.25)
Reverse rubrerythrin
Phosphoribosylaminoimidazolecarboxamide formyltransferase (EC 2.1.2.3)
Electron transport complex protein RnfB
no significant homology.
Flavoprotein
[FeFe] hydrogenase (EC 1.12.7.2)
Rubrerythrin
Activator of (R)-2-hydroxyglutaryl-CoA dehydratase
"Dihydroorotate dehydrogenase (NAD(+)), catalytic subunit (EC 1.3.1.14)"
Iron-sulfur flavoprotein
"L,L-diaminopimelate aminotransferase (EC 2.6.1.83)"
"RNA methyltransferase, TrmA family"
Flavodoxin
Flavoredoxin<jupyter_text>Observe clusters to make sure they make sense<jupyter_code>clusters = fs_pipeline['cluster'].clusters_
clusters.sort(key=lambda cluster: len(cluster), reverse=True)
for i, cluster in enumerate(clusters):
    print('cluster', i)
    for pgfam_ind in cluster:
        pgfam = k_best_pgfams[pgfam_ind]
        print(pgfam, pgfams[pgfam])
    print('-----------------')<jupyter_output>cluster 0
PGF_00000939 Uncharacterized protein YhiN
PGF_00008469 Glutathione synthetase (EC 6.3.2.3)
PGF_00019484 Membrane alanine aminopeptidase N (EC 3.4.11.2)
PGF_00031127 Periplasmic thiol:disulfide interchange protein DsbA
PGF_00054476 Succinate dehydrogenase hydrophobic membrane anchor protein
PGF_00057016 "Tol biopolymer transport system, TolR protein"
PGF_00057028 Tol-Pal system peptidoglycan-associated lipoprotein PAL
PGF_00064460 Ubiquinone biosynthesis regulatory protein kinase UbiB
PGF_00064685 Uncharacterized ferredoxin-like protein YfhL
PGF_00415581 23S rRNA (uridine(2552)-2'-O)-methyltransferase (EC 2.1.1.166)
PGF_00417755 3-demethylubiquinol 3-O-methyltransferase (EC 2.1.1.64) @ 2-polyprenyl-6-hydroxyphenyl methylase (EC 2.1.1.222)
PGF_00423927 4-hydroxybenzoate polyprenyltransferase (EC 2.5.1.39)
PGF_00648054 Tol-Pal system protein TolQ
PGF_00649271 Cell division protein BolA
PGF_00906364 Hydroxyacylglutathione hydrolase (EC 3.1.2.6)
PGF_01785347 Adenosine (5')-pentaph[...]<jupyter_text>Correlation entire train<jupyter_code>X_train_250_feats_corr_mat = X_train_corr_matrix_k450.loc[cluster_pgfams, cluster_pgfams]
X_train_250_feats_corr_mat.head()
corr_clust_matrix = X_train_250_feats_corr_mat
corr_all_matrix = X_train_corr_matrix_k450

corr_clust = [i for i in corr_clust_matrix.values.flatten() if i != 1]
corr_all = [i for i in corr_all_matrix.values.flatten() if i != 1]

dataset_names = dataset_names = [f'{len(corr_clust_matrix)} features', f'{len(corr_all_matrix)} features']

fig_path = figs_output / 'feature_correlation_entire_train.pdf'

create_correlation_plot(corr_clust, corr_all, dataset_names, out_path=fig_path)<jupyter_output><empty_output><jupyter_text>Final WSPC modelA random forest model, which includes the selected genes as features<jupyter_code>WSPC, results = train_and_predict(X_train_dataset_raw, y_train_dataset, X_test_raw, y_test, features=cluster_pgfams)
predictions, predictions_probs = predict(X_test_raw, WSPC)
predict_df = pd.DataFrame({'predictions':predictions, 'predictions_probs':predictions_probs}, index=X_test_raw.index)
predict_df.head()
save = False

if save:
    #save predictions

    predict_df.to_csv(path/'WSPC_test_preds'/'wspc_test_preds.csv')
    
    #save model
    out_path = r'WSPC Model'
    model_filename = 'WSPC_model'
    util.save_model(out_path, model_filename, WSPC)<jupyter_output><empty_output><jupyter_text>Novel species All classified species in train<jupyter_code>train_species = taxa.get_genomes_species(train_dataset.data.index)
len(train_species)
train_genera = taxa.get_genomes_tax(train_dataset.data.index, taxa.GENUS)
len(train_genera)<jupyter_output><empty_output><jupyter_text>Genera in test<jupyter_code>test_genera = taxa.get_genomes_tax(test_dataset.data.index, taxa.GENUS)
len(test_genera)<jupyter_output><empty_output><jupyter_text>All classified species in test<jupyter_code>test_species = taxa.get_genomes_species(test_dataset.data.index)
len(test_species)
test_species_not_in_train = taxa.get_novel_species(train_dataset.data.index, test_dataset.data.index)
len(test_species_not_in_train)
test_novel_genome_ids = taxa.get_genomes_with_tax(test_dataset.data.index, test_species_not_in_train, taxa.SPECIES)
novel_species_results = predict_and_print_results(X_test_raw[test_novel_genome_ids], y_test[test_novel_genome_ids], WSPC)<jupyter_output>false_positive: 7,total NHPs: 31
false_negative: 6,total HPs: 32
BACC: 0.79
Sensitivity: 0.81
Specificity: 0.77
AUPR: 0.89
AUROC: 0.87<jupyter_text>Plot Results<jupyter_code>ENTIRE_TEST = 'Entire Test'
NOVEL_SPECIES = 'Novel Species'
results_for_plot = {ENTIRE_TEST: results, NOVEL_SPECIES: novel_species_results}
bar_groups = [BACC, SEN, SPE, AUPR, AUROC]
bars = [ENTIRE_TEST, NOVEL_SPECIES]

fig_path = figs_output / 'wspc_test_results.pdf'
create_grouped_barplots(results_for_plot, bar_groups, bars, out_path=fig_path, fig_size=(8, 5), bar_width=0.4, legend_loc='upper left', y_lim=(0.1, 1.199), x_ticks_font_size=16)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09810040705563093, 'fraction_numerical': 0.025237449118046134, 'mean_word_length': 9.027891156462585, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 12, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '529273', 'license': 'MIT', 'n_tokens_mistral': 5961, 'n_tokens_neox': 5615, 'n_words': 845, 'path': 'Notebooks/WSPC_train_pipeline.ipynb', 'repo_name': 'shakedna1/wspc_rep', 'size': '14740'}