maomlab/CryptoCEN-ExpressionScatter

1import datasets2import streamlit as st3import numpy as np4import pandas as pd5import altair as alt6 7st.set_page_config(layout='wide')8 9# parse out gene_ids from URL query args to it's possible to link to this page10query_params = st.query_params11if "gene_id_1" in query_params.keys():12    gene_id_1 = query_params["gene_id_1"]13else:14    gene_id_1 = "CNAG_04365"15 16if "gene_id_2" in query_params.keys():17    gene_id_2 = query_params["gene_id_2"]18else:19    gene_id_2 = "CNAG_04222"20 21    22 23 24st.markdown("""25# CryptoCEN Expression Scatter26**CryptoCEN** is a co-expression network for *Cryptococcus neoformans* built on 1,524 RNA-seq runs across 34 studies.27A pair of genes are said to be co-expressed when their expression is correlated across different conditions and28is often a marker for genes to be involved in similar processes. 29 30To Cite:31MJ O'Meara, JR Rapala, CB Nichols, C Alexandre, B Billmyre, JL Steenwyk, A Alspaugh,32TR O'Meara CryptoCEN: A Co-Expression Network for Cryptococcus neoformans reveals33novel proteins involved in DNA damage repair34* Code available at https://github.com/maomlab/CalCEN/tree/master/vignettes/CryptoCEN35* Full network and dataset: https://huggingface.co/datasets/maomlab/CryptoCEN36 37## Plot scatter plot expression for a pair of genes across studies.38Put in the ``CNAG_#####`` gene_id for two genes.39""")40 41h99_transcript_annotations = datasets.load_dataset(42    path = "maomlab/CryptoCEN",43    data_files = {"h99_transcript_annotations": "h99_transcript_annotations.tsv"})44h99_transcript_annotations = h99_transcript_annotations["h99_transcript_annotations"].to_pandas()45 46estimated_expression_meta = datasets.load_dataset(47    path = "maomlab/CryptoCEN",48    data_files = {"estimated_expression_meta": "Data/estimated_expression_meta.tsv"})49estimated_expression_meta = estimated_expression_meta["estimated_expression_meta"].to_pandas()50 51estimated_expression = datasets.load_dataset(52    path = "maomlab/CryptoCEN",53    data_files = {"estimated_expression": "estimated_expression_matrix.parquet"})54estimated_expression = estimated_expression["estimated_expression"].to_pandas()55 56#DEBUG57print(f"estimated_expression shape: {estimated_expression.shape}")58 59col1, col2, col3, padding = st.columns(spec = [0.2, 0.2, 0.2, 0.4])60with col1:61    gene_id_1 = st.text_input(62        label = "Gene ID 1",63        value = f"{gene_id_1}",64        max_chars = 10,65        help = "CNAG Gene ID e.g. CNAG_04365")66 67with col2:68    gene_id_2 = st.text_input(69        label = "Gene ID 2",70        value = f"{gene_id_2}",71        max_chars = 10,72        help = "CNAG Gene ID e.g. CNAG_04222")73 74# check the user input75try:76    cnag_id_1 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_1]["cnag_id"].values[0]77    gene_symbol_1 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_1]["gene_symbol"].values[0]78    description_1 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_1]["description"].values[0]79except:80    st.error(f"Unable to locate cnag_id for Gene ID 1: {gene_id_1}, it should be of the form 'CNAG_######'")81 82try:83    cnag_id_2 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_2]["cnag_id"].values[0]84    gene_symbol_2 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_2]["gene_symbol"].values[0]85    description_2 = h99_transcript_annotations.loc[h99_transcript_annotations["gene_id"] == gene_id_2]["description"].values[0]86except:87    st.error(f"Unable to locate cnag_id for Gene ID 2: {gene_id_2}, it should be of the form 'CNAG_######'")88 89chart_data = pd.DataFrame({90    "gene_id_1": gene_id_1,91    "gene_id_2": gene_id_2,92    "expression_1": estimated_expression.loc[h99_transcript_annotations["gene_id"] == gene_id_1].to_numpy()[0],93    "expression_2": estimated_expression.loc[h99_transcript_annotations["gene_id"] == gene_id_2].to_numpy()[0],94    "log_expression_1": np.log10(estimated_expression.loc[h99_transcript_annotations["gene_id"] == gene_id_1].to_numpy()[0] + 1),95    "log_expression_2": np.log10(estimated_expression.loc[h99_transcript_annotations["gene_id"] == gene_id_2].to_numpy()[0] + 1),96    "run_accession": estimated_expression.columns})97chart_data = chart_data.merge(98    right = estimated_expression_meta,99    on = "run_accession")100 101with col3:102    st.text('') # help alignment with input box103    st.download_button(104        label = "Download data as TSV",105        data = chart_data.to_csv(sep ='\t').encode('utf-8'),106        file_name = f"CryptoCEN_expression_{gene_id_1}_vs_{gene_id_2}.tsv",107        mime = "text/csv")108 109 110st.markdown(f"""111#### Gene 1:112* *Gene ID*: [{gene_id_1}](https://fungidb.org/fungidb/app/record/gene/{gene_id_1})113{'* *Gene Symbol*:' + gene_symbol_1 if gene_symbol_1 is not None else ''}114* *Description*: {description_1}115* *Top [Co-Expressed Partners](https://huggingface.co/spaces/maomlab/CryptoCEN-TopHits?gene_id={gene_id_1})*116 117#### Gene 2:118* *Gene ID*: [{gene_id_2}](https://fungidb.org/fungidb/app/record/gene/{gene_id_2})119{'* *Gene Symbol*:' + gene_symbol_2 if gene_symbol_2 is not None else ''}120* *Description*: {description_2}121* *Top [Co-Expressed Partners](https://huggingface.co/spaces/maomlab/CryptoCEN-TopHits?gene_id={gene_id_2})*122""")123 124chart = (125   alt.Chart(126       chart_data,127       width = 750,128       height = 750)129   .mark_circle()130   .encode(131       x=alt.X("log_expression_1", title=f"Log10[{gene_id_1}+1] Expression"),132       y=alt.Y("log_expression_2", title=f"Log10[{gene_id_2}+1] Expression"),133       color=alt.Color("study_accession", title="Study Accession"),134       tooltip=["run_accession", "study_accession"]))135 136st.altair_chart(137    chart)138 139 140