Add generateLphyScripts process for creating lphy scripts

nextflow/simulation.nf

@@ -133,6 +133,67 @@ process runLCUBE {
     """
 }
 
+process generateLphyScripts {
+    // Note: long-term this needs to be updated to include a docker image for reproducibility purposes
+    //        for now, set up your own R environment with the necessary packages.
+
+    input:
+    path seqs
+    val seed
+    val prefix
+
+    output:
+    path "${prefix}-${seed}.lphy", emit: lphy_script
+
+    shell:
+    """
+    cat > !{prefix}-!{seed}.lphy << EOL
+    data {
+    // Specify options for reading the date information from file
+    options = {ageDirection="forward", ageRegex=".*_.*_(\d*\.\d+|\d+\.\d*)$"};
+    D = readNexus(file="data/simulated.nex", options=options);
+
+    // Retrieve number of sites in alignment and taxa names
+    L = D.nchar();
+    taxa = D.taxa();
+
+    // Extract Trait information and count # of traits
+    D_trait = extractTrait(taxa=taxa, sep="_", i=0);
+    K = D_trait.canonicalStateCount();
+
+    // Specify the data type for traits and number of pairwise combos
+    dim = K*(K-1)/2;
+    dataType = D_trait.dataType();
+    }
+    model {
+    // Model evolutionary history
+    Π ~ Dirichlet(conc=[2.0, 2.0, 2.0, 2.0]);      // Nucleotide Frequency prior
+    κ ~ LogNormal(meanlog=1.0, sdlog=1.25);        // Transition/transversion ratio prior
+    γ ~ LogNormal(meanlog=0.0, sdlog=2.0);         // Shape parameter for discretized gamma
+    r ~ DiscretizeGamma(shape=γ, ncat=4, 
+            replicates=L);                           // Site rates
+    Θ ~ LogNormal(meanlog=0.0, sdlog=1.0);         // Effective population size prior
+    ψ ~ Coalescent(taxa=taxa, theta=Θ);            // Coalescent time scaled phylogenetic tree prior
+    D ~ PhyloCTMC(Q=hky(kappa=κ, freq=Π), mu=0.004, 
+            siteRates=r, tree=ψ);                    // MCMC for evolution of alignment
+
+    // Model the geographic history
+    // modeled via P(t) = e^(Dt)
+    π_trait ~ Dirichlet(conc=rep(element=3.0, times=K));   // Base frequencies, assuming symmetry
+    I ~ Bernoulli(minSuccesses=dim-2, p=0.5, 
+            replicates=dim);                                 // Determines which rates are 0. 
+                                                            // This, plus select, implements BSSVS.
+    R_trait ~ Dirichlet(conc=rep(element=1.0, times=dim)); // Off-diagonal entries of the rate matrix.
+    Q_trait = generalTimeReversible(rates=select(x=R_trait, 
+                indicator=I), freq=π_trait);               // intantaneous rate matrix
+    μ_trait ~ LogNormal(meanlog=0, sdlog=1.25);            // migration events per unit time
+    D_trait ~ PhyloCTMC(L=1, Q=Q_trait, dataType=dataType, 
+                mu=μ_trait, tree=ψ);                       // MCMC for discrete locations
+    }
+    EOL
+    """
+}
+
 process runBEAST {
     // Note: long-term this needs to be updated for a BEAST 2 docker image with the requisite lphy package, but for now, set up your own BEAST 2 environment.