<jupyter_start><jupyter_text>Imports<jupyter_code>%load_ext autoreload
%autoreload 2

import glob
import os
import re

from joblib import Parallel, delayed
import kmermaid_utils
import numpy as np
import pandas as pd
import scanpy as sc
import seaborn as sns
import sig_utils
import sourmash
from tqdm import tqdm

pd.options.display.max_colwidth = 500
pd.options.display.max_columns = 50

from sourmash_search_utils import read_and_process_sourmash_search_results<jupyter_output><empty_output><jupyter_text>Constants (ksizes, scaled)<jupyter_code>ls -lha /home/olga/data_lg/data_sm_copy/immune-evolution/kmer-signatures
sig_outdir_base = '/home/olga/data_lg/data_sm_copy/immune-evolution/kmer-signatures'

bat_test_dir = os.path.join(sig_outdir_base, "3--test-bat")

aggregated_results_dir = os.path.join(bat_test_dir, '4--aggregated-results')
! mkdir -p $aggregated_results_dir
SHARED_CELLTYPES = set(['Alveolar Epithelial Type 2',
 'B cell',
 'Capillary',
 'Dendritic',
 'Fibroblast',
 'Macrophage',
 'Monocyte',
 'Natural Killer T cell',
 'Smooth Muscle and Myofibroblast',
 'T cell'])
SHARED_CELLTYPES<jupyter_output><empty_output><jupyter_text>Def describe<jupyter_code>def describe(df, random=False):
    print(df.shape)
    print("--- First 5 entries ---")
    display(df.head())
    if random:
        print('--- Random subset ---')
        display(df.sample(5))<jupyter_output><empty_output><jupyter_text>Remove rogue TQDM instances<jupyter_code>def remove_rogue_tqdm():
    import tqdm
    try:
        tqdm._instances.clear()
    except AttributeError:
        pass<jupyter_output><empty_output><jupyter_text>Read one2one 5had<jupyter_code>one2one = sc.read(
#     "/home/phoenix/data_lg/data_sm_copy/code/immune-evolution/h5ads/human-lemur-mouse-bat__lung_only_copy.h5ad"
    "/home/olga/data_lg/data_sm_copy/immune-evolution/h5ads/human-lemur-mouse-bat/human-lemur-mouse-bat__lung_only.h5ad"
)
one2one.obs = one2one.obs.reset_index().set_index('cell_id')

one2one.obs.head()<jupyter_output><empty_output><jupyter_text>Analyze Sourmash search  Functions / tools to add cell type ontologies  Search directory output  Column renamer<jupyter_code>cols = ["cell_ontology_class", "narrow_group", "broad_group", "compartment_group"]

mouse_renamer = {col: f"mouse_{col}"  for col in cols}
bat_renamer = {col: f"bat_{col}" for col in cols}<jupyter_output><empty_output><jupyter_text>Read output CSVs -- Parallelized<jupyter_code># from sourmash_search_utils import read_and_process_sourmash_search_results


# bat_kmermaid_base = (
#     "/home/olga/data_lg/data_sm_copy/immune-evolution/pipeline-results/bat/kmermaid/ksize_8"
# )
# bat_search_dir = os.path.join(
#     bat_kmermaid_base,
#     "search_bat_in_mouse_no_ribosome_no_containment_no_dissociation__per_celltype_remove_common_hashes__remove_refseq_ribosomal",
# #     "search_bat_in_mouse_no_ribosome_no_containment_no_dissociation__per_celltype_remove_common_hashes",
# )

# total = 612648

# globber = os.path.join(bat_search_dir, "alphabet-dayhoff__ksize-39__scaled-10", "*.csv")

# # dfs = Parallel(n_jobs=96)(
# #     delayed(read_and_process_sourmash_search_results)(
# #         csv, adata=one2one, column_renamer=bat_renamer, query_species='bat'
# #     )
# #     for csv in tqdm(glob.iglob(globber), total=total)
# # )

# # search_results = pd.concat(dfs, ignore_index=True)
# # print(search_results.shape)
# # search_results.head()
from sourmash_search_utils import add_ontologies<jupyter_output><empty_output><jupyter_text>Read output CSVs<jupyter_code>bat_test_dir


globber = os.path.join(bat_test_dir, "3-*", "*", "*.csv")

total = sum(1 for _ in glob.iglob(globber))
total
dna_count = sum(1 for x in glob.iglob(globber) if "DNA" in x)
dna_count

dfs = []

celltype_col = "broad_group"
query_species = "bat"
query_cell_id_col = f"{query_species}_cell_id"
query_celltype_col = f"{query_species}_{celltype_col}"

sbt_species = "mouse"
sbt_celltype_col = f"{sbt_species}_{celltype_col}"


for csv in tqdm(glob.iglob(globber), total=total):
    df = read_and_process_sourmash_search_results(
        csv,
        adata=one2one,
        column_renamer=bat_renamer,
        query_species="bat",
        shared_celltypes=SHARED_CELLTYPES,
    )

    df["database_type"] = os.path.dirname(csv).split("/")[-2].split("--")[-1]

    dfs.append(df)
# break
search_results = pd.concat(dfs, ignore_index=True)
print(search_results.shape)
search_results.head()
globber = os.path.join(bat_test_dir, "3-*", "*", "*.csv")

total = sum(1 for _ in glob.iglob(globber))
search_results.sketch_id.value_counts()
csv
dna_csv = '/home/olga/data_lg/data_sm_copy/immune-evolution/kmer-signatures/3--test-bat/3-3--mouse-search-results--min-10-percent-shared-kmers-remove-common-kmers/alphabet-dayhoff__ksize-51__scaled-10/bat2-LU__AAACCTGAGGTACTCT.csv'

dna_csv
csv
total
search_results.search_method.value_counts()
search_results.sketch_id.value_counts()
describe(search_results)
n_results_per_sketch_id = search_results.groupby(['sketch_id', 'alphabet', 'ksize']).size()
n_results_per_sketch_id.name = 'n_results'
n_results_per_sketch_id = n_results_per_sketch_id.reset_index()
describe(n_results_per_sketch_id)
n_results_per_sketch_id.nlargest(5, 'n_results')
n_results_per_sketch_id.loc[n_results_per_sketch_id.n_results.argmax()]<jupyter_output><empty_output><jupyter_text>Constants (`classification_groupby`, `celltype_col`, `query_organism`, `scoring_groupby`)<jupyter_code>celltype_col = "broad_group"
sbt_organism = "mouse"
query_organism = "bat"

classification_groupby = [
    f"{query_organism}_cell_id",
    "database_type",
    "search_method",
    "alphabet",
    "ksize",
]

predicted_celltype_col = f"{sbt_organism}_{celltype_col}"
ground_truth_celltype_col = f"{query_organism}_{celltype_col}"

scoring_groupby = ["database_type", "search_method", "alphabet", "ksize"]<jupyter_output><empty_output><jupyter_text>Add unifying column of predicted and ground truth celltypes<jupyter_code>search_results['groundtruth_celltype'] = search_results[ground_truth_celltype_col]
search_results['predicted_celltype'] = search_results[predicted_celltype_col]<jupyter_output><empty_output><jupyter_text>write results to parquet file<jupyter_code>%%time

parquet_fn = "sourmash-search-results.parquet"
parquet = os.path.join(aggregated_results_dir, parquet_fn)
search_results.to_parquet(parquet)

parquet<jupyter_output>CPU times: user 3.48 s, sys: 540 ms, total: 4.02 s
Wall time: 4.09 s<jupyter_text>Predict cell type with top hit  Get top N<jupyter_code>%%time

predicted_cells_top_hit = search_results.groupby(
    classification_groupby, as_index=False, group_keys=False
).apply(lambda x: x.nlargest(1, 'similarity'))

predicted_cells_top_hit = predicted_cells_top_hit.set_index(f"{query_organism}_cell_id")

print(predicted_cells_top_hit.shape)
predicted_cells_top_hit.head()<jupyter_output>(216998, 17)
CPU times: user 16min 48s, sys: 39.9 s, total: 17min 28s
Wall time: 17min 20s<jupyter_text>Add unifying column of predicted and ground truth celltypes<jupyter_code>predicted_cells_top_hit['groundtruth_celltype'] = predicted_cells_top_hit[ground_truth_celltype_col]
predicted_cells_top_hit['predicted_celltype'] = predicted_cells_top_hit[predicted_celltype_col]<jupyter_output><empty_output><jupyter_text>Write top N to file<jupyter_code>%%time

parquet_fn = "sourmash-search-results--top-hit.parquet"
parquet = os.path.join(aggregated_results_dir, parquet_fn)
predicted_cells_top_hit.to_parquet(parquet)
parquet<jupyter_output><empty_output><jupyter_text>Get ideal ksizes per alphabet sizes<jupyter_code>alphabet_sizes = pd.Series(dict(dayhoff=6, protein=20))  # hp=2,
np.log10(alphabet_sizes.pow(21 / 3))
ideal_ksizes = np.ceil(21 * np.log(4) / np.log(alphabet_sizes))
ideal_ksizes
np.log10(alphabet_sizes.pow(ideal_ksizes))
k_nuc_sizes = np.arange(21, 91, 3)
k_nuc_sizes
k_aa_sizes = k_nuc_sizes/3
k_aa_sizes
dna_entropy = pd.Series(k_nuc_sizes * np.log10(4), name='entropy', index=k_nuc_sizes)
dna_entropy = dna_entropy.reset_index()
dna_entropy = dna_entropy.rename(columns={'index': 'k_nuc'})
dna_entropy['alphabet'] = 'dna'
describe(dna_entropy)
aa_nuc_alphabet_sizes = alphabet_sizes.copy()
aa_nuc_alphabet_sizes['dna'] = 4
alphabet_complexity = alphabet_sizes.to_frame().apply(
    lambda x: pd.Series(k_aa_sizes * np.log10(x.iloc[0]), name=x.name), axis=1
)
alphabet_complexity.columns = k_nuc_sizes
alphabet_complexity = alphabet_complexity.unstack().reset_index()
alphabet_complexity = alphabet_complexity.rename(
    columns={"level_0": "k_nuc", "level_1": "alphabet", 0: "entropy"}
)
alphabet_complexity = pd.concat([alphabet_complexity, dna_entropy])
alphabet_complexity["k_aa"] = (alphabet_complexity["k_nuc"] / 3).astype(int)
alphabet_complexity["sigma"] = alphabet_complexity.alphabet.map(aa_nuc_alphabet_sizes)
alphabet_complexity["alphabet_type"] = alphabet_complexity.alphabet.map(
    lambda x: "nuc" if x == "dna" else "aa"
)
alphabet_complexity["alphabet_label"] = alphabet_complexity.apply(
    lambda x: "{alphabet} $|\Sigma_{{{alphabet_type}}}| = {sigma}$".format(**x), axis=1
)
describe(alphabet_complexity)
sorted(alphabet_complexity.alphabet_label.unique())
alphabet_order = [
    "dna $|\\Sigma_{nuc}| = 4$",
    "protein $|\\Sigma_{aa}| = 20$",
    "dayhoff $|\\Sigma_{aa}| = 6$",
#     "hp $|\\Sigma_{aa}| = 2$",
]

g = sns.catplot(
    x="k_aa",
    data=alphabet_complexity,
    y="entropy",
    hue="alphabet_label",
    kind="point",
    size=3,
    aspect=1.75,
    hue_order=alphabet_order,
)

entropy_min = alphabet_complexity.query(
    'k_nuc == 63 and alphabet == "hp"'
).entropy.iloc[0]
entropy_max = alphabet_complexity.query(
    'k_nuc == 60 and alphabet == "dayhoff"'
).entropy.iloc[0]

for ax in g.axes.flat:
    xmin, xmax = ax.get_xlim()
    ax.fill_between(
        (xmin, xmax), entropy_min, entropy_max, color="lightgrey", zorder=-1
    )
    ax.axhline(entropy_min, color="grey", linestyle="--", zorder=-1)
    ax.axhline(entropy_max, color="grey", linestyle="--", zorder=-1)

g.set(
    ylabel=r"$k \times \log_{10} \left( | \Sigma | \right)$",
    xlim=(xmin, xmax),
    ylim=(0, 30),
    xlabel=r"$k_{aa}$ size",
)
figure_folder = (
    "/home/olga/googledrive/ImmuneEvolution/analysis/alphabet_size_vs_complexity"
)
! mkdir -p $figure_folder
pdf = os.path.join(figure_folder, "alphabet_size_vs_complexity.pdf")
g.savefig(pdf)
alphabet_complexity.query('k_nuc == 39')
alphabet_complexity.query('k_nuc == 51')
alphabet_complexity.query()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.131375726952045, 'fraction_numerical': 0.014300695967203738, 'mean_word_length': 6.503576537911302, 'pattern_counts': {'":': 2, '<': 59, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '419608', 'license': 'MIT', 'n_tokens_mistral': 4274, 'n_tokens_neox': 4072, 'n_words': 578, 'path': 'notebooks/326_read_bat2mouse_reog_search_results.ipynb', 'repo_name': 'czbiohub/scrnaseq-for-the-99-percent', 'size': '10489'}