Pseudomonas validation

configs/config_pseudomonas_33245.tsv

@@ -8,6 +8,9 @@ normalized_compendium_filename	"/home/alexandra/Documents/Data/Generic_expressio
 shared_genes_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns/shared_genes_pseudomonas.pickle"
 scaler_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns/scaler_transform_pseudomonas.pickle"
 rank_genes_by	"logFC"
+reference_gene_filename	"GAPE_proportions.txt"
+reference_gene_name_col	"gene id"
+reference_rank_col	"prop DEGs"
 pathway_DB_filename	"https://raw.githubusercontent.com/greenelab/adage/master/Node_interpretation/pseudomonas_KEGG_terms.txt"
 gsea_statistic	'log2FoldChange'
 rank_pathways_by	"padj"

configs/config_pseudomonas_pao1.tsv

@@ -0,0 +1,28 @@
+local_dir	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/"
+dataset_name	"pseudomonas_analysis"
+raw_template_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/raw_pseudomonas_template_data.tsv"
+processed_template_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/processed_pseudomonas_template_data.tsv"
+raw_compendium_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/raw_pseudomonas_compendium_data.tsv"
+processed_compendium_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/processed_pseudomonas_compendium_data.tsv"
+normalized_compendium_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/normalized_pseudomonas_compendium_data.tsv"
+shared_genes_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/shared_genes_pseudomonas.pickle"
+scaler_filename	"/home/alexandra/Documents/Data/Generic_expression_patterns_pao1/scaler_transform_pseudomonas.pickle"
+rank_genes_by	"logFC"
+reference_gene_filename	"GAPE_proportions.txt"
+reference_gene_name_col	"gene id"
+reference_rank_col	"prop DEGs"
+pathway_DB_filename	"https://raw.githubusercontent.com/greenelab/adage/master/Node_interpretation/pseudomonas_KEGG_terms.txt"
+gsea_statistic	'log2FoldChange'
+rank_pathways_by	"padj"
+NN_architecture	"NN_2500_30_pao1"
+learning_rate	0.001
+batch_size	10
+epochs	100
+kappa	0.01
+intermediate_dim	2500
+latent_dim	30
+epsilon_std	1.0
+validation_frac	0.25
+project_id	"E-GEOD-33245"
+metadata_colname	'ml_data_source'
+num_simulated	25

generic_expression_patterns_modules/process.py

@@ -358,15 +358,12 @@ def map_recount2_data(
                 ofh.write(sample_id + "\t" + "\t".join(output_values) + "\n")
 
 
-# TO DO: decide if not needed
 def process_raw_template_pseudomonas(
     processed_compendium_filename,
     project_id,
     dataset_name,
     metadata_colname,
-    sample_id_metadata_filename,
     raw_template_filename,
-    processed_template_filename,
 ):
     """
     Create processed pseudomonas template data file based on
@@ -388,23 +385,6 @@ def process_raw_template_pseudomonas(
 
     template_data.to_csv(raw_template_filename, sep="\t")
 
-    sample_ids_to_drop = set()
-    if os.path.exists(sample_id_metadata_filename):
-        # Read in metadata and get samples to be dropped:
-        metadata = pd.read_csv(
-            sample_id_metadata_filename, sep="\t", header=0, index_col=0
-        )
-        sample_ids_to_drop = set(metadata[metadata["processing"] == "drop"].index)
-
-    # Write the processed pseudomonas template output file on disk
-    with open(raw_template_filename) as ifh, open(
-        processed_template_filename, "w"
-    ) as ofh:
-        for idx, line in enumerate(ifh):
-            sample_id = line.split("\t")[0]
-            if idx == 0 or sample_id not in sample_ids_to_drop:
-                ofh.write(line)
-
 
 def normalize_compendium(
     mapped_filename, normalized_filename, scaler_filename,
@@ -444,17 +424,21 @@ def process_raw_compendium_pseudomonas(
     """
     Create processed pseudomonas compendium data file based on raw compendium
     data file (`raw_filename`), and normalize the processed compendium.
+
+    Note: This function was designed to processed data from the pseudomonas
+    compendium defined in the ADAGE paper
+    (https://msystems.asm.org/content/1/1/e00025-15).
     """
 
     # Create processed pseudomonas compendium data file
     raw_compendium = pd.read_csv(raw_filename, header=0, index_col=0, sep="\t")
 
-    if raw_compendium.shape != (950, 5549):
+    if raw_compendium.shape[1] != 5549:
         processed_compendium = raw_compendium.T
     else:
         processed_compendium = raw_compendium
 
-    assert processed_compendium.shape == (950, 5549)
+    assert processed_compendium.shape[1] == 5549
 
     # Save transformed compendium data
     processed_compendium.to_csv(processed_filename, sep="\t")

generic_expression_patterns_modules/ranking.py

@@ -19,6 +19,7 @@
 from scipy import stats
 import pandas as pd
 import seaborn as sns
+import matplotlib.pyplot as plt
 from sklearn.preprocessing import MinMaxScaler
 
 
@@ -422,7 +423,9 @@ def get_shared_rank_scaled(
     ]
 
     # Get correlation
-    r, p, ci_low, ci_high = spearman_ci(0.95, shared_rank_scaled_df, 1000, data_type)
+    r, p, ci_low, ci_high = spearman_ci(
+        0.95, shared_rank_scaled_df, 1000, ref_rank_col, data_type
+    )
 
     correlations = {"r": r, "p": p, "ci_low": ci_low, "ci_high": ci_high}
 
@@ -467,9 +470,15 @@ def compare_gene_ranking(
         marginal_kws={"color": "white"},
     )
 
-    fig.set_axis_labels(
-        "SOPHIE", "DE prior (Crow et. al. 2019)", fontsize=14, fontname="Verdana"
-    )
+    if ref_rank_col == "DE_Prior_Rank":
+        fig.set_axis_labels(
+            "SOPHIE", "DE prior (Crow et. al. 2019)", fontsize=14, fontname="Verdana"
+        )
+    elif ref_rank_col == "prop DEGs":
+        fig.set_axis_labels(
+            "SOPHIE", "GAPE (Stanton lab, 2020)", fontsize=14, fontname="Verdana"
+        )
+    plt.colorbar()
 
     fig.savefig(
         output_figure_filename,
@@ -480,7 +489,7 @@ def compare_gene_ranking(
         dpi=300,
     )
 
-    return correlations
+    return correlations, shared_gene_rank_scaled_df
 
 
 def compare_pathway_ranking(summary_df, reference_filename, output_figure_filename):
@@ -585,7 +594,7 @@ def add_pseudomonas_gene_name_col(summary_gene_ranks, base_dir):
     return summary_gene_ranks
 
 
-def spearman_ci(ci, gene_rank_df, num_permutations, data_type):
+def spearman_ci(ci, gene_rank_df, num_permutations, ref_rank_col, data_type):
     """
     Returns spearman correlation score and confidence interval
 
@@ -597,11 +606,13 @@ def spearman_ci(ci, gene_rank_df, num_permutations, data_type):
         Dataframe containing the our rank and Crow et. al. rank
     num_permutations: int
         The number of permutations to estimate the confidence interval
+    ref_rank_col: str
+        Name of column header containing reference ranks of genes
     data_type: 'DE' or 'GSA'
     """
     if data_type.lower() == "de":
         r, p = stats.spearmanr(
-            gene_rank_df["Rank (simulated)"], gene_rank_df["DE_Prior_Rank"]
+            gene_rank_df["Rank (simulated)"], gene_rank_df[ref_rank_col]
         )
     elif data_type.lower() == "gsa":
         r, p = stats.spearmanr(
@@ -615,7 +626,7 @@ def spearman_ci(ci, gene_rank_df, num_permutations, data_type):
 
         if data_type.lower() == "de":
             r_perm, p_perm = stats.spearmanr(
-                sample["Rank (simulated)"], sample["DE_Prior_Rank"]
+                sample["Rank (simulated)"], sample[ref_rank_col]
             )
         elif data_type.lower() == "gsa":
             r_perm, p_perm = stats.spearmanr(

human_cancer_analysis/2_identify_generic_genes_pathways.ipynb

@@ -1151,7 +1151,7 @@
     "\n",
     "figure_filename = f\"gene_ranking_{col_to_rank_genes}.svg\"\n",
     "\n",
-    "ranking.compare_gene_ranking(\n",
+    "corr, shared_ranking = ranking.compare_gene_ranking(\n",
     "    summary_gene_ranks,\n",
     "    DE_prior_filename,\n",
     "    ref_gene_col,\n",

human_cancer_analysis/nbconverted/2_identify_generic_genes_pathways.py

@@ -517,7 +517,7 @@ def shift_template_experiment_with_metadatafile(
 
 figure_filename = f"gene_ranking_{col_to_rank_genes}.svg"
 
-ranking.compare_gene_ranking(
+corr, shared_ranking = ranking.compare_gene_ranking(
     summary_gene_ranks,
     DE_prior_filename,
     ref_gene_col,

human_general_analysis/2_identify_generic_genes_pathways.ipynb

@@ -1046,7 +1046,7 @@
     "\n",
     "figure_filename = f\"gene_ranking_{col_to_rank_genes}.svg\"\n",
     "\n",
-    "ranking.compare_gene_ranking(\n",
+    "corr, shared_ranking = ranking.compare_gene_ranking(\n",
     "    summary_gene_ranks,\n",
     "    DE_prior_filename,\n",
     "    ref_gene_col,\n",

human_general_analysis/nbconverted/2_identify_generic_genes_pathways.py

@@ -310,7 +310,7 @@
 
 figure_filename = f"gene_ranking_{col_to_rank_genes}.svg"
 
-ranking.compare_gene_ranking(
+corr, shared_ranking = ranking.compare_gene_ranking(
     summary_gene_ranks,
     DE_prior_filename,
     ref_gene_col,

other_enrichment_methods/nbconverted/1_simulate_data.py

@@ -1,16 +1,15 @@
-
 # coding: utf-8
 
 # # Simulate gene expression data
-# 
+#
 # This notebook simulates gene expression data that can then be plugged into different enrichment methods
 
 # In[1]:
 
 
-get_ipython().run_line_magic('load_ext', 'autoreload')
-get_ipython().run_line_magic('load_ext', 'rpy2.ipython')
-get_ipython().run_line_magic('autoreload', '2')
+get_ipython().run_line_magic("load_ext", "autoreload")
+get_ipython().run_line_magic("load_ext", "rpy2.ipython")
+get_ipython().run_line_magic("autoreload", "2")
 
 import os
 import sys
@@ -19,6 +18,7 @@
 import pickle
 
 from rpy2.robjects import pandas2ri
+
 pandas2ri.activate()
 
 from ponyo import utils, simulate_expression_data
@@ -45,36 +45,28 @@
 
 # Load params
 local_dir = params["local_dir"]
-dataset_name = params['dataset_name']
-NN_architecture = params['NN_architecture']
-num_runs = params['num_simulated']
-project_id = params['project_id']
-metadata_col_id = params['metadata_colname']
-mapped_template_filename = params['mapped_template_filename']
-processed_template_filename = params['processed_template_filename']
-normalized_compendium_filename = params['normalized_compendium_filename']
-scaler_filename = params['scaler_filename']
-col_to_rank_genes = params['rank_genes_by']
-col_to_rank_pathways = params['rank_pathways_by']
-statistic = params['gsea_statistic']
-count_threshold = params['count_threshold']
+dataset_name = params["dataset_name"]
+NN_architecture = params["NN_architecture"]
+num_runs = params["num_simulated"]
+project_id = params["project_id"]
+metadata_col_id = params["metadata_colname"]
+mapped_template_filename = params["mapped_template_filename"]
+processed_template_filename = params["processed_template_filename"]
+normalized_compendium_filename = params["normalized_compendium_filename"]
+scaler_filename = params["scaler_filename"]
+col_to_rank_genes = params["rank_genes_by"]
+col_to_rank_pathways = params["rank_pathways_by"]
+statistic = params["gsea_statistic"]
+count_threshold = params["count_threshold"]
 
 # Load metadata file with grouping assignments for samples
 sample_id_metadata_filename = os.path.join(
-    base_dir,
-    dataset_name,
-    "data",
-    "metadata",
-    f"{project_id}_process_samples.tsv"
+    base_dir, dataset_name, "data", "metadata", f"{project_id}_process_samples.tsv"
 )
 
 # Load metadata file with grouping assignments for samples
 metadata_filename = os.path.join(
-    base_dir,
-    dataset_name,
-    "data",
-    "metadata",
-    f"{project_id}_groups.tsv"
+    base_dir, dataset_name, "data", "metadata", f"{project_id}_groups.tsv"
 )
 
 # Load pickled file
@@ -87,22 +79,18 @@
 
 # Output files
 gene_summary_filename = os.path.join(
-    base_dir, 
-    dataset_name, 
-    f"generic_gene_summary_{project_id}.tsv"
+    base_dir, dataset_name, f"generic_gene_summary_{project_id}.tsv"
 )
 
 pathway_summary_filename = os.path.join(
-    base_dir, 
-    dataset_name, 
-    f"generic_pathway_summary_{project_id}.tsv"
+    base_dir, dataset_name, f"generic_pathway_summary_{project_id}.tsv"
 )
 
 
 # ### Simulate experiments using selected template experiment
-# 
+#
 # Workflow:
-# 
+#
 # 1. Get the gene expression data for the selected template experiment
 # 2. Encode this experiment into a latent space using the trained VAE model
 # 3. Linearly shift the encoded template experiment in the latent space
@@ -113,11 +101,11 @@
 
 
 # Simulate multiple experiments
-# This step creates the following files in "<local_dir>/pseudo_experiment/" directory:           
+# This step creates the following files in "<local_dir>/pseudo_experiment/" directory:
 #   - selected_simulated_data_SRP012656_<n>.txt
 #   - selected_simulated_encoded_data_SRP012656_<n>.txt
 #   - template_normalized_data_SRP012656_test.txt
-# in which "<n>" is an integer in the range of [0, num_runs-1] 
+# in which "<n>" is an integer in the range of [0, num_runs-1]
 os.makedirs(os.path.join(local_dir, "pseudo_experiment"), exist_ok=True)
 for run_id in range(num_runs):
     simulate_expression_data.shift_template_experiment(
@@ -129,12 +117,12 @@
         scaler,
         local_dir,
         base_dir,
-        run_id
+        run_id,
     )
 
 
 # ## Process template and simulated experiments
-# 
+#
 # * Remove samples not required for comparison
 # * Make sure ordering of samples matches metadata for proper comparison
 # * Make sure values are cast as integers for using DESeq
@@ -145,31 +133,29 @@
 
 if not os.path.exists(sample_id_metadata_filename):
     sample_id_metadata_filename = None
-    
+
 stats.process_samples_for_DESeq(
-        mapped_template_filename,
-        metadata_filename,
-        processed_template_filename,
-        count_threshold,
-        sample_id_metadata_filename,
-    )
+    mapped_template_filename,
+    metadata_filename,
+    processed_template_filename,
+    count_threshold,
+    sample_id_metadata_filename,
+)
 
 for i in range(num_runs):
     simulated_filename = os.path.join(
-        local_dir,
-        "pseudo_experiment",
-        f"selected_simulated_data_{project_id}_{i}.txt"
+        local_dir, "pseudo_experiment", f"selected_simulated_data_{project_id}_{i}.txt"
     )
     out_simulated_filename = os.path.join(
         local_dir,
         "pseudo_experiment",
-        f"selected_simulated_data_{project_id}_{i}_processed.txt"
+        f"selected_simulated_data_{project_id}_{i}_processed.txt",
     )
     stats.process_samples_for_DESeq(
         simulated_filename,
         metadata_filename,
         out_simulated_filename,
         count_threshold,
         sample_id_metadata_filename,
-)
+    )