Repository 'lifelines_km_cph_tool'
hg clone https://toolshed.g2.bx.psu.edu/repos/fubar/lifelines_km_cph_tool

Changeset 1:232b874046a7 (2023-08-10)
Previous changeset 0:dd49a7040643 (2023-08-09) Next changeset 2:dd5e65893cb8 (2023-08-10)
Commit message:
Uploaded
modified:
lifelines_tool/README.md
lifelines_tool/lifelineskmcph.xml
lifelines_tool/plotlykm.py
lifelines_tool/test-data/readme_sample
removed:
lifelines_tool/run_log.txt
b
diff -r dd49a7040643 -r 232b874046a7 lifelines_tool/README.md
--- a/lifelines_tool/README.md Wed Aug 09 11:12:16 2023 +0000
+++ b/lifelines_tool/README.md Thu Aug 10 07:15:22 2023 +0000
b
@@ -2,7 +2,7 @@
 
 ## Galaxy tool to run failure time models using lifelines
 
-## Install to your Galaxy server from the toolshed - search for lifelines_tool owned by fubar2
+## Install to your Galaxy server from the toolshed - search for lifelines_km_cph_tool owned by fubar2
 
 ### More at https://lazarus.name/demo/
 
@@ -25,7 +25,7 @@
 Should work with any tabular data with the required columns - time and status for observations.
 
 Issues to https://github.com/fubar2/lifelines_tool please.
-Autogenerated so pull requests are possibly meaningless but regeneration of a new version is easy so please tell me what is needed.
+Autogenerated so pull requests are possibly meaningless but regeneration of a new version should work.
 
 ## Tool made with the Galaxy ToolFactory: https://github.com/fubar2/galaxy_tf_overlay
 The current release includes this and a generic tabular version, and a java .jar wrapper in a history where the generating
b
diff -r dd49a7040643 -r 232b874046a7 lifelines_tool/lifelineskmcph.xml
--- a/lifelines_tool/lifelineskmcph.xml Wed Aug 09 11:12:16 2023 +0000
+++ b/lifelines_tool/lifelineskmcph.xml Thu Aug 10 07:15:22 2023 +0000
[
@@ -1,6 +1,6 @@
 <tool name="lifelineskmcph" id="lifelineskmcph" version="0.01">
   <!--Source in git at: https://github.com/fubar2/galaxy_tf_overlay-->
-  <!--Created by toolfactory@galaxy.org at 09/08/2023 17:43:16 using the Galaxy Tool Factory.-->
+  <!--Created by toolfactory@galaxy.org at 10/08/2023 15:48:43 using the Galaxy Tool Factory.-->
   <description>Lifelines KM and optional Cox PH models</description>
   <requirements>
     <requirement version="1.5.3" type="package">pandas</requirement>
@@ -40,8 +40,9 @@
 # km models for https://github.com/galaxyproject/tools-iuc/issues/5393
 # test as
 # python plotlykm.py --input_tab rossi.tab --htmlout "testfoo" --time "week" --status "arrest" --title "test" --image_dir images --cphcol="prio,age,race,paro,mar,fin"
+# Ross Lazarus July 2023
+import argparse
 
-import argparse
 import os
 import sys
 
@@ -51,15 +52,22 @@
 
 import pandas as pd
 
-# Ross Lazarus July 2023
 
+def trimlegend(v):
+    """
+    for int64 quintiles - must be ints - otherwise get silly legends with long float values
+    """
+    for i, av in enumerate(v):
+        x = int(av)
+        v[i] = str(x)
+    return v
 
 kmf = lifelines.KaplanMeierFitter()
 cph = lifelines.CoxPHFitter()
 
 parser = argparse.ArgumentParser()
 a = parser.add_argument
-a('--input_tab', default='', required=True)
+a('--input_tab', default='rossi.tab', required=True)
 a('--header', default='')
 a('--htmlout', default="test_run.html")
 a('--group', default='')
@@ -75,6 +83,7 @@
 df = pd.read_csv(args.input_tab, sep='\t')
 NCOLS = df.columns.size
 NROWS = len(df.index)
+QVALS = [.2, .4, .6, .8] # for partial cox ph plots
 defaultcols = ['col%d' % (x+1) for x in range(NCOLS)]
 testcols = df.columns
 if len(args.header.strip()) > 0:
@@ -106,14 +115,11 @@
     names = []
     times = []
     events = []
-    rmst = []
     for name, grouped_df in df.groupby(args.group):
         T = grouped_df[args.time]
         E = grouped_df[args.status]
         gfit = kmf.fit(T, E, label=name)
         kmf.plot_survival_function(ax=ax)
-        rst = lifelines.utils.restricted_mean_survival_time(gfit)
-        rmst.append(rst)
         names.append(str(name))
         times.append(T)
         events.append(E)
@@ -124,30 +130,50 @@
         results = lifelines.statistics.logrank_test(times[0], times[1], events[0], events[1], alpha=.99)
         print('Logrank test for %s - %s vs %s\n' % (args.group, names[0], names[1]))
         results.print_summary()
-    elif ngroup > 1:
-        fig, ax = plt.subplots(nrows=ngroup, ncols=1, sharex=True)
-        for i, rst in rmst:
-            lifelines.plotting.rmst_plot(rst, ax=ax)
-        fig.savefig(os.path.join(args.image_dir,'RMST_%s.png' % args.title))
 else:
     kmf.fit(df[args.time], df[args.status])
     kmf.plot_survival_function(ax=ax)
     ax.set_title(args.title)
     fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title))
+    print('#### No grouping variable, so no log rank or other Kaplan-Meier statistical output is available')
 if len(args.cphcols) > 0:
     fig, ax = plt.subplots()
-    ax.set_title('Cox PH model: %s' % args.title)
+    ax.set_title('Cox-PH model: %s' % args.title)
     cphcols = args.cphcols.strip().split(',')
     cphcols = [x.strip() for x in cphcols]
     notfound = sum([(x not in df.columns) for x in cphcols])
     if notfound > 0:
         sys.stderr.write('## CRITICAL USAGE ERROR (not a bug!): One or more requested Cox PH columns %s not found in supplied column header %s' % (args.cphcols, df.columns))
         sys.exit(6)
+    colsdf = df[cphcols]
     print('### Lifelines test of Proportional Hazards results with %s as covariates on %s' % (', '.join(cphcols), args.title))
-    cphcols += [args.time, args.status]
-    cphdf = df[cphcols]
+    cutcphcols = [args.time, args.status] + cphcols
+    cphdf = df[cutcphcols]
+    ucolcounts = colsdf.nunique(axis=0)
     cph.fit(cphdf, duration_col=args.time, event_col=args.status)
     cph.print_summary()
+    for i, cov in enumerate(colsdf.columns):
+         if ucolcounts[i] > 10:
+             v = pd.Series.tolist(cphdf[cov].quantile(QVALS))
+             vdt = df.dtypes[cov]
+             if vdt == 'int64':
+                 v = trimlegend(v)
+             axp = cph.plot_partial_effects_on_outcome(cov, cmap='coolwarm', values=v)
+             axp.set_title('Cox-PH %s quintile partials: %s' % (cov,args.title))
+             figr = axp.get_figure()
+             oname = os.path.join(args.image_dir,'%s_CoxPH_%s.%s' % (args.title, cov, args.image_type))
+             figr.savefig(oname)
+         else:
+             v = pd.unique(cphdf[cov])
+             v = [str(x) for x in v]
+             try:
+                 axp = cph.plot_partial_effects_on_outcome(cov, cmap='coolwarm', values=v)
+                 axp.set_title('Cox-PH %s partials: %s' % (cov,args.title))
+                 figr = axp.get_figure()
+                 oname = os.path.join(args.image_dir,'%s_CoxPH_%s.%s' % (args.title, cov, args.image_type))
+                 figr.savefig(oname)
+             except:
+                 pass
     cphaxes = cph.check_assumptions(cphdf, p_value_threshold=0.01, show_plots=True)
     for i, ax in enumerate(cphaxes):
         figr = ax[0].get_figure()
b
diff -r dd49a7040643 -r 232b874046a7 lifelines_tool/plotlykm.py
--- a/lifelines_tool/plotlykm.py Wed Aug 09 11:12:16 2023 +0000
+++ b/lifelines_tool/plotlykm.py Thu Aug 10 07:15:22 2023 +0000
[
@@ -2,8 +2,9 @@
 # km models for https://github.com/galaxyproject/tools-iuc/issues/5393
 # test as
 # python plotlykm.py --input_tab rossi.tab --htmlout "testfoo" --time "week" --status "arrest" --title "test" --image_dir images --cphcol="prio,age,race,paro,mar,fin"
+# Ross Lazarus July 2023
+import argparse
 
-import argparse
 import os
 import sys
 
@@ -13,15 +14,22 @@
 
 import pandas as pd
 
-# Ross Lazarus July 2023
 
+def trimlegend(v):
+    """
+    for int64 quintiles - must be ints - otherwise get silly legends with long float values
+    """
+    for i, av in enumerate(v):
+        x = int(av)
+        v[i] = str(x)
+    return v
 
 kmf = lifelines.KaplanMeierFitter()
 cph = lifelines.CoxPHFitter()
 
 parser = argparse.ArgumentParser()
 a = parser.add_argument
-a('--input_tab', default='', required=True)
+a('--input_tab', default='rossi.tab', required=True)
 a('--header', default='')
 a('--htmlout', default="test_run.html")
 a('--group', default='')
@@ -37,6 +45,7 @@
 df = pd.read_csv(args.input_tab, sep='\t')
 NCOLS = df.columns.size
 NROWS = len(df.index)
+QVALS = [.2, .4, .6, .8] # for partial cox ph plots
 defaultcols = ['col%d' % (x+1) for x in range(NCOLS)]
 testcols = df.columns
 if len(args.header.strip()) > 0:
@@ -57,62 +66,79 @@
     else:
         colsok = (args.time in defaultcols) and (args.status in defaultcols)
         if colsok:
-            sys.stderr.write('replacing first row of data derived header %s with %s' % (testcols, defaultcols))
+            print('Replacing first row of data derived header %s with %s' % (testcols, defaultcols))
             df.columns = defaultcols
         else:
             sys.stderr.write('## CRITICAL USAGE ERROR (not a bug!): time %s and status %s do not match anything in the file header, supplied header or automatic default column names %s' % (args.time, args.status, defaultcols))
-print('## Lifelines tool starting.\nUsing data header =', df.columns, 'time column =', args.time, 'status column =', args.status)
+print('## Lifelines tool\nInput data header =', df.columns, 'time column =', args.time, 'status column =', args.status)
 os.makedirs(args.image_dir, exist_ok=True)
 fig, ax = plt.subplots()
 if args.group > '':
     names = []
     times = []
     events = []
-    rmst = []
     for name, grouped_df in df.groupby(args.group):
         T = grouped_df[args.time]
         E = grouped_df[args.status]
         gfit = kmf.fit(T, E, label=name)
         kmf.plot_survival_function(ax=ax)
-        rst = lifelines.utils.restricted_mean_survival_time(gfit)
-        rmst.append(rst)
         names.append(str(name))
         times.append(T)
         events.append(E)
+    ax.set_title(args.title)
+    fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title))
     ngroup = len(names)
     if  ngroup == 2: # run logrank test if 2 groups
         results = lifelines.statistics.logrank_test(times[0], times[1], events[0], events[1], alpha=.99)
-        print(' vs '.join(names), results)
+        print('Logrank test for %s - %s vs %s\n' % (args.group, names[0], names[1]))
         results.print_summary()
-    elif ngroup > 1:
-        fig, ax = plt.subplots(nrows=ngroup, ncols=1, sharex=True)
-        for i, rst in rmst:
-            lifelines.plotting.rmst_plot(rst, ax=ax)
-        fig.savefig(os.path.join(args.image_dir,'RMST_%s.png' % args.title))
 else:
     kmf.fit(df[args.time], df[args.status])
     kmf.plot_survival_function(ax=ax)
-fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title))
+    ax.set_title(args.title)
+    fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title))
+    print('#### No grouping variable, so no log rank or other Kaplan-Meier statistical output is available')
 if len(args.cphcols) > 0:
     fig, ax = plt.subplots()
+    ax.set_title('Cox-PH model: %s' % args.title)
     cphcols = args.cphcols.strip().split(',')
     cphcols = [x.strip() for x in cphcols]
     notfound = sum([(x not in df.columns) for x in cphcols])
     if notfound > 0:
         sys.stderr.write('## CRITICAL USAGE ERROR (not a bug!): One or more requested Cox PH columns %s not found in supplied column header %s' % (args.cphcols, df.columns))
         sys.exit(6)
+    colsdf = df[cphcols]
     print('### Lifelines test of Proportional Hazards results with %s as covariates on %s' % (', '.join(cphcols), args.title))
-    cphcols += [args.time, args.status]
-    cphdf = df[cphcols]
+    cutcphcols = [args.time, args.status] + cphcols
+    cphdf = df[cutcphcols]
+    ucolcounts = colsdf.nunique(axis=0)
     cph.fit(cphdf, duration_col=args.time, event_col=args.status)
     cph.print_summary()
+    for i, cov in enumerate(colsdf.columns):
+         if ucolcounts[i] > 10: # a hack - assume categories are sparse - if not imaginary quintiles will have to do
+             v = pd.Series.tolist(cphdf[cov].quantile(QVALS))
+             vdt = df.dtypes[cov]
+             if vdt == 'int64':
+                 v = trimlegend(v)
+             axp = cph.plot_partial_effects_on_outcome(cov, cmap='coolwarm', values=v)
+             axp.set_title('Cox-PH %s quintile partials: %s' % (cov,args.title))
+             figr = axp.get_figure()
+             oname = os.path.join(args.image_dir,'%s_CoxPH_%s.%s' % (args.title, cov, args.image_type))
+             figr.savefig(oname)
+         else:
+             v = pd.unique(cphdf[cov])
+             v = [str(x) for x in v]
+             try:
+                 axp = cph.plot_partial_effects_on_outcome(cov, cmap='coolwarm', values=v)
+                 axp.set_title('Cox-PH %s partials: %s' % (cov,args.title))
+                 figr = axp.get_figure()
+                 oname = os.path.join(args.image_dir,'%s_CoxPH_%s.%s' % (args.title, cov, args.image_type))
+                 figr.savefig(oname)
+             except:
+                 pass
     cphaxes = cph.check_assumptions(cphdf, p_value_threshold=0.01, show_plots=True)
     for i, ax in enumerate(cphaxes):
         figr = ax[0].get_figure()
         titl = figr._suptitle.get_text().replace(' ','_').replace("'","")
         oname = os.path.join(args.image_dir,'CPH%s.%s' % (titl, args.image_type))
         figr.savefig(oname)
-
-
-
-
b
diff -r dd49a7040643 -r 232b874046a7 lifelines_tool/run_log.txt
--- a/lifelines_tool/run_log.txt Wed Aug 09 11:12:16 2023 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
[
@@ -1,107 +0,0 @@
-## Lifelines tool starting.
-Using data header = Index(['Unnamed: 0', 'week', 'arrest', 'fin', 'age', 'race', 'wexp', 'mar',
-       'paro', 'prio'],
-      dtype='object') time column = week status column = arrest
-### Lifelines test of Proportional Hazards results with prio, age, race, paro, mar, fin as covariates on test
-<lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations>
-             duration col = 'week'
-                event col = 'arrest'
-      baseline estimation = breslow
-   number of observations = 432
-number of events observed = 114
-   partial log-likelihood = -659.00
-         time fit was run = 2023-08-09 00:18:43 UTC
-
----
-            coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
-covariate                                                                                                         
-prio        0.10       1.10       0.03             0.04             0.15                 1.04                 1.16
-age        -0.06       0.94       0.02            -0.10            -0.02                 0.90                 0.98
-race        0.32       1.38       0.31            -0.28             0.92                 0.75                 2.52
-paro       -0.09       0.91       0.20            -0.47             0.29                 0.62                 1.34
-mar        -0.48       0.62       0.38            -1.22             0.25                 0.30                 1.29
-fin        -0.38       0.68       0.19            -0.75            -0.00                 0.47                 1.00
-
-            cmp to     z      p   -log2(p)
-covariate                                 
-prio          0.00  3.53 <0.005      11.26
-age           0.00 -2.95 <0.005       8.28
-race          0.00  1.04   0.30       1.75
-paro          0.00 -0.46   0.65       0.63
-mar           0.00 -1.28   0.20       2.32
-fin           0.00 -1.98   0.05       4.40
----
-Concordance = 0.63
-Partial AIC = 1330.00
-log-likelihood ratio test = 32.77 on 6 df
--log2(p) of ll-ratio test = 16.39
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-The ``p_value_threshold`` is set at 0.01. Even under the null hypothesis of no violations, some
-covariates will be below the threshold by chance. This is compounded when there are many covariates.
-Similarly, when there are lots of observations, even minor deviances from the proportional hazard
-assumption will be flagged.
-
-With that in mind, it's best to use a combination of statistical tests and visual tests to determine
-the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
-and looking for non-constant lines. See link [A] below for a full example.
-
-<lifelines.StatisticalResult: proportional_hazard_test>
- null_distribution = chi squared
-degrees_of_freedom = 1
-             model = <lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations>
-         test_name = proportional_hazard_test
-
----
-           test_statistic    p  -log2(p)
-age  km              6.99 0.01      6.93
-     rank            7.40 0.01      7.26
-fin  km              0.02 0.90      0.15
-     rank            0.01 0.91      0.13
-mar  km              1.64 0.20      2.32
-     rank            1.80 0.18      2.48
-paro km              0.06 0.81      0.31
-     rank            0.07 0.79      0.34
-prio km              0.92 0.34      1.57
-     rank            0.88 0.35      1.52
-race km              1.70 0.19      2.38
-     rank            1.68 0.19      2.36
-
-
-1. Variable 'age' failed the non-proportional test: p-value is 0.0065.
-
-   Advice 1: the functional form of the variable 'age' might be incorrect. That is, there may be
-non-linear terms missing. The proportional hazard test used is very sensitive to incorrect
-functional forms. See documentation in link [D] below on how to specify a functional form.
-
-   Advice 2: try binning the variable 'age' using pd.cut, and then specify it in `strata=['age',
-...]` in the call in `.fit`. See documentation in link [B] below.
-
-   Advice 3: try adding an interaction term with your time variable. See documentation in link [C]
-below.
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-
-   Bootstrapping lowess lines. May take a moment...
-
-
----
-[A]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html
-[B]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it
-[C]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates
-[D]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form
-[E]  https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification
-
b
diff -r dd49a7040643 -r 232b874046a7 lifelines_tool/test-data/readme_sample
--- a/lifelines_tool/test-data/readme_sample Wed Aug 09 11:12:16 2023 +0000
+++ b/lifelines_tool/test-data/readme_sample Thu Aug 10 07:15:22 2023 +0000
b
@@ -22,7 +22,7 @@
    number of observations = 432
 number of events observed = 114
    partial log-likelihood = -659.00
-         time fit was run = 2023-08-09 07:43:37 UTC
+         time fit was run = 2023-08-10 05:49:04 UTC
 
 ---
             coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%