Mercurial > repos > bimib > marea

changeset 48:e4235b5231e4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bimib
date Sun, 23 Feb 2020 09:41:14 -0500
parents 3af9d394367c
children 2c2a11aa1e02
files Marea/marea.py Marea/marea.xml Marea/ras_generator.xml
diffstat 3 files changed, 172 insertions(+), 95 deletions(-) [+]
line wrap: on
line diff
--- a/Marea/marea.py	Wed Feb 19 10:44:52 2020 -0500
+++ b/Marea/marea.py	Sun Feb 23 09:41:14 2020 -0500
@@ -18,14 +18,17 @@
     parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
                                      description = 'process some value\'s'+
                                      ' genes to create a comparison\'s map.')
-    parser.add_argument('-rs', '--rules_selector', 
+    parser.add_argument('-cr', '--custom_rules', 
                         type = str,
-                        default = 'HMRcore',
-                        choices = ['HMRcore', 'Recon', 'Custom'], 
-                        help = 'chose which type of dataset you want use')
-    parser.add_argument('-cr', '--custom',
+                        default = 'false',
+                        choices = ['true', 'false'],
+                        help = 'choose whether to use custom rules')
+    parser.add_argument('-cc', '--custom_rule',
                         type = str,
-                        help='your dataset if you want custom rules')
+                        help='custom rules to use')
+    parser.add_argument('-cm', '--custom_map',
+                        type = str,
+                        help='custom map to use')
     parser.add_argument('-n', '--none',
                         type = str,
                         default = 'true',
@@ -55,13 +58,6 @@
     parser.add_argument('-ic', '--input_class', 
                         type = str, 
                         help = 'sample group specification')
-    parser.add_argument('-cm', '--custom_map', 
-                        type = str, 
-                        help = 'custom map')
-    parser.add_argument('-yn', '--yes_no', 
-                        type = str,
-                        choices = ['yes', 'no'],
-                        help = 'if make or not custom map')
     parser.add_argument('-gs', '--generate_svg',
                         type = str,
                         default = 'true',
@@ -551,7 +547,7 @@
         react = recon.reactions
         rules = [react[i].gene_reaction_rule for i in range(len(react))]
         ids = [react[i].id for i in range(len(react))]
-    except cb.io.sbml3.CobraSBMLError:
+    except cb.io.sbml.CobraSBMLError:
         try:
             data = (pd.read_csv(data, sep = '\t', dtype = str, engine='python')).fillna('')
             if len(data.columns) < 2:
@@ -693,8 +689,8 @@
             tmp_csv.to_csv(tab, sep = '\t', index = False, header = header)
             
             if create_svg or create_pdf:
-                if args.rules_selector == 'HMRcore' or (args.rules_selector == 'Custom'
-                                                        and args.yes_no == 'yes'):
+                if args.custom_rules == 'false' or (args.custom_rules == 'true'
+                                                        and args.custom_map != ''):
                     fix_map(tmp, core_map, threshold_P_V, threshold_F_C, max_F_C)
                     file_svg = 'result/' + i + '_vs_' + j + ' (SVG Map).svg'
                     with open(file_svg, 'wb') as new_map:
@@ -760,13 +756,6 @@
                         
     elif comparison == "onevsmany":
         for i, j in it.combinations(class_pat.keys(), 2):
-
-            if i != control and j != control:
-                print(str(control) + " " + str(i) + " " + str(j))
-                #Se è un confronto fra elementi diversi dal controllo, skippo
-                continue
-            
-            print('vado')
             tmp = {}
             count = 0
             max_F_C = 0
@@ -823,13 +812,9 @@
 
     if os.path.isdir('result') == False:
         os.makedirs('result')
-    
-    if args.rules_selector == 'HMRcore':        
-        recon = pk.load(open(args.tool_dir + '/local/HMRcore_rules.p', 'rb'))
-    elif args.rules_selector == 'Recon':
-        recon = pk.load(open(args.tool_dir + '/local/Recon_rules.p', 'rb'))
-    elif args.rules_selector == 'Custom':
-        ids, rules, gene_in_rule = make_recon(args.custom)
+        
+    if args.custom_rules == 'true':
+        ids, rules, gene_in_rule = make_recon(args.custom_rule)
 
     class_pat = {}
     
@@ -885,15 +870,13 @@
         if resolve_rules_float != None:
             class_pat = split_class(classes, resolve_rules_float)
     	
-    if args.rules_selector == 'Custom':
-        if args.yes_no == 'yes':
-            try:
-                core_map = ET.parse(args.custom_map)
-            except (ET.XMLSyntaxError, ET.XMLSchemaParseError):
-                sys.exit('Execution aborted: custom map in wrong format')
-        elif args.yes_no == 'no':
-            core_map = ET.parse(args.tool_dir + '/local/HMRcoreMap.svg')
-    else:       
+       
+    if args.custom_rules == 'true':
+        try:
+            core_map = ET.parse(args.custom_map)
+        except (ET.XMLSyntaxError, ET.XMLSchemaParseError):
+            sys.exit('Execution aborted: custom map in wrong format')
+    else:
         core_map = ET.parse(args.tool_dir+'/local/HMRcoreMap.svg')
         
     maps(core_map, class_pat, ids, args.pValue, args.fChange, create_svg, create_pdf, args.comparison, args.control)
--- a/Marea/marea.xml	Wed Feb 19 10:44:52 2020 -0500
+++ b/Marea/marea.xml	Sun Feb 23 09:41:14 2020 -0500
@@ -1,4 +1,4 @@
-<tool id="MaREA" name="Metabolic Reaction Enrichment Analysis" version="1.0.7">
+<tool id="MaREA" name="Metabolic Reaction Enrichment Analysis" version="1.0.8">
 	<macros>
 		<import>marea_macros.xml</import>
 	</macros>
@@ -28,6 +28,11 @@
                 ${data.input_name}
             #end for
             --comparison ${cond.comparis.comparison}
+            #if $cond.advanced.cond_map == 'true':
+            	--custom_rules true
+            	--custom_rule ${cond.advanced.cond_map.custom_rule}
+            	--custom_map ${cond.advanced.cond_map.custom_map}
+            #end if
             #if $cond.advanced.choice == 'true':
       	       --pValue ${cond.advanced.pValue}
       	       --fChange ${cond.advanced.fChange}
@@ -47,6 +52,11 @@
             #if $cond.comparis.comparison == 'onevsmany'
             	--control ${cond.comparis.controlgroup}
             #end if
+            #if $cond.advanced.cond_map == 'true':
+            	--custom_rules true
+            	--custom_rule ${cond.advanced.cond_map.custom_rule}
+            	--custom_map ${cond.advanced.cond_map.custom_map}
+            #end if
             #if $cond.advanced.choice == 'true':
       	        --pValue ${cond.advanced.pValue}
       	        --fChange ${cond.advanced.fChange}
@@ -90,11 +100,12 @@
 					<when value="false"></when>
 					<when value="true">
 						<conditional name="cond_map">
-						<param name="choice" type="boolean" checked="false" label="Use custom map?" help="Use this option only if you have generated RAS using a custom set of rules">
+						<param name="choice" type="boolean" checked="false" label="Use custom map and rules?" help="Use this option only if you have generated RAS using a custom set of rules">
 							<option value="false" selected="true">No</option>
 							<option value="true">Yes</option>
 						</param>
-						<when value="true">							
+						<when value="true">	
+							<param name="Custom_rule" argument="--custom_rule" type="data" format="tabular, csv, tsv, xml" label="Custom rules" />					
 							<param name="Custom_map" argument="--custom_map" type="data" format="xml, svg" label="custom-map.svg"/>
 						</when>
 					</conditional>
@@ -131,6 +142,7 @@
 							<option value="true">Yes</option>
 						</param>
 						<when value="true">							
+							<param name="Custom_rule" argument="--custom_rule" type="data" format="tabular, csv, tsv, xml" label="Custom rules" />					
 							<param name="Custom_map" argument="--custom_map" type="data" format="xml, svg" label="custom-map.svg"/>
 						</when>
 					</conditional>
@@ -155,71 +167,108 @@
 What it does
 -------------
 
-This tool analyzes RNA-seq dataset(s) as described in Graudenzi et al."`MaREA`_: Metabolic feature extraction, enrichment and visualization of RNAseq data" bioRxiv (2018): 248724.
+This tool analyzes and visualizes differences in the Reaction Activity Scores (RASs) of groups of samples, as computed by the Expression2RAS tool, of groups of samples.
 
 Accepted files are: 
-    - option 1) two or more RNA-seq datasets, each referring to samples in a given condition/class. The user can specify a label for each class (as e.g. "*classA*" and "*classB*");
-    - option 2) one RNA dataset and one class-file specifying the class/condition each sample belongs to.
+    - option 1) two or more RAS datasets, each referring to samples in a given group. The user can specify a label for each group (as e.g. "classA" and "classB");
+    - option 2) one RAS dataset and one group-file specifying the group each sample belongs to.
+    
+RAS datasets format: tab-separated text files, reporting the RAS value of each reaction (row) for a given sample (column).
+
+Column header: sample ID.
+Raw header: reaction ID. 
 
 Optional files:
-    - custom GPR (Gene-Protein-Reaction) rules. Two accepted formats:
-
-	* (Cobra Toolbox and CobraPy compliant) xml of metabolic model;
-	* .csv file specifyig for each reaction ID (column 1) the corresponding GPR rule (column 2).
     - custom svg map. Graphical elements must have the same IDs of reactions. See HmrCore svg map for an example.
 
 The tool generates:
-    1) a tab-separated file: reporting fold-change and p-values of reaction activity scores (RASs) between a pair of conditions/classes;
-    2) a metabolic map file (downlodable as .svg): visualizing up- and down-regulated reactions between a pair of conditions/classes;
-    3) a log file (.txt).
+    - 1) a tab-separated file: reporting fold-change and p-values of reaction activity scores (RASs) between a pair of conditions/classes;
+    - 2) a metabolic map file (downloadable as .svg): visualizing up- and down-regulated reactions between a pair of conditions/classes;
+    - 3) a log file (.txt).
+    
+Output options:
+To calculate P-Values and Fold-Changes and to enrich maps, comparisons are performed for each possible pair of groups (default option ‘One vs One’).
 
-RNA-seq datasets format: tab-separated text files, reporting the expression level (e.g., TPM, RPKM, ...) of each gene (row) for a given sample (column). Header: sample ID.
-
-Class-file format: each row of the class-file reports the sample ID (column1) and the label of the class/condition the sample belongs to (column 2).
-
-To calculate P-Values and Fold-Changes and to generate maps, comparisons are performed for each possible pair of classes.
+Alternative options are:
+    - comparison of each group vs. the rest of samples (option ‘One vs Rest’)
+    - comparison of each group vs. a control group (option ‘One vs Control). If this option is selected the user must indicate the control group label.
 
 Output files will be named as classA_vs_classB. Reactions will conventionally be reported as up-regulated (down-regulated) if they are significantly more (less) active in class having label "classA".
 
-
 Example input
 -------------
 
-**"Custom Rules"** option:
+"RAS of group 1 + RAS of group 2 + ... + RAS of group N" option:
 
-Custom Rules Dastaset:
+RAS Dataset 1:
 
-@CUSTOM_RULES_EXEMPLE@
-
-**"RNAseq of group 1 + RNAseq of group 2 + ... + RNAseq of group N"** option:
++------------+----------------+----------------+----------------+ 
+| Reaction ID|   TCGAA62670   |   TCGAA62671   |   TCGAA62672   |  
++============+================+================+================+
+| r1642      |    0.523167    |    0.371355    |    0.925661    |  
++------------+----------------+----------------+----------------+    
+| r1643      |    0.568765    |    0.765567    |    0.456789    |    
++------------+----------------+----------------+----------------+    
+| r1640      |    0.876545    |    0.768933    |    0.987654    |  
++------------+----------------+----------------+----------------+
+| r1641      |    0.456788    |    0.876543    |    0.876542    |    
++------------+----------------+----------------+----------------+    
+| r1646      |    0.876543    |    0.786543    |    0.897654    |   
++------------+----------------+----------------+----------------+
 
-RNA-seq Dataset 1:						
-
-@DATASET_EXEMPLE1@
+RAS Dataset 2:
 
-RNA-seq Dataset 2:
++------------+----------------+----------------+----------------+ 
+| Reaction ID|   TCGAA62670   |   TCGAA62671   |   TCGAA62672   |  
++============+================+================+================+
+| r1642      |    0.523167    |    0.371355    |    0.925661    |  
++------------+----------------+----------------+----------------+    
+| r1643      |    0.568765    |    0.765567    |    0.456789    |    
++------------+----------------+----------------+----------------+    
+| r1640      |    0.876545    |    0.768933    |    0.987654    |  
++------------+----------------+----------------+----------------+
+| r1641      |    0.456788    |    0.876543    |    0.876542    |    
++------------+----------------+----------------+----------------+    
+| r1646      |    0.876543    |    0.786543    |    0.897654    |   
++------------+----------------+----------------+----------------+
 
-@DATASET_EXEMPLE2@
-
-**"RNAseq of all samples + sample group specification"** option:
+"RAS of all samples + sample group specification" option:
 
-RNA-seq Dataset:
+RAS Dataset:
 
-@DATASET_EXEMPLE1@
++------------+----------------+----------------+----------------+ 
+| Reaction ID|   TCGAA62670   |   TCGAA62671   |   TCGAA62672   |  
++============+================+================+================+
+| r1642      |    0.523167    |    0.371355    |    0.925661    |  
++------------+----------------+----------------+----------------+    
+| r1643      |    0.568765    |    0.765567    |    0.456789    |    
++------------+----------------+----------------+----------------+    
+| r1640      |    0.876545    |    0.768933    |    0.987654    |  
++------------+----------------+----------------+----------------+
+| r1641      |    0.456788    |    0.876543    |    0.876542    |    
++------------+----------------+----------------+----------------+    
+| r1646      |    0.876543    |    0.786543    |    0.897654    |   
++------------+----------------+----------------+----------------+
 
-Class-file:
+Group-file
 
-+------------+------------+   
-| Patient_ID |    class   |   
-+============+============+   
-| TCGAAA3529 |     MSI    |   
-+------------+------------+    
-| TCGAA62671 |     MSS    |    
-+------------+------------+    
-| TCGAA62672 |     MSI    |   
-+------------+------------+
++---------------+-----------+
+| Patient ID    |   Class   | 
++===============+===========+
+| TCGAAA3529    |    MSI    | 
++---------------+-----------+  
+| TCGAA62671    |    MSS    |    
++---------------+-----------+   
+| TCGAA62672    |    MSI    |
++---------------+-----------+
 
-|
+Advanced options
+----------------
+
+P-Value threshold: the threshold used for significance Kolmogorov-Smirnov (KS) test, to verify whether the distributions of RASs over the samples in two sets are significantly different
+
+Fold-Change threshold: threshold of the fold-change between the average RAS of two groups. Among the reactions that pass the KS test, only fold-change values larger than the indicated threshold will be visualized on the output metabolic map;
+
 
 .. class:: infomark
 
--- a/Marea/ras_generator.xml	Wed Feb 19 10:44:52 2020 -0500
+++ b/Marea/ras_generator.xml	Sun Feb 23 09:41:14 2020 -0500
@@ -1,4 +1,4 @@
-<tool id="MaREA RAS Generator" name="Expression2RAS" version="1.0.2">
+<tool id="MaREA RAS Generator" name="Expression2RAS" version="1.0.3">
     <description>- Reaction Activity Scores computation</description>
     <macros>
         <import>marea_macros.xml</import>
@@ -31,17 +31,6 @@
             </when>
             <when value="Custom">
                 <param name="Custom_rules" type="data" format="tabular, csv, tsv, xml" label="Custom rules" />
-                <conditional name="cond_map">
-                    <param name="yes_no" type="select" label="Custom map? (optional)">
-                        <option value="no" selected="true">no</option>
-                        <option value="yes">yes</option>
-                    </param>
-                    <when value="yes">
-                        <param name="Custom_map" argument="--custom_map" type="data" format="xml, svg" label="custom-map.svg"/>
-                    </when>
-                    <when value="no">
-                    </when>
-                </conditional>
             </when>
         </conditional>
         <param name="input" argument="--input" type="data" format="tabular, csv, tsv" label="Gene Expression dataset:" />
@@ -59,6 +48,62 @@
 
 What it does
 -------------
+
+This tool computes Reaction Activity Scores from gene expression (RNA-seq) dataset(s), as described in Graudenzi et al. Integration of transcriptomic data and metabolic networks in cancer samples reveals highly significant prognostic power. Journal of Biomedical Informatics, 2018, 87: 37-49.
+ 
+Accepted files:
+    - A gene expression dataset
+ 
+Format:
+Tab-separated text file reporting the normalized expression level (e.g., TPM, RPKM, ...) of each gene (row) for a given sample (column).
+Column header: sample ID.
+Row header: gene ID.
+ 
+ 
+Optional files:
+    - custom GPR (Gene-Protein-Reaction) rules. Two accepted formats:
+
+	* (Cobra Toolbox and CobraPy compliant) xml of metabolic model;
+	* .csv file specifyig for each reaction ID (column 1) the corresponding GPR rule (column 2).
+ 
+Computation option ‘(A and NaN) solved as (A)’:
+In case of missing expression value, referred to as NaN (Not a Number), for a gene joined with an AND operator in a given GPR rule, the rule ‘A and NaN’
+ 
+If YES is selected: the GPR will be solved as A.
+ 
+If NO is selected: the GPR will be disregarded tout-court (i.e., treated as NaN).
+
+Example input
+-------------
+
+Custom GPR rules:
+
++------------+--------------------------------------+   
+| id         |         rule (with entrez-id         |   
++============+======================================+   
+| r1642      |             155060 or 10357          |   
++------------+--------------------------------------+    
+| r1643      |        155060 or 100134869           |    
++------------+--------------------------------------+    
+| r1640      |     155060 and 100134869 or 10357    |   
++------------+--------------------------------------+
+
+RNA-seq dataset:
+
++------------+----------------+----------------+----------------+ 
+| Hugo_ID    |   TCGAA62670   |   TCGAA62671   |   TCGAA62672   |  
++============+================+================+================+
+| HGNC:24086 |    0.523167    |    0.371355    |    0.925661    |  
++------------+----------------+----------------+----------------+    
+| HGNC:24086 |    0.568765    |    0.765567    |    0.456789    |    
++------------+----------------+----------------+----------------+    
+| HGNC:9876  |    0.876545    |    0.768933    |    0.987654    |  
++------------+----------------+----------------+----------------+
+| HGNC:9     |    0.456788    |    0.876543    |    0.876542    |    
++------------+----------------+----------------+----------------+    
+| HGNC:23    |    0.876543    |    0.786543    |    0.897654    |   
++------------+----------------+----------------+----------------+
+
 ]]>
     </help>
     <expand macro="citations" />