Merge pull request #237 from PNNL-CompBio/split_broad_sanger

Cell Line Datasets Separated, Improved Error Handling for all Build Scripts

Author: sgosline

.dockerignore

@@ -4,5 +4,4 @@ coderdata/
 dataSummary/
 docs/
 candle_bmd/
-schema/
-build/local/
+build/local/

build/beatAML/GetBeatAML.py

@@ -174,10 +174,8 @@ def retrieve_drug_info(compound_name):
         return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan
 
     data = response.json()
-    #print(data)
     if "PropertyTable" in data:
         properties = data["PropertyTable"]["Properties"][0]
-        #print(properties)
         pubchem_id = properties.get('CID',np.nan)
         canSMILES = properties.get("CanonicalSMILES", np.nan)
         isoSMILES = properties.get("IsomericSMILES", np.nan)
@@ -259,9 +257,6 @@ def merge_drug_info(d_df,drug_map):
     pd.DataFrame
         The merged dataframe containing combined drug information.
     """
-    #print(drug_map)
-    #print(d_df.columns)
-    #print(d_df)
     print(d_df['isoSMILES'].dtype, drug_map['isoSMILES'].dtype)
     d_df['isoSMILES'] = d_df['isoSMILES'].astype(str)
     drug_map['isoSMILES'] = drug_map['isoSMILES'].astype(str)
@@ -337,10 +332,9 @@ def add_improve_id(previous_df, new_df):
     """
     if not previous_df.empty and 'improve_drug_id' in previous_df.columns:
         id_list = [int(val.replace('SMI_', '')) for val in previous_df['improve_drug_id'].tolist() if pd.notnull(val) and val.startswith('SMI_')]
-        max_id = max(id_list) if id_list else 0  # Default to 0 if the list is empty
+        max_id = max(id_list) if id_list else 0
     else:
-        max_id = 0  # Default value if the DataFrame is empty or doesn't have the column
-    # max_id = max([int(val.replace('SMI_', '')) for val in previous_df['improve_drug_id'].tolist() if pd.notnull(val) and val.startswith('SMI_')])
+        max_id = 0
     # Identify isoSMILES in the new dataframe that don't exist in the old dataframe
     unique_new_smiles = set(new_df['isoSMILES']) - set(previous_df['isoSMILES'])
     # Identify rows in the new dataframe with isoSMILES that are unique and where improve_drug_id is NaN
@@ -370,24 +364,9 @@ def map_exp_to_improve(exp_path):#df,improve_map_file):
     pd.DataFrame
         Mapped dataframe with 'improve_sample_id' added and 'sample_id' removed.
     """
-    mapped_df = pd.read_csv(exp_path,sep='\t')        # Map sample_id to improve_sample_id
-    #mapped_df = pd.merge(df, improve[['other_id', 'improve_sample_id']], left_on='sample_id', right_on='other_id', how='left')
-    #mapped_df.drop(columns=['sample_id', 'other_id'], inplace=True)
-    #mapped_df.insert(0, 'improve_sample_id', mapped_df.pop('improve_sample_id'))
+    mapped_df = pd.read_csv(exp_path,sep='\t')
     mapped_df['source'] = 'synapse'
     mapped_df['study'] = 'BeatAML'
-    #mapped_df= mapped_df.rename(columns={'Drug':'improve_sample_id',
-    #                                     'IC50':'ic50',
-    #                          'EC50':'ec50',
-    #                         'EC50se':'ec50se',
-    #                         'Einf':'einf',
-     #                        'HS':'hs',
-      #                       'AAC1':'aac1',
-       #                      'AUC1':'auc1',
-        #                     'DSS1':'dss1',
-         #                    'R2fit':'r2fit'
-          #                   }
-          #          )
     return mapped_df
 
 
@@ -445,12 +424,21 @@ def map_and_combine(df, data_type, entrez_map_file, improve_map_file, map_file=N
         mapped_df.rename(columns={"hgvsc": "mutation"}, inplace=True)
         mapped_df.rename(columns={"labId": "sample_id"}, inplace=True)
         mapped_df.rename(columns={"Entrez_Gene_Id": "entrez_id"}, inplace=True)
-        
-    elif data_type == "mutation":
-        df = df[['dbgap_sample_id','hgvsc', 'hgvsp', 'gene', 'variant_classification','t_vaf', 'refseq', 'symbol']]
-        mapped_df = df.merge(genes, left_on='symbol', right_on='gene_symbol', how='left').reindex(
-                        columns=['hgvsc', 'entrez_id', "dbgap_sample_id","variant_classification"])
 
+        variant_mapping = {
+            'frameshift_variant': 'Frameshift_Variant',
+            'missense_variant': 'Missense_Mutation',
+            'stop_gained': 'Nonsense_Mutation',
+            'inframe_deletion': 'In_Frame_Del',
+            'protein_altering_variant': 'Protein_Altering_Variant',
+            'splice_acceptor_variant': 'Splice_Site',
+            'splice_donor_variant': 'Splice_Site',
+            'start_lost': 'Start_Codon_Del',
+            'inframe_insertion': 'In_Frame_Ins',
+            'stop_lost': 'Nonstop_Mutation'
+        }
+
+        mapped_df['variant_classification'] = mapped_df['variant_classification'].map(variant_mapping)
 
     elif data_type == "proteomics":
         mapped_ids['sampleID'] = mapped_ids['sampleID'].str.split('_').apply(lambda x: x[2])
@@ -473,7 +461,6 @@ def map_and_combine(df, data_type, entrez_map_file, improve_map_file, map_file=N
             inplace=True
         )
 
-
     mapped_df = pd.merge(mapped_df, improve[['other_id', 'improve_sample_id']], 
                          left_on='sample_id', 
                          right_on='other_id',
@@ -482,7 +469,7 @@ def map_and_combine(df, data_type, entrez_map_file, improve_map_file, map_file=N
     mapped_df['source'] = 'synapse'
     mapped_df['study'] = 'BeatAML'
 
-    final_dataframe = mapped_df.dropna()#pd.dropna(mapped_df,0)
+    final_dataframe = mapped_df.dropna()
     return final_dataframe
 
 
@@ -659,8 +646,6 @@ def generate_drug_list(drug_map_path,drug_path):
 
 
             t_df = pd.read_csv('tpm_'+transcriptomics_file, sep = '\t')
-           # t_df.index = t_df.stable_id#display_label
-#            t_df = t_df.iloc[:, 4:]
             t_df = t_df.reset_index().rename(columns={'stable_id': 'Gene'})
             t_df = pd.melt(t_df, id_vars=['Gene'], var_name='sample_id', value_name='transcriptomics')
             print(improve_map_file)
@@ -724,7 +709,5 @@ def generate_drug_list(drug_map_path,drug_path):
             exp_res = map_exp_to_improve(drug_path)
             exp_res.to_csv("/tmp/beataml_experiments.tsv", index=False, sep='\t')
 
-            #drug_map_path = retrieve_figshare_data("https://figshare.com/ndownloader/files/43112314?private_link=0ea222d9bd461c756fb0")
-
 #        print("Finished Pipeline")
 

build/beatAML/build_drugs.sh

@@ -1,2 +1,10 @@
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running GetBeatAML.py with token and drugFile $1"
 python GetBeatAML.py --token $SYNAPSE_AUTH_TOKEN --drugs --drugFile $1
+
+echo "Running build_drug_desc.py..."
 python build_drug_desc.py --drugtable /tmp/beataml_drugs.tsv --desctable /tmp/beataml_drug_descriptors.tsv.gz

build/beatAML/build_exp.sh

@@ -1 +1,7 @@
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running GetBeatAML.py with token and curSamples $1 and drugFile $2."
 python GetBeatAML.py --exp --token $SYNAPSE_AUTH_TOKEN --curSamples $1 --drugFile $2

build/beatAML/build_omics.sh

@@ -1 +1,7 @@
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running GetBeatAML.py with token, curSamples $2, and genes $1."
 python GetBeatAML.py --token $SYNAPSE_AUTH_TOKEN --omics --curSamples $2 --genes $1

build/beatAML/build_samples.sh

@@ -1 +1,7 @@
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running GetBeatAML.py with token and prevSamples $1."
 python GetBeatAML.py --token $SYNAPSE_AUTH_TOKEN --samples --prevSamples $1

build/broad_sanger/03a-nci60Drugs.py

@@ -122,9 +122,21 @@ def main():
     merged = pl.concat([mdf,namedf],how='horizontal').select(['SMILES','pubchem_id','nscid','lower_name'])
     melted = merged.melt(id_vars=['SMILES','pubchem_id'],value_vars=['nscid','lower_name']).select(['SMILES','pubchem_id','value']).unique()
     melted.columns = ['canSMILES','pubchem_id','chem_name']
-    if newdf.shape[0]>0:
-        newdf = newdf.join(melted,on='canSMILES',how='inner').select(res.columns)
-        res = pl.concat([res,newdf],how='vertical')
+
+    if newdf.shape[0] > 0:
+        res = res.with_columns([
+            pl.col("InChIKey").cast(pl.Utf8),
+            pl.col("formula").cast(pl.Utf8),
+            pl.col("weight").cast(pl.Utf8)
+        ])
+        newdf = newdf.with_columns([
+            pl.col("InChIKey").cast(pl.Utf8),
+            pl.col("formula").cast(pl.Utf8),
+            pl.col("weight").cast(pl.Utf8)
+        ])
+    
+        newdf = newdf.join(melted, on='canSMILES', how='inner').select(res.columns)
+        res = pl.concat([res, newdf], how='vertical')
     res.write_csv(opts.output,separator='\t')
 
 if __name__=='__main__':

build/broad_sanger/05_separate_datasets.py

@@ -0,0 +1,83 @@
+import gc
+import polars as pl 
+
+
+
+def main():
+
+    datasets_to_process = ["CCLE", "CTRPv2", "PRISM", "GDSCv1", "GDSCv2", "FIMM", "gCSI", "NCI60"]
+    omics_datatypes = ["transcriptomics","proteomics", "copy_number","mutations"] # csv 
+    samples_datatypes = ["samples"] #csv
+
+    drugs_datatypes = ["drugs", "drug_descriptors"] # tsv
+
+
+    dataset_sources = {
+        "CCLE": ["Broad"],
+        "CTRPv2": ["Broad"],
+        "PRISM": ["Broad"],
+        "GDSCv1": ["Sanger"],
+        "GDSCv2": ["Sanger"],
+        "FIMM": ["Broad"],
+        "gCSI": ["Broad"],  # gCSI generates its own omics data but it is comparable to CCLE. In future, retrive gCSI omics.
+        "NCI60": ["Broad"]
+    }
+
+    for dataset in datasets_to_process:
+        exp = pl.read_csv("broad_sanger_experiments.tsv", separator="\t") # Keeping memory down, so I will not be making copies.
+        exp = exp.filter(pl.col("study") == dataset)
+
+        # Extract information to separate out datasets
+        exp_improve_sample_ids = exp["improve_sample_id"].unique().to_list()
+        exp_improve_drug_ids = exp["improve_drug_id"].unique().to_list()
+
+        # Write Filtered Experiments File to TSV. Then delete it from memory.
+        exp_filename = f"/tmp/{dataset}_experiments.tsv".lower()
+        exp.write_csv(exp_filename, separator="\t")
+        del exp
+        gc.collect()
+
+
+        #Filter Samples files, write to file, delete from mem.
+        for samples in samples_datatypes:
+            samples_filename_in = f"broad_sanger_{samples}.csv"
+            samples_filename_out = f"/tmp/{dataset}_{samples}.csv".lower()
+            samples_df = pl.read_csv(samples_filename_in)
+            samples_df = samples_df.filter(pl.col("improve_sample_id").is_in(exp_improve_sample_ids))
+            samples_df.write_csv(samples_filename_out) #csv
+            del samples_df
+            gc.collect()
+
+        #One by one, filter other Omics files, write to file, delete from mem.
+        for omics in omics_datatypes:
+            omics_filename_in = f"broad_sanger_{omics}.csv"
+            omics_filename_out = f"/tmp/{dataset}_{omics}.csv".lower()
+            omics_df = pl.read_csv(omics_filename_in)
+            omics_df = omics_df.filter(pl.col("improve_sample_id").is_in(exp_improve_sample_ids))
+            omics_df = omics_df.filter(pl.col("source").is_in(dataset_sources[dataset]))
+            omics_df.write_csv(omics_filename_out) #csv
+            del omics_df
+            gc.collect()
+
+
+        #One by one, filter other Drugs files, write to file, delete from mem.
+        for drugs in drugs_datatypes:
+            drugs_filename_in = f"broad_sanger_{drugs}.tsv"
+            drugs_filename_out = f"/tmp/{dataset}_{drugs}.tsv".lower()
+            if drugs == "drug_descriptors":
+                drugs_df = pl.read_csv(drugs_filename_in,separator="\t",
+                                       dtypes={"improve_drug_id": pl.Utf8,
+                                                         "structural_descriptor": pl.Utf8,
+                                                         "descriptor_value": pl.Utf8}
+                                      )
+
+            else:
+                drugs_df = pl.read_csv(drugs_filename_in,separator="\t")
+
+            drugs_df = drugs_df.filter(pl.col("improve_drug_id").is_in(exp_improve_drug_ids))
+            drugs_df.write_csv(drugs_filename_out,separator="\t") #tsv
+            del drugs_df
+            gc.collect()
+            
+if __name__ == "__main__":
+    main()

build/broad_sanger/build_drugs.sh

@@ -1,3 +1,15 @@
-/opt/venv/bin/python 03a-nci60Drugs.py 
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running 03a-nci60Drugs.py..."
+/opt/venv/bin/python 03a-nci60Drugs.py
+
+echo "Running 03-createDrugFile.R..."
 Rscript 03-createDrugFile.R CTRPv2,GDSC,gCSI,PRISM,CCLE,FIMM
-/opt/venv/bin/python build_drug_desc.py --drugtable /tmp/broad_sanger_drugs.tsv --desctable /tmp/broad_sanger_drug_descriptors.tsv.gz
+
+echo "Running build_drug_desc.py..."
+/opt/venv/bin/python build_drug_desc.py \
+  --drugtable /tmp/broad_sanger_drugs.tsv \
+  --desctable /tmp/broad_sanger_drug_descriptors.tsv.gz

build/broad_sanger/build_exp.sh

@@ -1 +1,7 @@
-/opt/venv/bin/python 04-drug_dosage_and_curves.py --drugfile $2 --curSampleFile $1
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running 04-drug_dosage_and_curves.py with drugfile $2 and curSampleFile $1"
+/opt/venv/bin/python 04-drug_dosage_and_curves.py --drugfile $2 --curSampleFile $1