veg / flea-pipeline Goto Github PK

• The current nesting in files like sequences.json and trees.json assumes longitudinal sampling. Switch to a flat, general format.
• Replace frequencies.json with coordinates.json.
• Make a run overview json file.

BioPython doesn't like some codons with X in them. For instance, translating "XAA" fails but "NAA" gives "X".

...if a timepoint has fewer than 3 consensus sequences.

In 6479e2e I added an exception if this condition occurs, to make the problem more obvious. Really, we should fix the script, though.

• parse file
• save task runtimes. boxplot or similar visualization
• save task start/stop times. plot time line.

Their input files are not correct

• number of changes per sequence
• location of the changes
• change in diversity of set of sequences

...but now we use a seperate file for copynumbers.

If the user provides an alignment, the pipeline will fail. They need to also provide a copynumber file. Also have the option to assume a copynumber of 1 for each sequence.

Even though it has three dependencies, sometimes compute_copynumber does not wait for them to complete.

• Newest version of pysam fails to build in the FLEA virtual environment
• flea-env's copy of bam2msa.py does not work. Temporarily using system version.
• submit pull request for our version of BioExt to veg/BioExt.
• write some tests for bealign.
• fix bealign gap extension penalty bug. It should be converted to a float.

So insertions relative to reference structure are visible in PV.

poly-A filtering and alignment to reference steps may shorten a CCS so that it no longer passes the minimum length filter.

Detect this with a new length filter after shift correction, which will ensure the sequences are de-gapped.

Result checks currently raise generic exceptions, such as when a result file is empty. Those should be caught and reported by the Result class.

Record the start and stop time of each task. Use them to generate benchmarks and identify slow tasks.

Julia dependencies are not currently defined.

Current dependencies:

• ArgParse
• Dates
• JSON
• DataStructures

After quality filtering, but before alignment.

Running FUBAR on the PC76 data resulted in different-sized arrays for each timepoint in rates.json. Each should have the same number of columns, equal to the number of columns in the multiple sequence alignment.

This causes short consensus sequences with long insertions or deletions to get overestimated copy numbers.

A long poly-A head can make the tail trimmer identify the whole sequence as tail, and vice-versa. Use a five-state HMM and Viterbi to implement trimming correctly.

• write HMM
• do not hardcode parameters
• rewrite in Julia or Python

Currently we are doing a full MSA for each cluster, only to generate the consensus sequence.

This will avoid errors like forgetting to pass length-filtered sequences downstream.

• quality filtering pipeline
• clustering pipeline
• copynumber pipeline
• msa pipeline

Name them after their cluster number, not after a random sequence in the cluster.

Log the following information:

• infiles
• stdout
• stderr
• external commands
• whether the task failed
• exceptions and tracebacks

Latest version takes gaps into account.

Currently the bad column summaries only show the top amino acid in each alignment, but these tend to be the same. Instead show the top three.

Could the qsub_files directory not exist when first task of the pipeline finishes?

Until we find a way to completely package the pipeline to include all third party dependencies, we will need to define all binaries used in the config for the sake of proper version control and explicit environment configuration.

• generate frequencies.json and turnover.json
• change file locations for evo history script
• add evo history script to pipeline
• change file locations for fubar script
• Add FUBAR script to generate rates.json
• rename subpipeline
• refactor main pipeline file
• remove old subpipeline
• add subdirectories for each subpipeline
• ensure new pipeline works for short amplicons, or at least have an option to disable parts that do not.
• update config file
• copynumbers in evo history script

For each step, show how many sequences were kept, and how many were discarded.

Deleting HQCSs during alignment editing should be reflected in the copynumbers, but currently it is not because they are done separately. Need to refactor so we do copynumber counting with the results of the edited alignment.

Package column-wise divergence and amino acid summaries into a json file.

Currently some tasks run locally. They are computationally light, but they should still be rewritten to run remotely.

The following can be done now, as the command-line scripts already exist:

• trim heads/tails
• shift correction
• translate
• backtranslate
• consensus

The rest will have to be rewritten as scripts.

• resolve translation of codon-aligned CCSs
• better metric: use chi-squared distance or earthmover's distance for each column's counts
• improve the plots: distribution of distances
• plot all time points combined: box plot for each column

The shift correction algorithm may split and correct within codons. Until we come up with a good general-purpose algorithm, at least we can do shift correction in codon-length chunks.

Currently assumes all input paths are absolute. Handle relative paths.

• print discarded sequences and the reason for the discard
• codon grouping with lookahead/behind
• use codon scoring matrix in out-of-frame regions.
• identify stutters
• possibly tolerate more errors in less conserved regions
• penalize early stop codons
• use base-wise quality scores.
• use bealign to do real codon alignment. possibly modify bealign to generate in-frame references.

• directory structure
• setup.py
• name all scripts in setup.py
• make scripts executable, but do not rely on doctests, because nose ignores excutables.
• do not find scripts with file
• update flea.config.tpl.

If the entire sequence is poly-A or poly-T, it should get dropped.

This should not happen, because the nucleotide alignment from which it is generated is a backtranslated amino-acid alignment. Possibly the backtranslation method is wrong.

The majority of sequences in the earliest timepoint have a 7 bp insertion after HXB2 coordinate 132. For some reason this insertion does not appear in the MRCA.