luca model change implemented

.devcontainer/devcontainer.json

@@ -1,6 +1,7 @@
 {
     // Build development container from pre-existing image
-    "image": "rachelicr/env-pisca-v1:latest",
+    //"image": "rachelicr/env-pisca-v1:latest",
+    "image": "icrsc/pisca-box:latest",
 
     // Set *default* container specific settings.json values on container create.
     "settings": {},

pisca-box-vue/app/static/phyfum-ages.csv

@@ -1,9 +1,9 @@
 taxon,age
-taxon1,45
-taxon2,45
-taxon3,45
-taxon4,45
+taxon1,41
+taxon2,42
+taxon3,43
+taxon4,44
 taxon5,45
-taxon6,45
-taxon7,45
-taxon8,45
+taxon6,46
+taxon7,47
+taxon8,48

pisca-box-vue/libs/cls_datatype.py

@@ -101,19 +101,19 @@ def all_logs(self,priors,operators):
         #create list from priors and operators
         logs = []
         for pr in prs:
-            prr = [pr]
-            if "|" in pr:
-                prr = pr.split("|")
+            prr = pr.split("|")            
             for pri in prr:
-                if pri not in logs:
-                    logs.append(pri)
+                pris = pri.split(",")
+                for prii in pris:
+                    if pri not in logs:
+                        logs.append(pri)
         for op in ops:
-            opp = [op]
-            if "|" in op:
-                opp = op.split("|")
+            opp = op.split("|")            
             for opi in opp:
-                if opi not in logs:
-                    logs.append(opi)        
+                opis = opi.split(",")
+                for opi in opis:
+                    if opi not in logs:
+                        logs.append(opi)        
 
         #for lg in self.default_logs:
         #    if lg not in logs and "." not in lg:

pisca-box-vue/libs/cls_dt_acna.py

@@ -24,9 +24,7 @@ def _default_operators(self):
         ops.append(['wideExchange','treeModel','3.0','','',''])
         ops.append(['wilsonBalding','treeModel','3.0','','',''])
         ops.append(['scaleOperator','treeModel.rootHeight','5.0','0.75','',''])
-        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])
-        ops.append(['scaleOperator','luca_branch','1.0','0.2','',''])
-        ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights','5.0','0.75','',''])
+        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])        
         return ops
     ##############################################################
     def _default_priors(self):

pisca-box-vue/libs/cls_dt_biallelic.py

@@ -24,9 +24,7 @@ def _default_operators(self):
         ops.append(['wideExchange','treeModel','3.0','','',''])
         ops.append(['wilsonBalding','treeModel','3.0','','',''])
         ops.append(['scaleOperator','treeModel.rootHeight','5.0','0.75','',''])
-        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])
-        ops.append(['scaleOperator','luca_branch','1.0','0.2','',''])        
-        ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights','5.0','0.75','',''])                
+        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])        
         ops.append(['scaleOperator','biallelicBinary.demethylation','0.25','0.25','',''])
         ops.append(['scaleOperator','biallelicBinary.homozygousMethylation','0.25','0.25','',''])
         ops.append(['scaleOperator','biallelicBinary.homozygousDemethylation','0.25','0.25','',''])

pisca-box-vue/libs/cls_dt_cnv.py

@@ -26,9 +26,7 @@ def _default_operators(self):
         ops.append(['wideExchange','treeModel','3.0','','',''])
         ops.append(['wilsonBalding','treeModel','3.0','','',''])
         ops.append(['scaleOperator','treeModel.rootHeight','5.0','0.75','',''])
-        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])
-        ops.append(['scaleOperator','luca_branch','1.0','0.2','',''])        
-        ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights','5.0','0.75','',''])                
+        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])        
         return ops    
     ##############################################################
     def _default_priors(self):

pisca-box-vue/libs/cls_dt_phyfum.py

@@ -26,9 +26,7 @@ def _default_operators(self):
         ops.append(['wideExchange','treeModel','3.0','','',''])
         ops.append(['wilsonBalding','treeModel','3.0','','',''])
         ops.append(['scaleOperator','treeModel.rootHeight','5.0','0.75','',''])
-        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])
-        ops.append(['scaleOperator','luca_branch','1.0','0.2','',''])
-        ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights','5.0','0.75','',''])
+        ops.append(['uniformOperator','treeModel.internalNodeHeights','30.0','','',''])        
         return ops
     ##############################################################
     def _default_priors(self):

pisca-box-vue/libs/cls_operators.py

@@ -4,7 +4,7 @@
 # ruff: noqa: F541
 
 class Operators(object):
-    def __init__(self,demographic,clocks,default_ops):                        
+    def __init__(self,demographic,luca_model,clocks,default_ops):                        
         self.ops = default_ops
         # operator, paramater, weight, scaleFactor, size, gaussian
         if demographic == "constant size":
@@ -15,6 +15,11 @@ def __init__(self,demographic,clocks,default_ops):
         if clocks["type"] == "random local clock":
             self.ops.append(['scaleOperator','localClock.relativeRates','15','0.75','',''])
             self.ops.append(['bitFlipOperator ','localClock.changes', '15','','',''])                
+        if luca_model == "Neoplastic progression":
+            self.ops.append(['scaleOperator','luca_branch','1.0','0.2','',''])
+            self.ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights,luca_branch','5.0','0.75','',''])
+        else:
+            self.ops.append(['upDownOperator',  'clock.rate|treeModel.allInternalNodeHeights','5.0','0.75','',''])
     ### PUBLIC INTERFACE #########        
     #---------------------------------------------------------------------------------
     def get_operators(self):
@@ -57,12 +62,21 @@ def _get_scaleOperator(self,op):
     #---------------------------------------------------------------------------------
     def _get_upDownOperator(self,op):
         operator, paramater, weight, scaleFactor, size, gaussian = op[0],op[1],op[2],op[3],op[4],op[5]
-        up_p = paramater.split("|")[0]
-        down_p = paramater.split("|")[1]
+        up_ps = paramater.split("|")[0].split(",")
+        down_ps = paramater.split("|")[1].split(",")
         op = ""
         op += f'\t<upDownOperator scaleFactor="{scaleFactor}" weight="{weight}">\n'
-        op += f'\t\t<up><parameter idref="{up_p}"/></up>\n'
-        op += f'\t\t<down><parameter idref="{down_p}"/></down>\n'
+
+        op += f'\t\t<up>\n'
+        for up_p in up_ps:        
+            op += f'\t\t\t<parameter idref="{up_p}"/>\n'
+        op += f'\t\t</up>\n'
+
+        op += f'\t\t<down>\n'
+        for down_p in down_ps:        
+            op += f'\t\t\t<parameter idref="{down_p}"/>\n'
+        op += f'\t\t</down>\n'
+
         op += f'\t</upDownOperator>\n'
         return op
     #---------------------------------------------------------------------------------

pisca-box-vue/libs/cls_xml.py

@@ -47,9 +47,11 @@ def get_xml(self):
         xml += self.fasta.get_cna_model()
         xml += self._add_comment("SITE MODEL")
         xml += self.fasta.get_site_model()
-        lh_val,lb_val,lb_up,lb_low = self.lucas["height"],self.lucas["branch"],self.lucas["upper"],self.lucas["lower"]
+        luca_model = self.lucas["model"]
+        hvalue,hupper,hlower = self.lucas["hvalue"],self.lucas["hupper"],self.lucas["hlower"]
+        bvalue,bupper,blower = self.lucas["bvalue"],self.lucas["bupper"],self.lucas["blower"]
         xml += self._add_comment("TREE LIKELIHOOD")
-        xml += self._get_tree_likelihood(lh_val,lb_val,lb_up,lb_low, self.clocks)
+        xml += self._get_tree_likelihood(luca_model,hvalue,hupper,hlower,bvalue,bupper,blower, self.clocks)
         xml += self._add_comment("OPERATORS")
         xml += self.operators.get_operators()
         xml += self._add_comment("MCMC and PRIORS")
@@ -188,25 +190,39 @@ def _get_rates(self):
         rts += '</rateCovarianceStatistic>\n'
         return rts                
     #----------------------------------
-    def _get_tree_likelihood(self,lh_val,lb_val,lb_up,lb_low,clocks):        
+    def _get_tree_likelihood(self,luca_model,hvalue,hupper,hlower,bvalue,bupper,blower,clocks):
         cen = ""
-        cen += '<cenancestorTreeLikelihood id="treeLikelihood" useAmbiguities="false">\n'
+
+        if luca_model == "Neoplastic progression":
+            cen += '<cenancestorTreeLikelihood id="treeLikelihood" useAmbiguities="false" heightRules="false" >\n'
+        else:
+            cen += '<cenancestorTreeLikelihood id="treeLikelihood" useAmbiguities="false" heightRules="true" >\n'
+
         cen += '\t<patterns idref="patterns"/>\n'
         cen += '\t<treeModel idref="treeModel"/>\n'
         cen += '\t<siteModel idref="siteModel"/>\n'
+
         cen += '\t<cenancestorHeight>\n'
-        cen += f'\t\t<parameter id="luca_height" value="{lh_val}"/>\n'
+        if luca_model == "Neoplastic progression":
+            cen += f'\t\t<parameter id="luca_height" value="{hvalue}" upper="{hupper}" lower="{hlower}" />\n'
+        else:
+            cen += f'\t\t<parameter id="luca_height" value="{hvalue}"" />\n'
         cen += '\t</cenancestorHeight>\n'
+
         cen += '\t<cenancestorBranch>\n'
-        cen += f'\t\t<parameter id="luca_branch" value="{lb_val}" upper="{lb_up}" lower="{lb_low}"/>\n'
+        if luca_model == "Neoplastic progression":            
+            cen += f'\t\t<parameter id="luca_branch" value="{bvalue}" upper="{bupper}" lower="{blower}" />\n'
+        else:
+            cen += f'\t\t<parameter id="luca_branch" value="{bvalue}" upper="{bupper}" lower="{blower}" />\n'
         cen += '\t</cenancestorBranch>\n'
+
         if clocks['type'] == "strict clock":
             cen += '\t<strictClockCenancestorBranchRates idref="branchRates"/>\n'
         elif clocks['type'] == "random local clock":
             cen += '\t<randomLocalClockModelCenancestor idref="branchRates"/>\n'
         cen += '</cenancestorTreeLikelihood>\n'
-        return cen        
-    #----------------------------------        
+        return cen
+    #----------------------------------
     def _get_report(self):
         rp = ""
         rp += '<report>\n'

pisca-box-vue/libs_pages/page_beauti.py

@@ -133,38 +133,51 @@ def add_widgets(include_header):
             with tabLuca:
                 st.subheader("Luca Height and Branch")
                 with st.container():
-                    cols = st.columns([1,1,3])
-                    with cols[0]:
-                        st.write("luca-upper")
-                    with cols[1]:
-                        min_age = st.number_input(label="luca-upper",value=min_age,label_visibility="collapsed")
-                    with cols[2]:
-                        st.caption("This is first branch height, or minimum age")
-                    cols = st.columns([1,1,3])
-                    with cols[0]:
-                        st.write("luca-height")
-                    with cols[1]:
-                        max_age = st.number_input(label="luca-height",value=max_age,label_visibility="collapsed")
-                    with cols[2]:
-                        st.caption("This is total tree height, or maximum age")
 
-                    cols = st.columns([5,1])
-                    with cols[0]:
-                        col0s = st.columns([4,1])
-                        with col0s[0]:
-                            lb_val = st.slider("luca-branch",
-                                            value = float(round(min_age/2)),
-                                            min_value = float(0),
-                                            max_value = float(min_age),
-                                            help="This can be between 0 and the root node, or minimum age")
-
-
-
+                    luca_model = st.radio("Luca model",["Neoplastic progression","Normal tissue"])
+                    st.write("---")
+                    hlower,hupper = -1,-1
+                    blower,bvalue=0,1
+                    if luca_model == "Neoplastic progression":
+                        height_rules = "false"
+                        cols = st.columns(3)
+                        with cols[0]:
+                            hlower = st.number_input("height lower DOLS-DOFD",value=max_age-min_age)
+                            blower = 0
+                            st.write("Branch lower")
+                            st.write(f"{blower}")
+                        with cols[1]:
+                            hvalue = st.number_input("height value DOLS-DOFD",value=max_age-min_age)
+                            bvalue = st.number_input("Branch value", value=bvalue)
+                        with cols[2]:
+                            hupper = st.number_input("Height upper DOLS-DOB",value=max_age)
+                            bupper = st.number_input("Branch upper DOLS-DOB",value=min_age)
+                    elif luca_model == "Normal tissue":
+                        height_rules = "true"
+                        cols = st.columns(3)
+                        with cols[0]:
+                            st.write("Height lower")
+                            st.write("NA")
+                            blower=0
+                            st.write("Branch lower")
+                            st.write("0")
+                        with cols[1]:
+                            hvalue = st.number_input("Height value DOLS-DOB", value=max_age)
+                            bvalue = st.number_input("Branch value", value=bvalue)
+                        with cols[2]:
+                            st.write("Height upper")
+                            st.write("NA")
+                            bupper = st.number_input("branch upper DOFS-DOB",value=min_age)                                                                                                                                                                                                                                                                                                                                                          
+                    st.write("---")                                                                                                                                                                                                                                                                                                                
                 lucas = {}
-                lucas["height"] = max_age
-                lucas["branch"] = lb_val
-                lucas["lower"] = 0.0
-                lucas["upper"] = min_age
+                lucas["model"] = luca_model
+                lucas["hvalue"] = hvalue
+                lucas["hupper"] = hupper
+                lucas["hlower"] = hlower
+                lucas["bvalue"] = bvalue
+                lucas["bupper"] = bupper
+                lucas["blower"] = blower
+
 
             ### TREES ########################################################
             with tabTrees:
@@ -195,38 +208,25 @@ def add_widgets(include_header):
                 prrs = prs.Priors(demographic,dt_obj.prs)
 
                 prior_types = ['oneOnX', 'logNormal', 'normal','exponential','uniform','laplace']
-                prh,prb = [],[]
-                with st.expander("View or change luca priors"):
-                    luca_priors = st.checkbox("Luca priors",value=True,help="Do you want to set priors on the luca branches?")
-                    prh = ['uniform','luca_height',str(lb_val),str(max_age),'','','','','','','']
-                    prb = ['uniform','luca_branch',str(lb_val),str(max_age),'','','','','','','']
-                    if luca_priors:
-                        cols = st.columns([1,1,1])
+                pri = []                                 
+                pri = ['uniform','luca_height',str(hlower),str(hupper),'','','','','','','']
+                use_prior = True
+                if luca_model == "Normal tissue":
+                    use_prior = False                                            
+                with st.expander("View or change luca height priors"):
+                    luca_priors = st.checkbox("Include luca-height priors",value=use_prior,help="Do you want to set priors on the luca branches?")                                        
+                    if luca_priors: 
+                        cols = st.columns(3)
                         with cols[0]:
-                            prior_type = st.selectbox('prior type', prior_types,key="luca",index=4)
-                            prior_type += "Prior"
+                            plower = st.number_input("lower",value=hlower)
                         with cols[1]:
-                            inc_height = st.checkbox("luca-height prior",value=True)
-                            inc_branch = st.checkbox("luca-prior prior",value=False)
-                        with cols[2]:
-                            luca_val = st.number_input(label="luca prior",value=round(lb_val,4))
-
-                        prh = [prior_type,'luca_height',str(luca_val),str(max_age),'','','','','','','']
-                        prb = [prior_type,'luca_branch',str(luca_val),str(max_age),'','','','','','','']
-                        if inc_height:
-                            prrs.update_one_prior(False,'luca_height',prh)
-                        else:
-                            prrs.update_one_prior(True,'luca_height',prh)
-
-                        if inc_branch:
-                            prrs.update_one_prior(False,'luca_branch',prb)
-                        else:
-                            prrs.update_one_prior(True,'luca_branch',prb)
-
+                            pupper = st.number_input("upper",value=hupper)                        
+                        pri = ['uniform','luca_height',str(plower),str(pupper),'','','','','','','']
+                        prrs.update_one_prior(False,'luca_height',pri)
                     else:
-                        prrs.update_one_prior(True,'luca_height',prh)
-                        prrs.update_one_prior(True,'luca_branch',prb)
-
+                        prrs.update_one_prior(True,'luca_height',[])
+
+                                                                          
                 if dtyp == 'biallelic':
                     with st.expander("View or change biallelic binary priors"):
                         bb_priors = st.checkbox("Biallelic binary priors",value=dtyp=='biallelic',help="Do you want to set priors on the biallelic binary?")
@@ -313,7 +313,7 @@ def add_widgets(include_header):
 
                 with st.expander("View or change all operators"):
                     ### OPERATORS #############################################################
-                    operators = ops.Operators(demographic,clocks,dt_obj.ops)
+                    operators = ops.Operators(demographic,luca_model,clocks,dt_obj.ops)
                     edited_df = st.data_editor(operators.get_as_dataframe(),num_rows="dynamic",use_container_width=True)
                     operators.update_from_dataframe(edited_df)