<jupyter_start><jupyter_text>Compare generic genesThe goal of this notebook is to compare the generic genes found using the same template experiment run two times and 2 different recount2 template experiments.<jupyter_code>%load_ext autoreload
%autoreload 2

import os
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from ponyo import utils
# Read in config variables
base_dir = os.path.abspath(os.path.join(os.getcwd(), "../"))
config_filename = os.path.abspath(
    os.path.join(base_dir, "configs", "config_human_general.tsv")
)

params = utils.read_config(config_filename)

local_dir = params["local_dir"]

project_id1 = "SRP012656"
project_id2 = "SRP061689"
# Get data directory containing gene summary data
data_dir = os.path.join(base_dir, "human_general_analysis")

# Get gene ranking files
gene_ranking_filename1 = os.path.join(
    data_dir, f"generic_gene_summary_{project_id1}.tsv"
)
gene_ranking_filename1_run2 = os.path.join(
    data_dir, f"generic_gene_summary_{project_id1}_run2.tsv"
)
gene_ranking_filename2 = os.path.join(
    data_dir, f"generic_gene_summary_{project_id2}.tsv"
)

# Get template data
template_filename1 = os.path.join(
    data_dir, "data", f"processed_recount2_template_{project_id1}.tsv"
)
template_filename2 = os.path.join(
    data_dir, "data", f"processed_recount2_template_{project_id2}.tsv"
)<jupyter_output><empty_output><jupyter_text>Correlation between rankings between same experimentHere we compare gene ranking after running SOPHIE 2 times using the same template experiment but different seeds.<jupyter_code># Load gene ranking
gene_ranking_summary1 = pd.read_csv(
    gene_ranking_filename1, sep="\t", index_col=0, header=0
)
gene_ranking_summary1_run2 = pd.read_csv(
    gene_ranking_filename1_run2, sep="\t", index_col=0, header=0
)
# Get simulated ranking
gene_ranking1 = gene_ranking_summary1["Rank (simulated)"].rename("Rank 1")
gene_ranking1_run2 = gene_ranking_summary1_run2["Rank (simulated)"].rename("Rank 2")

# Combine ranking
gene_ranking_same_combined = pd.concat([gene_ranking1, gene_ranking1_run2], axis=1)
print(gene_ranking_same_combined.shape)
gene_ranking_same_combined.head()
# Check for NAs
gene_ranking_same_combined[pd.isnull(gene_ranking_same_combined).any(axis=1)]
# Plot correlation between ranking
fig = sns.jointplot(
    data=gene_ranking_same_combined,
    x="Rank 1",
    y="Rank 2",
    kind="hex",
    marginal_kws={"color": "white"},
)

fig.set_axis_labels(
    f"Ranking in {project_id1}",
    f"Ranking in {project_id1} different runs",
    fontsize=14,
    fontname="Verdana",
)
plt.colorbar()

output_figure_filename = "concordance_between_same_recount2_templates.svg"
fig.savefig(
    output_figure_filename,
    format="svg",
    bbox_inches="tight",
    transparent=True,
    pad_inches=0,
    dpi=300,
)<jupyter_output>/home/alexandra/anaconda3/envs/generic_expression/lib/python3.7/site-packages/matplotlib/font_manager.py:1238: UserWarning: findfont: Font family ['Verdana'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))<jupyter_text>**Takeaway:*** Running SOPHIE twice using the same template experiment will generate 2 different sets of simulated experiments.* Since the template experiment is the same, these 2 sets of simulated experiments will have the same experimental design structure/biological context* As expected, the concordance is very high especially for high ranked and low ranked genes. The genes in the middle rank are more sensitive to changes so you don't get as clear of a signal compared to the extreme ranked genes.  Correlation between rankings between 2 different experimentsHere we compare gene ranking generated by SOPHIE using 2 different template experiments.<jupyter_code># Load gene ranking
gene_ranking_summary2 = pd.read_csv(
    gene_ranking_filename2, sep="\t", index_col=0, header=0
)
# Get simulated ranking
gene_ranking1 = gene_ranking_summary1["Rank (simulated)"].rename("Rank 1")
gene_ranking2 = gene_ranking_summary2["Rank (simulated)"].rename("Rank 2")

# Combine ranking
gene_ranking_diff_combined = pd.concat([gene_ranking1, gene_ranking2], axis=1)
print(gene_ranking_diff_combined.shape)
gene_ranking_diff_combined.head()
# Check for NAs
gene_ranking_diff_combined[pd.isnull(gene_ranking_diff_combined).any(axis=1)]
# Plot correlation between ranking
fig = sns.jointplot(
    data=gene_ranking_diff_combined,
    x="Rank 1",
    y="Rank 2",
    kind="hex",
    marginal_kws={"color": "white"},
)

fig.set_axis_labels(
    f"Ranking in {project_id1}",
    f"Ranking in {project_id2}",
    fontsize=14,
    fontname="Verdana",
)

plt.colorbar()

output_figure_filename = "concordance_between_diff_recount2_templates.svg"
fig.savefig(
    output_figure_filename,
    format="svg",
    bbox_inches="tight",
    transparent=True,
    pad_inches=0,
    dpi=300,
)<jupyter_output><empty_output><jupyter_text>**Takeaway:*** Looks like there is good concordance between highly ranked genes (i.e. generic genes)* By comparison if we run SOPHIE using two different template experiments, there are genes in the off-diagonal regions that might indicate that there are generic within the given context of the specific experiment.* In general, the genes in the middle rank are more sensitive to changes so you don't get as clear of a signal compared to the highest rank genes.  Examine gene expression data<jupyter_code># Read expression data
template_1 = pd.read_csv(template_filename1, sep="\t", index_col=0, header=0)
template_2 = pd.read_csv(template_filename2, sep="\t", index_col=0, header=0)
# Get concordance genes
concordant_genes = list(
    gene_ranking_diff_combined[
        (gene_ranking_diff_combined["Rank 1"] > 15000)
        & (gene_ranking_diff_combined["Rank 2"] > 15000)
    ].index
)

# Get disconcordant genes
discordant_genes = set(gene_ranking_diff_combined.index).difference(concordant_genes)
# Distribution of concordant genes in template experiment 1
template1_mean = template_1.mean()

print(
    "Percent concordant genes with 0 expression in template 1:",
    len(template1_mean[concordant_genes].loc[template1_mean[concordant_genes] == 0])
    / len(template1_mean[concordant_genes]),
)

print(
    "Percent nonzero concordant genes in template 1:",
    len(
        template1_mean[concordant_genes].loc[
            (template1_mean[concordant_genes] > 0)
            & (template1_mean[concordant_genes] < 1000)
        ]
    )
    / len(template1_mean[concordant_genes]),
)

f1 = sns.distplot(template_1.mean()[concordant_genes], kde=False)
f1.set_title(f"Expression of concordant genes in {project_id1}")
f1.set_xlabel("log(gene expression)")
f1.set_ylabel("log(count)")
f1.set(xscale="log", yscale="log")
# Distribution of concordant genes in template experiment 2
template2_mean = template_2.mean()
print(
    "Percent concordant genes with 0 expression in template 2:",
    len(template2_mean[concordant_genes].loc[template2_mean[concordant_genes] == 0])
    / len(template2_mean[concordant_genes]),
)

print(
    "Percent nonzero concordant genes in template 2:",
    len(
        template2_mean[concordant_genes].loc[
            (template2_mean[concordant_genes] > 0)
            & (template2_mean[concordant_genes] < 1000)
        ]
    )
    / len(template2_mean[concordant_genes]),
)

# There are more 0 expressed genes in this template experiment
f2 = sns.distplot(template_2.mean()[concordant_genes], kde=False)
f2.set_title(f"Expression of concordant genes in {project_id2}")
f2.set_xlabel("log(gene expression)")
f2.set_ylabel("log(count)")
f2.set(xscale="log", yscale="log")
# Distribution of discordant gense in template experiment 1
template1_mean = template_1.mean()

print(
    "Percent discordant genes with 0 expression in template 1:",
    len(template1_mean[discordant_genes].loc[template1_mean[discordant_genes] == 0])
    / len(template1_mean[discordant_genes]),
)

print(
    "Percent nonzero discordant genes in template 1:",
    len(
        template1_mean[discordant_genes].loc[
            (template1_mean[discordant_genes] > 0)
            & (template1_mean[discordant_genes] < 1000)
        ]
    )
    / len(template1_mean[discordant_genes]),
)

print(
    len(template1_mean[discordant_genes].loc[template1_mean[discordant_genes] > 0])
    / len(template1_mean[discordant_genes])
)
f3 = sns.distplot(template_1.mean()[discordant_genes], kde=False)
f3.set_title(f"Expression of discordant genes in {project_id1}")
f3.set_xlabel("log(gene expression)")
f3.set_ylabel("log(count)")
f3.set(xscale="log", yscale="log")
# Distribution of discordant genes in template experiment 2
template2_mean = template_2.mean()

print(
    "Percent discordant genes with 0 expression in template 2:",
    len(template2_mean[discordant_genes].loc[template2_mean[discordant_genes] == 0])
    / len(template2_mean[discordant_genes]),
)

print(
    "Percent nonzero discordant genes in template 2:",
    len(
        template2_mean[discordant_genes].loc[
            (template2_mean[discordant_genes] > 0)
            & (template2_mean[discordant_genes] < 1000)
        ]
    )
    / len(template2_mean[discordant_genes]),
)

f4 = sns.distplot(template_2.mean()[discordant_genes], kde=False)
f4.set_title(f"Expression of discordant genes in {project_id2}")
f4.set_xlabel("log(gene expression)")
f4.set_ylabel("log(count)")
f4.set(xscale="log", yscale="log")<jupyter_output>Percent discordant genes with 0 expression in template 2: 0.059811578033965536
Percent nonzero discordant genes in template 2: 0.5682409817776125
{'fraction_non_alphanumeric': 0.10355648535564854, 'fraction_numerical': 0.028242677824267783, 'mean_word_length': 5.460135135135135, 'pattern_counts': {'":': 2, '<': 19, '<?xml version=': 0, '>': 22, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '397663', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3394, 'n_tokens_neox': 3134, 'n_words': 729, 'path': 'compare_experiments/compare_generic_genes.ipynb', 'repo_name': 'jjc2718/generic-expression-patterns', 'size': '9560'}