<jupyter_start><jupyter_text>Acute myeloid leukemia mouse model heatmap_Author: <NAME> adapted this analysis based on this [R refine.bio-examples notebook](https://alexslemonade.github.io/refinebio-examples/03-rnaseq/clustering_rnaseq_01_heatmap.html)_   _Date: October 2021_   Purpose of the analysisIn this analysis, we will use this [acute myeloid leukemia sample dataset](https://www.refine.bio/experiments/SRP070849) from [Shih et al, 2017](https://pubmed.ncbi.nlm.nih.gov/28193779/) and pre-processed by [refinebio](https://www.refine.bio/).The data that we downloaded from [refine.bio](https://www.refine.bio/) for this analysis has 19 samples (obtained from 19 acute myeloid leukemia (AML) model mice), containing RNA-sequencing results for types of AML under controlled treatment conditions.This [dataset can be downloaded from this page on refine.bio](https://www.refine.bio/experiments/SRP070849). We will download it already [processed and quantile normalized](http://docs.refine.bio/en/latest/main_text.htmlrefine-bio-processed-refinebio-processedibadge).   Import libraries<jupyter_code># We need this for creating directories
import os
from pathlib import Path

# For session info
import session_info

# We need this for downloading data from refine.bio
import pyrefinebio

# We will use these for manipulating our data
import pandas as pd
import numpy as np
import random

# We will use this for making the heatmap and these for the color codes
import seaborn as sns

sns.set_theme(color_codes=True)

# Import custom color key function
import util.color_key as color_key

import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Set the seed<jupyter_code>random.seed(1234)<jupyter_output><empty_output><jupyter_text>Set up analysis folders<jupyter_code>dirnames = ["data", "results", "plots"]

# Create the directory if it doesn't exist
for dirname in dirnames:
    if not os.path.isdir(dirname):
        os.mkdir(dirname)<jupyter_output><empty_output><jupyter_text>This notebook will be downloading the data for the project ID 'SRP070849' from refine.bio. By switching the line below, the rest of the notebook where the the experiment ID is referenced will be updated.<jupyter_code># Declare project ID
id = "SRP070849"
# Define the file path to the data directory
data_dir = Path(f"data/{id}")

# Declare the file path to the gene expression matrix file
data_file = data_dir.joinpath(f"{id}.tsv")

# Declare the file path to the metadata file
# inside the directory saved as `data_dir`
metadata_file = data_dir.joinpath(f"metadata_{id}.tsv")<jupyter_output><empty_output><jupyter_text>Import and set up dataFirst we will download our dataset from refine.bio using their API client. But we will only download if the files don't already exist locally.<jupyter_code>if not data_file.exists() or not metadata_file.exists():

    print(f"Downloading {id} from refine.bio")

    pyrefinebio.create_token(agree_to_terms=True, save_token=False)

    pyrefinebio.download_dataset(
        data_dir.joinpath("dataset.zip"),
        "<EMAIL>",
        experiments=[f"{id}"],
        extract=True,
    )<jupyter_output><empty_output><jupyter_text>Now we can read in our expression data and metadata with `pandas`.<jupyter_code># Read in metadata TSV file
metadata = pd.read_csv(metadata_file, sep="\t")
# Read in data TSV file
expression_df = pd.read_csv(data_file, sep="\t")<jupyter_output><empty_output><jupyter_text>For our seaborn heatmap, we will what our gene IDs to be set as the index of our data frame.<jupyter_code># Set the Gene column as the index
expression_df.set_index("Gene", inplace=True)<jupyter_output><empty_output><jupyter_text>Let's take a look at the metadata object that we read into the R environment.<jupyter_code>metadata.head(5)<jupyter_output><empty_output><jupyter_text>Now let's ensure that the metadata and data are in the same sample order.<jupyter_code>print(metadata["refinebio_accession_code"].tolist() == expression_df.columns.tolist())<jupyter_output>True<jupyter_text>Choose genes of interestFor this example, we will sort genes by variance and select genes in the upper quartile, but there are many alternative criterion by which you may want to sort your genes, e.g. fold change, t-statistic, membership in a particular gene ontology, so on.<jupyter_code># Calculate the variance for each gene
expression_df["variance"] = expression_df.var(axis=1, skipna=True)

# Find the upper quartile for these data
upper_quartile = expression_df["variance"].quantile([0.90]).values

# Filter the data choosing only genes whose variances are in the upper quartile
df_by_var = expression_df[expression_df.variance > float(upper_quartile)]

# Drop the variance column we calculated
df_by_var = df_by_var.drop(columns="variance")
# Check how many genes we are left with
print(len(df_by_var.index))<jupyter_output>4125<jupyter_text>Write top variance genes to a TSV file<jupyter_code># Write df_by_var to tsv file
df_by_var.to_csv("results/top_90_var_genes.tsv", sep="\t")<jupyter_output><empty_output><jupyter_text>Prepare metadata for annotationFrom the accompanying [paper](https://pubmed.ncbi.nlm.nih.gov/28193779/), we know that the mice with `IDH2` mutant AML were treated with vehicle or AG-221 (the first small molecule in-vivo inhibitor of IDH2 to enter clinical trials) and the mice with `TET2` mutant AML were treated with vehicle or 5-Azacytidine (Decitabine, hypomethylating agent). [@Shih2017]<jupyter_code># Recode refinebio_title variable into something more useful
# Split wherever there is a '-' or '.'
exp_group = metadata["refinebio_title"].str.split("-|\\.", expand=True, n=1)
# Now we'll store that first column from the string split back in the metadata as a column
metadata["exp_group"] = exp_group[0]
# Print to check
print(metadata["exp_group"])
# Color code our two variables
refinebio_treatment_colors = color_key.make_color_key(metadata["refinebio_treatment"])
exp_groups_colors = color_key.make_color_key(metadata["exp_group"])<jupyter_output><empty_output><jupyter_text>Make them into one data frame we can provide to the heatmap for annotation.<jupyter_code>color_key_df = pd.concat(
    [refinebio_treatment_colors["color_key"], exp_groups_colors["color_key"]],
    axis=1,
    names=["treatment", "exp_group"],
)

color_key_df = color_key_df.set_index(df_by_var.columns)<jupyter_output><empty_output><jupyter_text>Print out our annotation data frame to make sure it looks right.<jupyter_code># Check that this is what we think it is
color_key_df<jupyter_output><empty_output><jupyter_text>Create and store the annotated heatmap object<jupyter_code>heatmap = sns.clustermap(
    df_by_var,
    cmap="mako",
    col_colors=color_key_df,
    dendrogram_ratio=(0, 0.2),
    cbar_pos=(-0.1, 0.2, 0.03, 0.5),
)

legend1 = color_key.make_legend(exp_groups_colors["color_key_dict"])
legend2 = color_key.make_legend(refinebio_treatment_colors["color_key_dict"])

# Make a legend 1
plt.legend(
    legend1,
    exp_groups_colors["color_key_dict"],
    title="Experiment \nGroup",
    bbox_to_anchor=(4, 1.3),
)

# Make a legend 2
heatmap.ax_row_dendrogram.legend(
    legend2,
    refinebio_treatment_colors["color_key_dict"],
    title="Treatment",
    bbox_to_anchor=(0, 1.37),
)

plt.show()<jupyter_output>/opt/conda/envs/reproducible-python/lib/python3.8/site-packages/seaborn/matrix.py:654: UserWarning: Clustering large matrix with scipy. Installing `fastcluster` may give better performance.
  warnings.warn(msg)<jupyter_text>Now we need to set up a legend for the color bar.<jupyter_code># Save to png
heatmap.savefig(os.path.join("plots", "aml_heatmap.png"))<jupyter_output><empty_output><jupyter_text>Print out session info<jupyter_code>session_info.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09667868177136972, 'fraction_numerical': 0.0148043254376931, 'mean_word_length': 6.024412296564195, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 69, 'https://': 7, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '221009', 'license': 'MIT', 'n_tokens_mistral': 2666, 'n_tokens_neox': 2474, 'n_words': 723, 'path': 'python-examples/python-heatmap-chapt-8/make_heatmap.ipynb', 'repo_name': 'jhudsl/Reproducibility_in_Cancer_Informatics', 'size': '7788'}