Merge remote-tracking branch 'origin/main' into cdc_organoids

Author: alexandriai168

README.md

@@ -1,5 +1,7 @@
 ## Cancer Omics Drug Experiment Response Dataset 
 
+
+
 There is a recent explosion of deep learning algorithms that to tackle the computational problem of predicting drug treatment outcome from baseline molecular measurements. To support this,we have built a benchmark dataset that harmonizes diverse datasets to better assess algorithm performance.
 
 This package collects diverse sets of paired molecular datasets with corresponding drug sensitivity data. All data here is reprocessed and standardized so it can be easily used as a benchmark dataset for the 

build/beatAML/GetBeatAML.py

@@ -653,7 +653,7 @@ def generate_drug_list(drug_map_path,drug_path):
             # New Transcriptomics Data
             print("Starting Transcriptomics Data")
             ##first run conversion tool
-            os.system("python tpmFromCounts.py --counts "+transcriptomics_file)
+            os.system("python tpmFromCounts.py --counts {} --out_file {}".format(transcriptomics_file,'tpm_'+transcriptomics_file))
 
 
             t_df = pd.read_csv('tpm_'+transcriptomics_file, sep = '\t')

build/bladderpdo/00_createBladderPDOSampleFile.py

@@ -0,0 +1,50 @@
+import synapseclient
+import pandas as pd
+import numpy as np
+import argparse
+import os
+
+
+def get_bladder_pdo_samples(synLoginObject, maxval):
+    
+    # download from Synapse..
+    samples_syn = synLoginObject.get('syn64765486')
+    # and read the file
+    samples_df = pd.read_csv(samples_syn.path, sep="\t")
+
+    samples = samples_df[['Sample ID', 'Patient ID', 'Cancer Type Detailed', 'Sample Class']]
+    samples = samples.rename({"Sample ID" : 'other_id', 'Patient ID' : 'common_name', 'Cancer Type Detailed': 'cancer_type', 'Sample Class' : 'model_type'}, axis=1)
+
+    samples.loc[:,['species']] = 'Homo sapiens(Human)'
+    samples.loc[:,['other_id_source']] = 'Synapse'
+    samples.loc[:,['other_names'] ]= ''
+    samples.loc[:,['cancer_type']]=samples['cancer_type'].str.lower()
+    samples.loc[:, ['model_type']] = samples['model_type'].str.lower()
+
+    samples['improve_sample_id'] = range(maxval+1, maxval+1+samples.shape[0])
+
+    return samples
+
+
+if __name__ == "__main__":
+    
+    parser = argparse.ArgumentParser(description="This script handles downloading, processing and formatting of sample files for the Sarcoma PDO project into a single samplesheet")
+    
+    parser.add_argument('-t', '--token', type=str, help='Synapse Token')
+
+    parser.add_argument("-p", '--prevSamples', nargs="?", type=str, default ="", const  = "", help = "Use this to provide previous sample file, will run sample file generation")
+
+    args = parser.parse_args()
+   
+    print("Logging into Synapse")
+    PAT = args.token
+    synObject = synapseclient.login(authToken=PAT)
+
+    if (args.prevSamples):
+        prev_max_improve_id = max(pd.read_csv(args.prevSamples).improve_sample_id)
+    else: 
+        prev_max_improve_id = 0
+
+    bladder_pdo_samples = get_bladder_pdo_samples(synObject, prev_max_improve_id)
+
+    bladder_pdo_samples.to_csv("/tmp/bladderpdo_samples.csv", index=False)

build/bladderpdo/01_createBladderPDOOmicsFiles.py

@@ -0,0 +1,132 @@
+import synapseclient
+import pandas as pd
+import numpy as np
+import argparse
+import os
+import wget
+import gzip
+import subprocess
+import math
+
+def get_copy_call(a):
+    """
+    Heler Function - Determine copy call for a value.
+    """
+
+    if a is None:
+        return float('nan')
+
+    if math.isnan(a):
+        return float('nan')
+
+    a_val = math.log2(float(a)+0.000001)
+    if a_val < 0.5210507:
+        return 'deep del'
+    elif a_val < 0.7311832:
+        return 'het loss'
+    elif a_val < 1.214125:
+        return 'diploid'
+    elif a_val < 1.422233:
+        return 'gain'
+    else:
+        return 'amp'
+
+    return pd.Series([get_copy_call(a) for a in arr])
+
+def get_bladder_pdo_transcriptomics(GEO_id_link_table, samples, genes):
+
+    bladderpdo_url ='https://ftp.ncbi.nlm.nih.gov/geo/series/GSE103nnn/GSE103990/suppl/GSE103990_Normalized_counts.txt.gz'
+    transcriptomic_txt = wget.download(bladderpdo_url)
+    transcriptomics = pd.read_csv(transcriptomic_txt, compression='gzip', sep="\t")
+    subprocess.call (["/usr/bin/Rscript", "--vanilla", "obtainGSMidLink.R"])
+
+    GEO_ids_link = pd.read_csv("./gsmlinkDf.csv")
+    fpkm_totals = transcriptomics.iloc[:, 1:43].sum()
+    transcriptomics.iloc[:, 1:43] = transcriptomics.iloc[:, 1:43].div(fpkm_totals).mul(1e6)
+    transcriptomics['ensembl'] = transcriptomics['Unnamed: 0'].str.split("_", expand=True)[0]
+    mapped_df = transcriptomics.merge(genes[['entrez_id', 'other_id']].drop_duplicates(), left_on='ensembl', right_on='other_id', how='left')
+    # transform data to long format
+
+    mapped_df.drop('other_id', axis=1)
+    value_variables = transcriptomics.columns[transcriptomics.columns.str.contains("M")]
+    melted_txomics = mapped_df.melt(id_vars = "entrez_id", value_vars = value_variables, var_name='sample_name')
+    # use info from GEO to get Sample IDS
+    txomics_with_GEOid = melted_txomics.merge(GEO_ids_link, how = 'left', left_on = "sample_name", right_on='RNAid')
+    # use samplesheet to link sample_ids to improve ids
+    txomics_with_GEOid['sampleid'] = txomics_with_GEOid['sampleid'].str.replace("org", "Organoid_")
+    txomics_with_GEOid['sampleid'] = txomics_with_GEOid['sampleid'].str.replace("tumor", "Tumor")
+    txomics_with_improveid = txomics_with_GEOid.merge(samples, left_on="sampleid", right_on="other_id", how="left")
+    final_transcriptomics = txomics_with_improveid[['entrez_id', 'value', 'improve_sample_id']]
+    final_transcriptomics['source'] = "Gene Expression Omnibus"
+    final_transcriptomics['study'] = "Lee etal 2018 Bladder PDOs"
+    final_transcriptomics.rename({'value' : 'transcriptomics' })
+    # remove duplicates
+    toreturn = final_transcriptomics.drop_duplicates()
+
+    return toreturn
+
+def get_bladder_pdo_mutations(synObject, samples, genes):
+    print(samples.head)
+    mutations = synObject.get("syn64765525")
+    mutations_df = pd.read_csv(mutations.path, sep='\t')
+    mutations_df['mutation'] = mutations_df['HGVSc'].str.split(":", expand=True)[1]
+    #samplesheet = pd.read_csv(samples)
+    selectioncols_mutations = mutations_df[['Entrez_Gene_Id',"Variant_Classification", "Tumor_Sample_Barcode", "mutation"]]
+    merged_mutations = selectioncols_mutations.merge(samples, left_on="Tumor_Sample_Barcode", right_on="other_id", how="left")
+    merged_mutations_renamed = merged_mutations.rename({"Entrez_Gene_Id" : 'entrez_id', 'Variant_Classification' : "variant_classification"}, axis=1)
+    print(merged_mutations_renamed.head)
+    final_mutations = merged_mutations_renamed[['entrez_id', "mutation", "variant_classification", "improve_sample_id"]]
+    final_mutations['study'] = "Lee etal 2018 Bladder PDOs"
+    print(final_mutations.head)
+    return final_mutations
+
+def get_bladder_pdo_copynumber(synObject, samples, genes):
+    segfile = synObject.get("syn64765499")
+    segfile_df = pd.read_csv(segfile.path, sep='\t')
+
+    segfile_df.to_csv("bladder_segfile.csv")
+    subprocess.call (["/usr/bin/Rscript", "--vanilla", "CNV-segfile-annotation.R", "bladder_segfile.csv", "bladder_annotated_segfile.csv"])
+    copynumber = pd.read_csv("bladder_annotated_segfile.csv")
+    copynumber['copy_number'] = np.exp2(copynumber['score'].div(2))*2
+    copynumber['copy_call'] = [get_copy_call(a) for a in copynumber['copy_number']]
+    copynumber_with_improveids = copynumber.merge(samples, left_on='ID', right_on = 'other_id', how='left')
+    copynumber_with_correct_colnames = copynumber_with_improveids.rename({"ENTREZID":'entrez_id'}, axis=1)
+    final_copynumber = copynumber_with_correct_colnames[['entrez_id', 'improve_sample_id', 'copy_number', 'copy_call']]
+    final_copynumber['source'] = "Synapse"
+    final_copynumber['study'] = "Lee etal 2018 Bladder PDOs"
+
+    return final_copynumber
+
+
+
+
+if __name__ == "__main__":
+    print('in main')
+    parser = argparse.ArgumentParser(description="This script handles downloading, processing and formatting of omics data files for the Bladder PDO project")
+    parser.add_argument('-s', '--samples', help='Path to sample file',default=None)
+    parser.add_argument('-g', '--genes', help='Path to genes file', default = None)
+    parser.add_argument('-c', '--copy', help='Flag to capture copy number data', action='store_true', default=False)
+    parser.add_argument('-m', '--mutation', help='Flag to capture mutation data', action='store_true', default=False)
+    parser.add_argument('-e', '--expression', help='Flag to capture transcriptomic data', action='store_true', default=False)
+    parser.add_argument('-i', '--geolink', help=".csv file that is the output of 'CNV-segfile-anotation.R")
+    parser.add_argument('-t', '--token', help='Synapse token')
+
+    args = parser.parse_args()
+    print("Logging into Synapse")
+    PAT = args.token
+    synObject = synapseclient.login(authToken=PAT)
+    print('gene file is:')
+    print(args.genes)
+    print('sample file is :')
+    print(args.samples)
+    genes = pd.read_csv(args.genes)
+    samples = pd.read_csv(args.samples)
+
+    if args.expression:
+        get_bladder_pdo_transcriptomics(args.geolink, samples, genes).to_csv("/tmp/bladderpdo_transcriptomics.csv", index=False)
+
+    if args.mutation:
+        get_bladder_pdo_mutations(synObject, samples, genes).to_csv('/tmp/bladderpdo_mutations.csv', index=False)
+    
+    if args.copy:
+        get_bladder_pdo_copynumber(synObject, samples, genes).to_csv("/tmp/bladderpdo_copynumber.csv", index=False)

build/bladderpdo/02_createBladderPDODrugsFile.py

@@ -0,0 +1,54 @@
+import synapseclient
+import pandas as pd
+import numpy as np
+import argparse
+import os
+# for testing locally
+#from utils.pubchem_retrieval import update_dataframe_and_write_tsv
+# for building in docker
+from pubchem_retrieval import update_dataframe_and_write_tsv
+
+
+def create_bladder_pdo_drugs_file(synObject, prevDrugFilepath, outputPath):
+    bladder_dir = synObject.get('syn64765430')
+    filenames = list(synObject.getChildren(parent='syn64765430', includeTypes=['file']))
+    bladder_drugs = pd.DataFrame({'drugNames' : [str]})
+    # '-4' - there are 4 nondrug files in this directory. 
+    for i in range(len(filenames)-4):
+        bladder_drugs.loc[i,'drugNames'] = filenames[i]['name'].split(")")[1].split("(")[0].split(".")[0].strip()
+
+    # get unique drugs 
+    newdrugnames = bladder_drugs['drugNames'].unique()
+    # use helper functions in pubchem_retrieval.py 
+    alldrugs = []
+    if prevDrugFilepath is not None and prevDrugFilepath is not "":
+        prevdrugs = [pd.read_csv(t,sep='\t') for t in prevDrugFilepath.split(',')]
+        alldrugs = pd.concat(prevdrugs).drop_duplicates()
+
+        imps = alldrugs[alldrugs.chem_name.isin(newdrugnames)]
+        newdrugs = alldrugs[alldrugs.improve_drug_id.isin(imps.improve_drug_id)]
+        
+        ##write drugs
+        newdrugs.to_csv(outputPath, sep='\t', index=False)
+
+    if len(alldrugs)==0 or len(newdrugnames)>len(set(newdrugs.improve_drug_id)): #we have more names we didn't match
+        print('Missing drugs in existing file, querying pubchem')
+        update_dataframe_and_write_tsv(newdrugnames,outputPath)
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(description="This script handles downloading, processing and formatting of drug data files for the Lee Bladder PDO project")
+    parser.add_argument('-d', '--prevDrugFilePath', help='Path to a previous drug file for sarcpdo', nargs="?", default = None)
+    parser.add_argument('-o', '--outputPath', help='Output path for updated sarcpdo drug file', default = "/tmp/sarcpdo_drugs.tsv") 
+    parser.add_argument('-t', '--token', help='Synapse token')
+
+    args = parser.parse_args()
+    print("Logging into Synapse")
+    PAT = args.token
+    synObject = synapseclient.login(authToken=PAT)
+    if args.prevDrugFilePath:
+        previousDrugs = args.prevDrugFilePath
+    else:
+        previousDrugs = None
+    create_bladder_pdo_drugs_file(synObject, previousDrugs, args.outputPath)

build/bladderpdo/03_createBladderPDOExperimentFile.py

@@ -0,0 +1,72 @@
+import synapseclient
+import pandas as pd
+import argparse
+
+
+def get_bladder_pdo_experiments(synObject, samples, drugs):
+    # get list of syn id info
+    files = list(synObject.getChildren(parent='syn64765430', includeTypes=['file']))
+    # load sample sheet and format _ to .
+    # load conversion table and remove trailing _T
+    conversion_table = synObject.get(files[50]['id'])
+    conversion_table_df = pd.read_excel(conversion_table.path, header=None)
+    conversion_table_df[2] = conversion_table_df[1].str.rsplit("_", expand=True)[0]#.replace(".", "_")
+    conversion_table_df[3] = conversion_table_df[2].str.replace(".", "_")
+    #print(conversion_table_df.head)
+
+    # initiate empty pd.dat.frame
+    drug_df = pd.DataFrame()
+    # for each drug, 
+    for i in range(len(files)-4):
+        drug_table_syn =synObject.get(files[i]['id'])
+        drug_table = pd.read_csv(drug_table_syn.path, sep="\t")
+    # melt
+    # link to conversion table
+    # link to sample sheet
+    # Rename, add columns
+        melted_single_drug = drug_table.melt(id_vars = 'Unnamed: 0', value_vars = drug_table.columns[1:], var_name="sample")
+        melted_single_drug['linkID'] = melted_single_drug['sample'].str.split(".", expand=True)[0] 
+        drugdata = melted_single_drug.merge(conversion_table_df, left_on = 'linkID', right_on = 0, how='left')[['Unnamed: 0', 'value', 3]]
+        #print(drugdata.head)
+        drugdata_with_improvesample = drugdata.merge(samples, left_on = 3, right_on='common_name')
+    
+    #   print(drugdata_with_improvesample.head)
+        drugdata_with_improvesample = drugdata_with_improvesample[['Unnamed: 0', 'value', 'improve_sample_id']]
+        #print(drugdata_with_improvesample.columns)
+        drugdata_with_improvesample = drugdata_with_improvesample.rename({"Unnamed: 0" : "DOSE", 'value' : 'GROWTH'}, axis=1)
+        #print(drugdata_with_improvesample.columns)
+
+        selected_drugdata = drugdata_with_improvesample
+        selected_drugdata['chem_name'] = files[i]['name'].split(")")[1].split("(")[0].split(".")[0].strip().lower()
+        #print(selected_drugdata.head)
+        drugdata_with_both_improveIds = selected_drugdata.merge(drugs[['improve_drug_id', 'chem_name']], how='left')
+        final_drugdata = drugdata_with_both_improveIds[['DOSE', 'GROWTH', 'improve_sample_id', 'improve_drug_id']]
+        final_drugdata = final_drugdata.rename({'improve_drug_id' : "Drug"}, axis=1)
+        final_drugdata['study'] = "Lee etal 2018 Bladder PDOs"
+        final_drugdata['source'] = "Synapse"
+        final_drugdata['time'] = 6
+        final_drugdata['time_unit'] = 'days'
+        #print(final_drugdata.head)
+        # append to dataframe
+        drug_df = pd.concat([drug_df, final_drugdata])
+    
+    return drug_df
+
+
+if __name__ == "__main__": 
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-t', '--token', help='Synapse authentication token')
+    parser.add_argument('-s', '--curSampleFile', help='Sample mapping file for bladder pdo samples')
+    parser.add_argument('-d', '--drugfile', help='Drug mapping file for bladder pdo samples')
+    parser.add_argument('-o', '--output', default = '/tmp/bladderpdo_doserep.tsv',help='Output file to be read into curve fitting code')
+
+    args = parser.parse_args()
+    print("Logging into Synapse")
+    PAT = args.token
+    synObject = synapseclient.login(authToken=PAT)
+    drug_df = pd.read_csv(args.drugfile, sep='\t')
+    samples_df = pd.read_csv(args.curSampleFile)
+
+    doseresponse_data = get_bladder_pdo_experiments(synObject, samples_df, drug_df)
+    doseresponse_data.to_csv(args.output, sep='\t')
+

build/bladderpdo/CNV-segfile-annotation.R

@@ -0,0 +1,59 @@
+#!/usr/bin/env Rscript
+args = commandArgs(trailingOnly=TRUE)
+
+# If BiocManager not installed, install it
+if (!require("BiocManager", quietly = TRUE))
+  install.packages("BiocManager")
+BiocManager::install(version = "3.20")
+
+library(GenomicRanges)
+library(Homo.sapiens)
+
+
+# function to return gene mappings for each coordinate range (row) in the segfile
+splitColumnByOverlap <-
+  function(query, subject, column="ENTREZID", ...)
+  {
+    olaps <- findOverlaps(query, subject, ...)
+    f1 <- factor(subjectHits(olaps),
+                 levels=seq_len(subjectLength(olaps)))
+    splitAsList(mcols(query)[[column]][queryHits(olaps)], f1)
+  }
+
+segfile = read.csv(args[1])
+
+# create genomic ranges object from segfile for use with findOverlaps()
+gr <- GRanges(seqnames = Rle(paste0('chr', segfile$chrom)), 
+              ranges = IRanges(segfile$loc.start, end = segfile$loc.end), 
+              strand = Rle(c("-", "0", "+")[sign(segfile$loc.start) +2]), 
+              score = segfile$seg.mean, 
+              ID = segfile$ID)
+
+# create genes GRanges obj from database
+genes<- genes(Homo.sapiens, columns =c("ENTREZID"))
+
+geneHitsByRow <- splitColumnByOverlap(genes, gr, "ENTREZID")
+
+# create matrices of annotations and scores/patient IDs in segfile
+# with a catch if there are no gene hits found
+matrixresults <- list()
+noresults <- c()
+for (i in 1:length(geneHitsByRow)){
+  if(length(geneHitsByRow[[i]])>0){
+    
+    matrixresults[[i]] <- data.frame(ENTREZID = as.matrix(geneHitsByRow[[i]]), rep(mcols(gr[i,]), length(geneHitsByRow[[i]])))
+  }else{
+    noresults <- c(noresults, i)
+  matrixresults[[i]] <- data.frame(ENTREZID = NA, mcols(gr)[i,])
+}
+}
+# concatenate annotated matrices
+allCNV <- do.call(rbind, matrixresults)
+
+# drop NAs
+completeAllCNV <- allCNV[complete.cases(allCNV),]
+# aggregate scores for the same genes (that came from different regions) for the same patient ID 
+aggregatedAllCNV <- aggregate(completeAllCNV$score, by = list(ENTREZID =completeAllCNV$ENTREZID, ID = completeAllCNV$ID), FUN=mean)
+names(aggregatedAllCNV)[names(aggregatedAllCNV)=='x'] <- "score"
+# write results to csv for further processing in Python
+write.csv(aggregatedAllCNV, args[2], row.names=F, quote =F)