Merge pull request #54 from AlexandrovLab/u15

U15

Author: marcos-diazg

README.md

@@ -127,14 +127,15 @@ Analyze.cosmic_fit( samples,
 | **samples** | String | Path to input file for `input_type`:<ul><li>"matrix"</li><li>"seg:TYPE"</li></ul> Path to input folder for `input_type`:<ul><li>"vcf"</li></ul>|
 | **output** | String | Path to the output folder. |
  | **input_type** | String | The type of input:<br><ul><li>"matrix": used for table format inputs using a tab-separated file where the rows are mutation types and the columns are sample IDs.</li><li>"vcf": used for mutation calling file inputs (VCFs, MAFs or simple text files).</li><li>"seg:TYPE": used for a multi-sample segmentation file for copy number analysis. The accepted callers for TYPE are the following {"ASCAT", "ASCAT_NGS", "SEQUENZA", "ABSOLUTE", "BATTENBERG", "FACETS", "PURPLE", "TCGA"}. For example, when using segmentation file from BATTENBERG then set input_type to "seg:BATTENBERG".</li></ul> The default value is "matrix".|
-| **context_type**| String| Required context type if `input_type` is "vcf". `context_type` takes which context type of the input data is considered for assignment. Valid options include "96", "288", "1536", "DINUC", and "INDEL". The default value is "96".|
+| **context_type**| String| Required context type if `input_type` is "vcf". `context_type` takes which context type of the input data is considered for assignment. Valid options include "96", "288", "1536", "DINUC", and "ID". The default value is "96".|
 | **signatures** | String | Path to a tab delimited file that contains the signature table where the rows are mutation types and colunms are signature IDs. |
 | **genome_build** | String | The reference genome build. List of supported genomes: "GRCh37", "GRCh38", "mm9", "mm10" and "rn6". The default value is "GRCh37". If the selected genome is not in the supported list, the default genome will be used. |
 | **cosmic_version** | Float | Takes a positive float among 1, 2, 3, 3.1, 3.2 and 3.3. Defines the version of the COSMIC reference signatures. The default value is 3.3. |
 | **new_signature_thresh_hold**| Float | Parameter in cosine similarity to declare a new signature. Applicable for decompose_fit only. The default value is 0.8. |
 | **exclude_signature_subgroups** | List | Removes the signatures corresponding to specific subtypes for better fitting. The usage is given above. The default value is None. |
 | **exome** | Boolean | Defines if the exome renormalized signatures will be used. The default value is False. |
-| **export_probabilities** | Boolean | Defines if the probability matrix is created. The default value is True. |
+| **export_probabilities** | Boolean | Defines if the probability matrix per mutational context for all samples is created. The default value is True. |
+| **export_probabilities_per_mutation** | Boolean | Defines if the probability matrices per mutation for all samples are created. Only available when `input_type` is "vcf". The default value is False. |
 | **make_plots** | Boolean | Toggle on and off for making and saving all plots. The default value is True. |
 | **verbose** | Boolean | Prints statements. The default value is False. |
 

SigProfilerAssignment/Analyzer.py

@@ -1,11 +1,11 @@
 from SigProfilerAssignment import decomposition as decomp
 
-def decompose_fit(samples, output, signatures=None, signature_database=None,nnls_add_penalty=0.05, nnls_remove_penalty=0.01, initial_remove_penalty=0.05,genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,new_signature_thresh_hold=0.8,exclude_signature_subgroups=None,exome=False,input_type='matrix',context_type="96",export_probabilities=True):
+def decompose_fit(samples, output, signatures=None, signature_database=None,nnls_add_penalty=0.05, nnls_remove_penalty=0.01, initial_remove_penalty=0.05,genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,new_signature_thresh_hold=0.8,exclude_signature_subgroups=None,exome=False,input_type='matrix',context_type="96",export_probabilities=True, export_probabilities_per_mutation=False):
 
-                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= True,denovo_refit_option=False,cosmic_fit_option=False,devopts=devopts,new_signature_thresh_hold=new_signature_thresh_hold,exclude_signature_subgroups=exclude_signature_subgroups,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities)
+                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= True,denovo_refit_option=False,cosmic_fit_option=False,devopts=devopts,new_signature_thresh_hold=new_signature_thresh_hold,exclude_signature_subgroups=exclude_signature_subgroups,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities, export_probabilities_per_mutation=export_probabilities_per_mutation)
 
-def denovo_fit(samples,  output, signatures=None, signature_database=None,nnls_add_penalty=0.05,nnls_remove_penalty=0.01, initial_remove_penalty=0.05, genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,new_signature_thresh_hold=0.8,exome=False,input_type='matrix',context_type="96",export_probabilities=True):
-                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, new_signature_thresh_hold=new_signature_thresh_hold, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= False,denovo_refit_option=True,cosmic_fit_option=False,devopts=devopts,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities)
+def denovo_fit(samples,  output, signatures=None, signature_database=None,nnls_add_penalty=0.05,nnls_remove_penalty=0.01, initial_remove_penalty=0.05, genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,new_signature_thresh_hold=0.8,exome=False,input_type='matrix',context_type="96",export_probabilities=True, export_probabilities_per_mutation=False):
+                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, new_signature_thresh_hold=new_signature_thresh_hold, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= False,denovo_refit_option=True,cosmic_fit_option=False,devopts=devopts,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities, export_probabilities_per_mutation=export_probabilities_per_mutation)
 
-def cosmic_fit(samples, output, signatures=None, signature_database=None,nnls_add_penalty=0.05, nnls_remove_penalty=0.01, initial_remove_penalty=0.05,genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,exclude_signature_subgroups=None,exome=False,input_type='matrix',context_type="96",export_probabilities=True):
-                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= False,denovo_refit_option=False,cosmic_fit_option=True,devopts=devopts,exclude_signature_subgroups=exclude_signature_subgroups,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities)
+def cosmic_fit(samples, output, signatures=None, signature_database=None,nnls_add_penalty=0.05, nnls_remove_penalty=0.01, initial_remove_penalty=0.05,genome_build="GRCh37", cosmic_version=3.3, make_plots=True, collapse_to_SBS96=True,connected_sigs=True, verbose=False,devopts=None,exclude_signature_subgroups=None,exome=False,input_type='matrix',context_type="96",export_probabilities=True, export_probabilities_per_mutation=False):
+                decomp.spa_analyze(samples=samples,  output=output, signatures=signatures, signature_database=signature_database,nnls_add_penalty=nnls_add_penalty, nnls_remove_penalty=nnls_remove_penalty, initial_remove_penalty=initial_remove_penalty,genome_build=genome_build, cosmic_version=cosmic_version, make_plots=make_plots, collapse_to_SBS96=collapse_to_SBS96,connected_sigs=connected_sigs, verbose=verbose,decompose_fit_option= False,denovo_refit_option=False,cosmic_fit_option=True,devopts=devopts,exclude_signature_subgroups=exclude_signature_subgroups,exome=exome,input_type=input_type,context_type=context_type,export_probabilities=export_probabilities, export_probabilities_per_mutation=export_probabilities_per_mutation)

SigProfilerAssignment/decompose_subroutines.py

@@ -471,8 +471,8 @@ def signature_decomposition(signatures, mtype, directory, genome_build="GRCh37",
 #############################################################################################################
 def make_final_solution(processAvg, allgenomes, allsigids, layer_directory, m, index, allcolnames, process_std_error = "none", signature_stabilities = " ", \
                         signature_total_mutations= " ", signature_stats = "none",  cosmic_sigs=False, attribution= 0, denovo_exposureAvg  = "none", add_penalty=0.05, \
-                        remove_penalty=0.01, initial_remove_penalty=0.05, de_novo_fit_penalty=0.02, background_sigs=0, genome_build="GRCh37", sequence="genome", export_probabilities=True, \
-                        refit_denovo_signatures=True, collapse_to_SBS96=True, connected_sigs=True, pcawg_rule=False, verbose=False,make_plots = True):
+                        remove_penalty=0.01, initial_remove_penalty=0.05, de_novo_fit_penalty=0.02, background_sigs=0, genome_build="GRCh37", sequence="genome", export_probabilities=True, export_probabilities_per_mutation=False, \
+                        refit_denovo_signatures=True, collapse_to_SBS96=True, connected_sigs=True, pcawg_rule=False, verbose=False,make_plots = True, samples='./', input_type='matrix', denovo_refit_option=True):
 
     if processAvg.shape[0]==allgenomes.shape[0] and processAvg.shape[0] != 96:
         collapse_to_SBS96=False
@@ -774,15 +774,45 @@ def make_final_solution(processAvg, allgenomes, allsigids, layer_directory, m, i
         probability = probabilities(processAvg, exposureAvg, index, allsigids, allcolnames)
         probability=probability.set_index("Sample Names" )
 
-        if cosmic_sigs==False:
-            
+        if denovo_refit_option==True:
             if refit_denovo_signatures==True:
-                probability.to_csv(layer_directory+"/Activities"+"/"+"De_Novo_Mutation_Probabilities_refit.txt", "\t") 
+                probability.to_csv(layer_directory+"/Activities"+"/"+"De_Novo_MutationType_Probabilities_refit.txt", "\t") 
+            else:
+                probability.to_csv(layer_directory+"/Activities"+"/"+"De_Novo_MutationType_Probabilities.txt", "\t") 
+        if denovo_refit_option==False:
+            probability.to_csv(layer_directory+"/Activities"+"/"+"Decomposed_MutationType_Probabilities.txt", "\t") 
+    
+    if export_probabilities_per_mutation==True:
+        if export_probabilities==True:
+            if input_type=='vcf':
+                if m=='96' or m=='78' or m=='83':
+                    probability_per_mutation, samples_prob_per_mut = probabilities_per_mutation(probability, samples, m)
+
+                    if denovo_refit_option==True:
+                        if refit_denovo_signatures==True:
+                            ppm_file_name = "De_Novo_Mutation_Probabilities_refit"
+                            output_path_prob_per_mut = layer_directory+"/Activities"+"/"+ppm_file_name
+                        else:
+                            ppm_file_name = "De_Novo_Mutation_Probabilities"
+                            output_path_prob_per_mut = layer_directory+"/Activities"+"/"+ppm_file_name
+                    else:
+                        ppm_file_name = "Decomposed_Mutation_Probabilities"
+                        output_path_prob_per_mut = layer_directory+"/Activities"+"/"+ppm_file_name
+
+                    if not os.path.exists(output_path_prob_per_mut):
+                        os.makedirs(output_path_prob_per_mut)
+                    for matrix,sample in zip(probability_per_mutation, samples_prob_per_mut):
+                        matrix=matrix.set_index('Sample Names')
+                        matrix=matrix.sort_values(by=['Chr','Pos'])
+                        matrix.to_csv(layer_directory+"/Activities"+"/"+ ppm_file_name + "/" + ppm_file_name + "_" + sample + ".txt", "\t")
+                else:
+                    print('Probabilities per mutation are only calculated for SBS96, DBS78 and ID83 mutational contexts.')
             else:
-                probability.to_csv(layer_directory+"/Activities"+"/"+"De_Novo_Mutation_Probabilities.txt", "\t") 
-        if cosmic_sigs==True:
-            probability.to_csv(layer_directory+"/Activities"+"/"+"Decomposed_Mutation_Probabilities.txt", "\t") 
+                print('Probabilities per mutation are only calculated if input_type is "vcf".')
+        else:
+            print('Probabilities per mutation require to calculate probabilities per context type. Please re-run your analysis setting export_probabilites=True.')
 
+    # import pdb; pdb.set_trace()
 
     return exposures
 ################################################################### FUNCTION ONE ###################################################################
@@ -901,6 +931,59 @@ def probabilities(W, H, index, allsigids, allcolnames):
 
     return result
 
+
+################################################### Generation of probabilities for each processes given to A mutation ############################################
+def probabilities_per_mutation(probability_matrix, samples_path, m):  
+#
+    probability_matrix=probability_matrix.reset_index()
+#
+    if m=='96':
+        seqinfo_path = samples_path + '/output/vcf_files/SNV/'
+        interval_low = 3
+        interval_high = -1
+    if m=='78':
+        seqinfo_path = samples_path + '/output/vcf_files/DBS/'
+        interval_low = 4
+        interval_high = -2
+    if m=='83':
+        seqinfo_path = samples_path + '/output/vcf_files/ID/'
+        interval_low = 2
+        interval_high = 100
+#
+    seqinfo_files = os.listdir(seqinfo_path)
+    seqinfo_files.sort()
+#
+    all_mutations = pd.DataFrame()
+    for file in seqinfo_files:
+        try:
+            new = pd.read_csv(seqinfo_path + file, sep='\t',header=None)
+            all_mutations = pd.concat([all_mutations, new])
+        except (pd.errors.EmptyDataError):
+            pass
+    all_mutations[3] = all_mutations[3].str[interval_low:interval_high]
+    if m=='96' or m=='78':
+        del all_mutations[4]
+    else:
+        del all_mutations[6]
+        del all_mutations[5]
+        del all_mutations[4]
+
+    all_mutations.columns = ['Sample Names', 'Chr', 'Pos', 'MutationType']
+#
+    all_samples_mutations = [y for x, y in all_mutations.groupby('Sample Names')]
+#
+    prob_per_mut = []
+    sample_names = []
+    for sample_mutations in all_samples_mutations:
+        new = sample_mutations.merge(probability_matrix)
+        prob_per_mut.append(new)
+        sample_names.append(new['Sample Names'][0])
+#
+    result = [prob_per_mut, sample_names]
+#
+    return result
+
+
 def custom_signatures_plot(signatures, output):
     with PdfPages(output+'/Custom_Signature_Plots.pdf') as pdf:
         plt.figure(figsize=(10, 3))