Merge pull request #347 from PNNL-CompBio/bladder_pdo

Bladder PDO dataset addition

Author: jjacobson95

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
+        dose_resp_df = pd.concat([drug_df, final_drugdata])
+    
+    return dose_resp_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)

build/bladderpdo/README.md

@@ -0,0 +1,14 @@
+#!/bin/bash
+set -euo pipefail
+
+trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit 1' ERR
+
+echo "Running script with token and drugFile $1"
+# for running locally (from build directory):
+#python3 -m bladderpdo.02_createBladderPDODrugsFile --token $SYNAPSE_AUTH_TOKEN -d $1 -o ./bladderpdo/bladderpdo_drugs.tsv
+python3 02_createBladderPDODrugsFile.py --token $SYNAPSE_AUTH_TOKEN -d $1 -o /tmp/bladderpdo_drugs.tsv
+
+echo "Running build_drug_desc.py..."
+#for running locally: 
+#python3 utils/build_drug_desc.py --drugtable ./bladderpdo/bladderpdo_drugs.tsv --desctable ./bladderpdo/bladderpdo_drug_descriptors.tsv.gz
+python3 build_drug_desc.py --drugtable /tmp/bladderpdo_drugs.tsv --desctable /tmp/bladderpdo_drug_descriptors.tsv.gz