<jupyter_start><jupyter_text>Imports<jupyter_code>%load_ext autoreload
%autoreload 2
import itertools
import math
import pandas as pd
import json
import os
import glob

from tqdm import tqdm
import seaborn as sns


import screed
import sklearn

%matplotlib inline
from path_constants import QFO_EUKARYOTA_FOLDER, ORPHEUM_BENCHMARKING_FOLDER<jupyter_output><empty_output><jupyter_text>Get ids of "good" reads to use for classification  Subset to only reads from complete protein sequences -- *will not change with input species*could maybe write these ids to file<jupyter_code>busco_mammalia_uniprot_protein_starts_with_m = []

# Use only the busco mammalia proteins
with screed.open(
    os.path.join(QFO_EUKARYOTA_FOLDER, "UP000005640_9606__busco_mammlia_odbv10.fasta")
) as records:
    for record in records:
        if record["sequence"].startswith("M"):
            busco_mammalia_uniprot_protein_starts_with_m.append(
                record["name"].split()[0].split("|")[1]
            )
print(
    "busco_mammalia_uniprot_protein_starts_with_m",
    len(busco_mammalia_uniprot_protein_starts_with_m),
)


uniprot_dna_starts_with_atg = []
with screed.open(
    os.path.join(QFO_EUKARYOTA_FOLDER, "UP000005640_9606_DNA.fasta")
) as records:
    for record in records:
        if record["sequence"].startswith("ATG"):
            uniprot_dna_starts_with_atg.append(record["name"].split()[0].split("|")[1])
print("uniprot_dna_starts_with_atg", len(uniprot_dna_starts_with_atg))
busco_mammalia_startswith_m__and__dna_startswith_atg = set(busco_mammalia_uniprot_protein_starts_with_m).intersection(uniprot_dna_starts_with_atg)
len(busco_mammalia_startswith_m__and__dna_startswith_atg)<jupyter_output><empty_output><jupyter_text>Read gold standard reading frame file<jupyter_code>parquet = os.path.join(ORPHEUM_BENCHMARKING_FOLDER, 'true_reading_frames.parquet')

true_coding_frame = pd.read_parquet(parquet)
# Create just a series (single column) from this
true_coding_frame = true_coding_frame['is_coding']
true_coding_frame.head()<jupyter_output><empty_output><jupyter_text>Read coding score csvs  Add read_id_frame and is_coding for computing metrics<jupyter_code>def add_read_id_frame_and_is_coding(df):
    df["read_id_frame"] = (
        df.read_id.astype(str) + "__frame=" + df.translation_frame.astype(str)
    )
    df["is_coding"] = df["category"] == "Coding"
    return df<jupyter_output><empty_output><jupyter_text>Get human busco mammalia reads<jupyter_code>csv = os.path.join(
    QFO_EUKARYOTA_FOLDER,
    "busco_mammalia_human_uniprot_ids_in_qfo.csv",
)

human_busco_mammalia = pd.read_csv(csv)
print(human_busco_mammalia.shape)
human_busco_mammalia.head()<jupyter_output><empty_output><jupyter_text>Function to read coding scores CSVs consistently<jupyter_code>def read_translate_csvs(globber, species=None, 
                        human_busco_mammalia=human_busco_mammalia, 
                        uniprot_dna_starts_with_atg=uniprot_dna_starts_with_atg, 
                        true_coding_frame=true_coding_frame):
    dfs = []
    
    for filename in tqdm(glob.glob(globber)):
        dirname = os.path.dirname(filename)
        parent_dir = os.path.basename(dirname)
        molecule_, ksize_ = parent_dir.split('_')
        molecule = molecule_.split('-')[1]
        ksize = int(ksize_.split('-')[1])
        
        df = pd.read_csv(filename)
        df['alphabet'] = molecule
        df['ksize'] = ksize
        df['species'] = species 
        df = add_read_id_frame_and_is_coding(df)
        
        # Get only reads generated from busco mammalia data
        df['protein_id'] = df['read_id'].map(lambda x: x.split('/')[1].split(';')[0])
        df['uniprot_id'] = df['protein_id'].str.split('|').str[1]
        df = df.query('uniprot_id in @human_busco_mammalia.source__uniprot_id')
        
        # Get only reads from transcripts starting with ATG and no Ns, and only busco mammalia
        df = df.query('uniprot_id in @busco_mammalia_startswith_m__and__dna_startswith_atg')
        
        # Only use reads and frames in the gold standard data
        df = df.query('read_id_frame in @true_coding_frame.index')
        
        # Set the read id and frame as the row names
        df = df.set_index('read_id_frame')

        dfs.append(df)
    concatenated = pd.concat(dfs, ignore_index=False)
    concatenated = concatenated.sort_index()
    return concatenated<jupyter_output><empty_output><jupyter_text>Read Busco mammalia coding scores<jupyter_code>%%time
dfs = []
globber = os.path.join(
    ORPHEUM_PIPELINE_RESULTS_FOLDER, "nf-predictorthologs--busco-mammalia-*"
)
for subfolder in glob.glob(globber):
    basename = os.path.basename(subfolder)
    print(basename)
    species = basename.split("--")[-1]
    print(species)
    csvs = os.path.join(subfolder, "translate", "*", "*.csv")
    df = read_translate_csvs(
        csvs,
        species=species,
        human_busco_mammalia=human_busco_mammalia,
        uniprot_dna_starts_with_atg=uniprot_dna_starts_with_atg,
        true_coding_frame=true_coding_frame,
    )
    dfs.append(df)
coding_scores = pd.concat(dfs)
print(coding_scores.shape)
coding_scores.head()
for d in dfs:
    print(set(d['species']))
for d in dfs:
    d.species.unique()
coding_scores.species.value_counts()<jupyter_output>{'mouse'}
{'rabbit'}
{'bat'}
{'hedgehog'}
{'shrew'}
{'platypus'}
{'night_monkey'}
{'rhino'}
{'camel'}
{'baiji'}
{'chinchilla'}
{'macaque'}
{'human'}
{'tupaia'}
{'deer_mouse'}
{'goat'}
{'spalax'}
{'koala'}<jupyter_text>Write concatenated n frames per read, categorization to file<jupyter_code>coding_scores.to_parquet(
    os.path.join(ORPHEUM_PIPELINE_RESULTS_FOLDER, "coding_scores.parquet")
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11023759608665269, 'fraction_numerical': 0.006813417190775681, 'mean_word_length': 4.569066147859922, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '507696', 'license': 'MIT', 'n_tokens_mistral': 2200, 'n_tokens_neox': 2111, 'n_words': 340, 'path': 'notebooks/figure_2_04_read_orpheum_translate_mammalia_benchmarking_results.ipynb', 'repo_name': 'czbiohub/scrnaseq-for-the-99-percent', 'size': '5724'}