fredjaya / topofind

Finding topologies for MAST

TopoFind

Finding tree topologies from sequence alignments.

Dependencies

iqtree 2.2.3.hmmster
R 4.1.2
devtools 2.4.3
MixtureModelHMM
AMAS

Installation

git clone https://github.com/fredjaya/topofind.git
cd topofind
conda env create -f env.yml
conda activate topofind
Rscript install_pacakges.R

• Install IQ-TREE manually and add symlink to /bin

Usage

• Infer the (best number of) topologies from a sequence alignment (test1.fa):

TopoFind.py -a data/test1.fa

• Set the number of threads (e.g. 4) to use for IQ-TREE:

TopoFind.py -a data/test1.fa -nt 4

• Show all available options:

TopoFind.py -h  

• Run unit tests:

python3 -m unittest tests/*

topofind's Issues

second iteration of t1_r2 on split blocks overwrites the initial t1_r2 run

Need better way of prefix handling when writing files

sh argument list too long at t=8

.command.sh: line 3: /home/frederickjaya/GitHub/rec-mast/bin/update_run_names.py: Argument list too long

Heaps of ways to partition by t=8. Easy fix is to set ENV ulimit before running nextflow. Better fix is to output each new name in a new channel.

`t1_iqtree_per_split()` only makes one tree when passing `val prefix`

Error thrown in concatenate_trees_for_mast() where only one tree exists.

Only occurs when passing val prefix:

process t1_iqtree_per_split {
    publishDir "${params.out}/${prefix}", mode: "copy"
    
    input:
        val prefix
        path splitted_aln
        val nthreads

    output:
        path '*.model.gz'
        path '*.treefile'
        path '*.log'     
        path '*.iqtree'  

    script:
    """
    iqtree2 -s ${splitted_aln} -pre ${splitted_aln.simpleName} -nt ${nthreads}
    """

workflow {

    t1_iqtree_per_split(params.prefix, split_aln.out[0].flatten(), params.nthreads)

}

Replace MixtureModelHMM with HMMSTER

Add catch for n trees before second mast iteration

Remake trees using MAST-assigned partitions instead of +R2 assignments

When splitting the alignment to b+1 blocks, each occurrence of block splitting is handled separately so that only t+1 trees are input to MAST per step.

When running MAST with t>=4 trees as input, There will be block regions generated from upstream splits that can be reused to avoid rebuilding trees.

However, the sites used to construct single trees from a given block will be reassigned by the HMM after the MAST process, hence the sites are likely to differ between the ones used to construct the tree.

I think trees should be reconstructed with the post-HMM site assignment, even though the same block is shared. I hope this makes sense... (for others and future me)

Use `val params.aln_path` for the output prefix

Instead of calling the file and using ${params.aln.simpleName}.

Add scripts to process HMMSTER outputs

• Generate partition file
• Generate site assignments
• Contiguous block size histogram

Output report

First, add MAST BIC to existing models.

See: https://carpentries-incubator.github.io/workflows-nextflow/12-reporting/index.html#templates

---

workflow.onComplete {collectFile} might be an ok option to collate BICs regardless of where the pipeline terminates.

package ‘MixtureModelHMM’ is not available for this version of R

install.packages("MixtureModelHMM")
Installing package into ‘/Users/lixingguang/Library/R/x86_64/4.1/library’
(as ‘lib’ is unspecified)
Warning in install.packages :
package ‘MixtureModelHMM’ is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

BIC format error after t=5 mast runs

On concat.fa and set2_paragone data. Might be slurm/scheduling issue

multiprocessing.pool.RemoteTraceback:                                                          
"""                                                                                            
Traceback (most recent call last):                                                             
  File "/home/u1070770/.conda/envs/rec-mast/lib/python3.10/multiprocessing/pool.py", line 125, 
in worker                                                                                      
    result = (True, func(*args, **kwds))                                                       
  File "/home/u1070770/.conda/envs/rec-mast/lib/python3.10/multiprocessing/pool.py", line 51, i
n starmapstar          
    return list(itertools.starmap(args[0], args[1]))                                           
  File "/home/u1070770/rec-mast/run.py", line 133, in mast                                     
    bic=get_bic(iqtree_out_path)               
  File "/home/u1070770/rec-mast/run.py", line 48, in get_bic                                   
    return(float(bic))                         
ValueError: could not convert string to float: ''                                              
"""                                            
                                              

Add catch for lack of informative sites during tree construction

Tree construction fails when partitioned alignment has no informative sites i.e. partition is tiny. This isn't a bad thing, but a proper catch needs to be added to account for this.

Error:

Caused by:
  Process `iqtree_default (1)` terminated with an error exit status (2)

Command executed:

  iqtree2 -s test2_class_1-out.fas -pre test2_class_1-out -nt 1

Command exit status:
  2

Command error:
  ERROR: Only one state is observed in alignment

To reproduce:

nextflow run ~/GitHub/rec-mast/main.nf --aln_path test2.phy --out `pwd`

Test that multiple trees fit data better than one tree

• Make single-tree with best-fitting model
• Input to MAST
• Record BIC (to compare with t>1 models)
• Success

Account for duplicate topologies

When inputting identical topologies into MAST and using the fully unlinked model, sites will have different weights despite the same topology. Also no need to store the same topology again

Speed up program

Ideas to make the program finish faster.

Currently, only a single instance of either split_aln() or mast() processes are run at once. Slow, but easy to control resource usage and debugging. IQ-TREE threads can be defined through python CLI with -T.

• Use multiprocessing.pool
• Continue iterations with the best performing model (will only work when discarding duplicate topologies #14)

-Benchmark the use of -T AUTO vs. arbitrary number i.e. -T 8

• #15 (comment)

Depending on infrastructure, change --num_threads to be --total_threads across all processes
Same with memory usage

• Just use mp.pool but perhaps add flag, so it only runs on a cluster

give useful names for output files and iterations

A standard naming convention would be nice.

concat_class_1-out_class-2-out_class-1-out_concat_mast_tr is unwieldy.