Tutorial: applying strainFlye to the "SheepGut" dataset¶

Usage: strainFlye [OPTIONS] COMMAND [ARGS]...

  Pipeline for the analysis of rare mutations in metagenomes.

  Please consult https://github.com/fedarko/strainFlye if you have any
  questions, comments, etc. about strainFlye. Thank you for using this tool!

Options:
  -v, --version  Show the version and exit.
  -h, --help     Show this message and exit.

Commands:
  align   Align reads to contigs, and filter the resulting alignment.
  call    [+] Call mutations in contigs naïvely & compute diversity indices.
  fdr     [+] Estimate and fix FDRs for contigs' naïve mutation calls.
  spot    [+] Identify putative mutational hotspots or coldspots.
  smooth  [+] Create and assemble smoothed and virtual reads.
  link    [+] Create link graphs showing co-occurring alleles.
  matrix  [+] Create codon and amino acid mutation matrices.
  dynam   [+] Compute simple information about growth dynamics.
  utils   [+] Miscellaneous utility commands provided with strainFlye.

Using strainFlye version "0.2.0-dev".
--------
utils gfa-to-fasta @ 0.00s: Starting...
Input GFA file: /Poppy/mfedarko/misc-data/sheepgut_flye_big_2.8_graph.gfa
Output FASTA file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs.fasta
--------
utils gfa-to-fasta @ 16.41s: Done.
Output FASTA file contains 78,793 sequences.

Usage: strainFlye align [OPTIONS] READS...

  Align reads to contigs, and filter the resulting alignment.

  Files of reads should be in the FASTA or FASTQ formats; GZIP'd files are
  allowed.

  This command involves multiple steps, including:

    1) Align reads to contigs (using minimap2) to generate a SAM file
    2) Convert this SAM file to a sorted and indexed BAM file
    3) Filter overlapping supplementary alignments (OSAs) from this BAM file
    4) Filter partially-mapped reads from this BAM file

  Note that we only sort the alignment file once, although we do re-index it
  after the two filtering steps. This decision is motivated by
  https://www.biostars.org/p/131333/#131335.

Options:
  -c, --contigs PATH              FASTA file of contigs to which reads will be
                                  aligned.  [required]
  -g, --graph PATH                GFA 1-formatted file describing an assembly
                                  graph of the contigs. This is used in the
                                  "filter partially-mapped reads" step to make
                                  the filter less strict for reads mapped to
                                  adjacent contigs in the graph. This isn't
                                  required; if it isn't passed, adjacent
                                  contigs will not be considered in this
                                  filter.  [default: (no graph)]
  -p, --minimap2-params TEXT      Additional parameters to pass to minimap2,
                                  besides the contig and read information.
                                  Depending on the size of your dataset and
                                  the amount of memory your system has, you
                                  may want to adjust the -I parameter; see
                                  minimap2's documentation for details. Please
                                  note that we do not perform any validation
                                  on this string before passing it to minimap2
                                  (so if you are allowing users to run
                                  strainFlye through a web server, be careful
                                  about shell injection).  [default: -ax asm20
                                  --secondary=no -I 8g --MD]
  -o, --output-dir DIRECTORY      Directory to which an output BAM file
                                  (final.bam) and BAM index file will be
                                  written. Some temporary files will also be
                                  written to this directory.  [required]
  --rm-tmp-bam / --no-rm-tmp-bam  Remove temporary (before filtering, and
                                  after just filtering OSAs) BAM files, to
                                  reduce space requirements.  [default: rm-
                                  tmp-bam]
  --verbose / --no-verbose        Display extra details during alignment
                                  filtering.  [default: no-verbose]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.2.0-dev".
--------
align @ 0.00s: Starting...
Input file(s) of reads: 
1. /Poppy/mkolmogo/sheep_meta/data/sheep_poop_CCS_dedup.fastq.gz
2. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200203_184822.Q20.fasta.gz
3. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200211_145859.Q20.fasta.gz
4. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200213_024013.Q20.fasta.gz
5. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200214_212611.Q20.fasta.gz
6. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200220_195233.Q20.fasta.gz
7. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200222_075656.Q20.fasta.gz
8. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200223_200916.Q20.fasta.gz
9. /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200227_201603.Q20.fasta.gz
Input contig file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs.fasta
Input graph file: /Poppy/mfedarko/misc-data/sheepgut_flye_big_2.8_graph.gfa
minimap2 parameters: -ax asm20 --secondary=no -I 8g --MD
Remove temporary BAM files?: Yes
Verbose?: No
Output directory: /Poppy/mfedarko/sftests/tutorial-output/alignment
--------
align @ 0.00s: Loading and checking contig information...
align @ 19.89s: Looks good. The FASTA file describes 78,793 contigs.
--------
align @ 19.89s: Sanity-checking that the assembly graph and FASTA file describe the same contigs (based on names and lengths)...
align @ 38.12s: Everything looks good so far.
--------
align @ 38.12s: Running minimap2 --> samtools view --> samtools sort...
align @ 38.12s: Exact command we're running: minimap2 -ax asm20 --secondary=no -I 8g --MD /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs.fasta /Poppy/mkolmogo/sheep_meta/data/sheep_poop_CCS_dedup.fastq.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200203_184822.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200211_145859.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200213_024013.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200214_212611.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200220_195233.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200222_075656.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200223_200916.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200227_201603.Q20.fasta.gz | samtools view -b - | samtools sort - -o /Poppy/mfedarko/sftests/tutorial-output/alignment/sorted-unfiltered.bam
[M::mm_idx_gen::67.042*1.68] collected minimizers
[M::mm_idx_gen::84.530*1.94] sorted minimizers
[M::main::84.583*1.94] loaded/built the index for 78793 target sequence(s)
[M::mm_mapopt_update::104.604*1.76] mid_occ = 100
[M::mm_idx_stat] kmer size: 19; skip: 10; is_hpc: 0; #seq: 78793
[M::mm_idx_stat::113.809*1.70] distinct minimizers: 466400347 (83.22% are singletons); average occurrences: 1.330; average spacing: 5.502
[M::worker_pipeline::234.045*2.35] mapped 60083 sequences
[[Note from Marcus: trimming down minimap2 outputs here to reduce notebook length...]]
[M::worker_pipeline::41791.167*3.09] mapped 46553 sequences
[M::main] Version: 2.17-r941
[M::main] CMD: minimap2 -ax asm20 --secondary=no -I 8g --MD /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs.fasta /Poppy/mkolmogo/sheep_meta/data/sheep_poop_CCS_dedup.fastq.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200203_184822.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200211_145859.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200213_024013.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200214_212611.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200220_195233.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200222_075656.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200223_200916.Q20.fasta.gz /Poppy/mkolmogo/sheep_meta/data/ccs_sequel_II/m54337U_200227_201603.Q20.fasta.gz
[M::main] Real time: 41792.733 sec; CPU: 129283.952 sec; Peak RSS: 30.432 GB
[bam_sort_core] merging from 526 files and 1 in-memory blocks...
align @ 45,699.81s: Done running minimap2 --> samtools view --> samtools sort.
--------
align @ 45,699.81s: Indexing the sorted BAM...
align @ 46,270.42s: Done indexing the sorted BAM.
--------
align @ 46,270.42s: Filtering reads with overlapping supplementary alignments (OSAs)...
align @ 46,938.50s: Done with pass 1 of the OSA filter; moving on to pass 2...
align @ 50,284.20s: Done filtering reads with overlapping supplementary alignments.
--------
align @ 50,284.24s: Indexing the OSA-filtered BAM...
align @ 50,868.70s: Done indexing the OSA-filtered BAM.
--------
align @ 50,868.70s: Removing the before-filtering BAM and its index to save space...
align @ 50,876.99s: Done removing these files.
--------
align @ 50,876.99s: Filtering partially-mapped reads...
align @ 50,877.11s: Loading assembly graph...
align @ 50,894.56s: Loaded assembly graph.
align @ 56,868.92s: Done filtering out partially-mapped reads.
--------
align @ 56,870.69s: Indexing the final BAM...
align @ 57,332.50s: Done indexing the final BAM.
--------
align @ 57,332.50s: Removing the just-OSA-filtered BAM and its index to save space...
align @ 57,341.75s: Done removing these files.
--------
align @ 57,341.76s: Done.

Found 468 contigs with lengths ≥ 1,000,000 bp. Took 12.52 sec.
Found all 3 of the selected contigs (CAMP, BACT1, BACT2).

Usage: strainFlye call [OPTIONS] COMMAND [ARGS]...

  [+] Call mutations in contigs naïvely & compute diversity indices.

  Consider a position "pos" in a contig. Using the alignment, we can count how
  many reads have a (mis)match operation to "pos" with one of the four
  nucleotides (A, C, G, T; we ignore degenerate nucleotides in reads). We
  represent these four nucleotides' counts at pos as follows:

      N1 = # reads of the most-common aligned nucleotide at pos,
      N2 = # reads of the second-most-common aligned nucleotide at pos,
      N3 = # reads of the third-most-common aligned nucleotide at pos,
      N4 = # reads of the fourth-most-common aligned nucleotide at pos.

  (We break ties arbitrarily.)

  strainFlye supports two types of naïve mutation calling based on these
  counts: p-mutations and r-mutations. These are described below.

  p-mutations (naïve percentage-based mutation calling)
  -----------------------------------------------------

  This takes as input some percentage p in the range (0%, 50%]. Define
  freq(pos) = N2 / (N1 + N2 + N3 + N4). This value, which is inherently
  constrained to the range [0%, 50%], is an estimate of the mutation frequency
  of this position. We classify pos as a p-mutation if freq(pos) ≥ p, AND if
  N2 ≥ the --min-alt-pos parameter (an integer representing a minimum number
  of reads that must support the alternate nucleotide).

  r-mutations (naïve read-count-based mutation calling)
  -----------------------------------------------------

  This takes as input some integer r > 0. We classify pos as an r-mutation if
  N2 ≥ r.

  Diversity indices
  -----------------

  Later on in the pipeline, if we perform FDR estimation on these called
  mutations using the "strainFlye fdr" module, we will need to select a "decoy
  contig" with relatively few mutations. A contig with low diversity indices
  may represent a promising decoy contig; so, for the sake of convenience,
  both commands output this information.

Options:
  -h, --help  Show this message and exit.

Commands:
  p-mutation  Call p-mutations and compute diversity indices.
  r-mutation  Call r-mutations and compute diversity indices.

Usage: strainFlye call p-mutation [OPTIONS]

  Call p-mutations and compute diversity indices.

  The primary parameter for this command is the lower bound of p, defined by
  --min-p. The BCF output will include "mutations" for all positions that pass
  this (likely very low) threshold; this BCF can be filtered using the
  utilities contained in the "strainFlye fdr" module.

Options:
  -c, --contigs PATH              FASTA file of contigs in which to naïvely
                                  call mutations. All contigs in this FASTA
                                  file should also be contained in the BAM
                                  file; it's ok if the BAM file contains
                                  contigs not in this FASTA file (we'll ignore
                                  them).  [required]
  -b, --bam PATH                  Sorted and indexed BAM file representing an
                                  alignment of reads to contigs.  [required]
  --min-p INTEGER RANGE           Minimum value of p for which to call
                                  p-mutations. This is scaled up by 100 (i.e.
                                  the default of 50 corresponds to 50 / 100 =
                                  0.5%) in order to bypass floating-point
                                  precision issues.  [default: 50; 0<x<=5000]
  --min-alt-pos INTEGER RANGE     In order for us to call a p-mutation at a
                                  position, this position's alternate
                                  nucleotide must be supported by at least
                                  this many reads.  [default: 2; x>=1]
  --div-index-p-list TEXT         List of values of p for which we'll compute
                                  diversity indices. These should all be
                                  separated by commas; and, as with --min-p,
                                  these are scaled up by 100. Please don't use
                                  commas as thousands separators.  [default:
                                  50,100,200,500,1000,2500,5000]
  -m, --min-read-number FLOAT RANGE
                                  Parameter that impacts the minimum
                                  (mis)match coverage needed in order to
                                  consider "counting" a position / mutation
                                  towards the diversity index. Given a value
                                  of p (converted to the range (0, 0.5]), a
                                  position must have a coverage of at least
                                  (--min-read-number / p) in order to be
                                  "sufficiently covered" and thus counted
                                  towards the diversity index.  [default: 5;
                                  x>0]
  -o, --output-dir DIRECTORY      Directory to which an output BCF file
                                  (describing the called mutations), BCF index
                                  file, and diversity index TSV file will be
                                  written. Some temporary files will also be
                                  written to this directory.  [required]
  --verbose / --no-verbose        Display extra details while running.
                                  [default: no-verbose]
  -h, --help                      Show this message and exit.

Usage: strainFlye fdr [OPTIONS] COMMAND [ARGS]...

  [+] Estimate and fix FDRs for contigs' naïve mutation calls.

Options:
  -h, --help  Show this message and exit.

Commands:
  estimate  Estimate the FDRs of contigs' mutation calls.
  fix       Fix contigs' mutation calls' estimated FDRs to an upper limit.

Usage: strainFlye fdr estimate [OPTIONS]

  Estimate the FDRs of contigs' mutation calls.

  We do this using the target-decoy approach (TDA). Given a set of C contigs,
  we select a "decoy contig" with relatively few called mutations. We then
  compute a mutation rate for this decoy contig, and use this mutation rate
  (along with the mutation rates of the other C - 1 "target" contigs) to
  estimate the FDRs of all of these target contigs' mutation calls.

  We can produce multiple FDR estimates for a single target contig's calls by
  varying the p or r threshold used (from the --min-p or --min-r threshold
  used to generate the input BCF file, up to the --high-p or --high-r
  threshold given here). Using this information (and information about the
  numbers of mutations called per megabase), we can plot an FDR curve for a
  given target contig's mutation calls.

  This command accepts an input BCF file of p- or r-mutations; however, in
  general we recommend using p-mutations (rather than r-mutations) for FDR
  estimation / fixing. Please see the Supplemental Material section of the
  strainFlye paper named "Identifying mutations based solely on read counts"
  for details.

Options:
  -c, --contigs PATH              FASTA file of contigs for which we will
                                  estimate mutation calls' FDRs. All contigs
                                  in this FASTA file should also be contained
                                  in the BAM and BCF files. It's ok if the BAM
                                  file contains contigs not in this FASTA file
                                  (we'll ignore them), but all contigs
                                  described in the BCF file's header must also
                                  be in this FASTA file.  [required]
  --bam PATH                      Sorted and indexed BAM file representing an
                                  alignment of reads to contigs.  [required]
  --bcf PATH                      Indexed BCF file describing naïvely called
                                  p- or r-mutations in the FASTA file's
                                  contigs, produced by "strainFlye call".
                                  [required]
  -di, --diversity-indices PATH   TSV file describing the diversity indices of
                                  a set of contigs. Used to automatically
                                  select a decoy contig. This option is
                                  mutually exclusive with --decoy-contig.
                                  [default: (nothing)]
  -dc, --decoy-contig TEXT        Name of a specific contig to use as the
                                  decoy contig for FDR estimation. This option
                                  is mutually exclusive with --diversity-
                                  indices.  [default: (nothing)]
  -dctx, --decoy-contexts [Full|CP2|Tv|Nonsyn|Nonsense|CP2Tv|CP2Nonsyn|CP2Nonsense|TvNonsyn|TvNonsense|CP2TvNonsense|Everything]
                                  "Context-dependent" types of positions
                                  and/or mutations to which the computation of
                                  mutation rates in the decoy contig will be
                                  limited. The "Full" option will consider the
                                  entire decoy contig; all other options will
                                  limit the computation to certain types of
                                  positions and/or potential mutations. (CP2
                                  will focus on positions in the second codon
                                  position of genes predicted in the decoy
                                  contig using Prodigal; Nonsyn will focus on
                                  potential nonsynonymous mutations in these
                                  genes; Nonsense will focus on potential
                                  nonsense mutations in these genes; and Tv
                                  will focus on potential transversion
                                  mutations.) You can specify this option
                                  multiple times to generate multiple sets of
                                  FDR estimates. (You can also specify
                                  "Everything" to use all available contexts.)
                                  [default: CP2; required]
  -hp, --high-p INTEGER RANGE     (Only applies if the input BCF file
                                  describes p-mutations.) p-mutations with a
                                  mutation rate (freq(pos)) greater than or
                                  equal to this are considered "indisputable,"
                                  and will not be included in FDR estimation.
                                  Like the values of p used as parameters to
                                  "strainFlye call p-mutation", this is in the
                                  range (0, 5000], such that h = N corresponds
                                  to (N / 100)%. Corresponds to the
                                  "highFrequency" threshold mentioned in the
                                  paper.  [default: 500; 0<x<=5000]
  -hr, --high-r INTEGER RANGE     (Only applies if the input BCF file
                                  describes r-mutations.) r-mutations with an
                                  alternate nucleotide read coverage greater
                                  than or equal to this are considered
                                  indisputable, and will not be included in
                                  FDR estimation.  [default: 100; x>0]
  -dml, --decoy-min-length INTEGER RANGE
                                  Minimum length of a potential decoy contig
                                  (in bp). Only used if --diversity-indices is
                                  specified.  [default: 1000000; x>0]
  -dmac, --decoy-min-average-coverage FLOAT RANGE
                                  Minimum average (mis)match coverage of a
                                  potential decoy contig. Only used if
                                  --diversity-indices is specified.  [default:
                                  500; x>0]
  -o, --output-dir DIRECTORY      Directory to which TSV files describing the
                                  estimated FDRs and numbers of naïvely called
                                  mutations per megabase will be written. In
                                  all TSV files, rows correspond to target
                                  contigs and columns correspond to p or r
                                  values. We will generate one estimated FDR
                                  TSV file for every time --decoy-context is
                                  specified, and just one number-of-naïvely-
                                  called-mutations-per-megabase file.
                                  [required]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
fdr estimate @ 0.00s: Starting...
Input contig file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs_atleast_1Mbp.fasta
Input BAM file: /Poppy/mfedarko/sftests/tutorial-output/alignment/final.bam
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/call-p15/naive-calls.bcf
Input diversity indices file: /Poppy/mfedarko/sftests/tutorial-output/call-p15/diversity-indices.tsv
Manually-set decoy contig: None
Decoy contig context-dependent position / mutation type(s): ('Everything',)
High p threshold (only used if the BCF describes p-mutations): 500
High r threshold (only used if the BCF describes r-mutations): 100
Minimum length of a potential decoy contig (only used if diversity indices are specified): 1,000,000 bp
Minimum average coverage of a potential decoy contig (only used if diversity indices are specified): 500.0x
Output directory: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr-info
--------
fdr estimate @ 0.00s: Loading and checking FASTA, BAM, and BCF files...
fdr estimate @ 5.42s: The FASTA file describes 468 contig(s).
fdr estimate @ 5.50s: All FASTA contig(s) are included in the BAM file (this BAM file has 78,793 reference(s)).
fdr estimate @ 19.55s: The FASTA file's contig(s) and BCF file's contig(s) match.
fdr estimate @ 19.55s: Also, the input BCF file describes p-mutations (minimum p = 15).
fdr estimate @ 19.56s: The lengths of all contig(s) in the FASTA file match the corresponding lengths in the BAM and BCF files.
fdr estimate @ 19.56s: So far, these files seem good.
--------
fdr estimate @ 19.56s: Selecting a decoy contig based on the diversity indices...
fdr estimate @ 19.58s: Selected edge_6104 as the decoy contig.
fdr estimate @ 19.58s: Verified that this decoy contig is present in the FASTA, BAM, and BCF files.
fdr estimate @ 19.58s: (Sorry for doubting you.)
--------
fdr estimate @ 19.58s: Determining range of value(s) of p to consider...
fdr estimate @ 19.58s: We'll consider 485 value(s) of p: from 15 to 499.
fdr estimate @ 19.58s: p-mutations for p ≥ 500 will be considered "indisputable."
fdr estimate @ 19.58s: These "indisputable" mutations won't be included in the FDR estimation results.
--------
fdr estimate @ 19.58s: Computing mutation rates for edge_6104 at these threshold values, for each of the 11 decoy context(s)...
fdr estimate @ 21.67s: At least one of the decoy context(s) requires us to have gene predictions for edge_6104; running Prodigal...
-------------------------------------
PRODIGAL v2.6.3 [February, 2016]         
Univ of Tenn / Oak Ridge National Lab
Doug Hyatt, Loren Hauser, et al.     
-------------------------------------
Request:  Single Genome, Phase:  Training
Reading in the sequence(s) to train...1289244 bp seq created, 32.65 pct GC
Locating all potential starts and stops...30910 nodes
Looking for GC bias in different frames...frame bias scores: 2.68 0.19 0.12
Building initial set of genes to train from...done!
Creating coding model and scoring nodes...done!
Examining upstream regions and training starts...done!
-------------------------------------
Request:  Single Genome, Phase:  Gene Finding
Finding genes in sequence #1 (1289244 bp)...done!
fdr estimate @ 27.01s: Finished running Prodigal! It predicted 1,297 gene(s) in edge_6104.
fdr estimate @ 27.83s: At least one of the decoy context(s) requires us to know "unreasonable" positions (where reference and consensus disagree) in edge_6104; going through alignment...
fdr estimate @ 804.88s: Done: identified 0 unreasonable position(s) in edge_6104.
fdr estimate @ 804.88s: Computing mutation rates for decoy context "CP2" (1 / 11)...
fdr estimate @ 804.89s: Done computing "CP2" mutation rates for edge_6104.
fdr estimate @ 804.89s: Computing mutation rates for decoy context "CP2Nonsense" (2 / 11)...
fdr estimate @ 811.65s: Done computing "CP2Nonsense" mutation rates for edge_6104.
fdr estimate @ 811.65s: Computing mutation rates for decoy context "CP2Nonsyn" (3 / 11)...
fdr estimate @ 819.02s: Done computing "CP2Nonsyn" mutation rates for edge_6104.
fdr estimate @ 819.02s: Computing mutation rates for decoy context "CP2Tv" (4 / 11)...
fdr estimate @ 819.03s: Done computing "CP2Tv" mutation rates for edge_6104.
fdr estimate @ 819.03s: Computing mutation rates for decoy context "CP2TvNonsense" (5 / 11)...
fdr estimate @ 825.73s: Done computing "CP2TvNonsense" mutation rates for edge_6104.
fdr estimate @ 825.73s: Computing mutation rates for decoy context "Full" (6 / 11)...
fdr estimate @ 825.75s: Done computing "Full" mutation rates for edge_6104.
fdr estimate @ 825.75s: Computing mutation rates for decoy context "Nonsense" (7 / 11)...
fdr estimate @ 845.39s: Done computing "Nonsense" mutation rates for edge_6104.
fdr estimate @ 845.39s: Computing mutation rates for decoy context "Nonsyn" (8 / 11)...
fdr estimate @ 866.92s: Done computing "Nonsyn" mutation rates for edge_6104.
fdr estimate @ 866.92s: Computing mutation rates for decoy context "Tv" (9 / 11)...
fdr estimate @ 867.05s: Done computing "Tv" mutation rates for edge_6104.
fdr estimate @ 867.05s: Computing mutation rates for decoy context "TvNonsense" (10 / 11)...
fdr estimate @ 886.23s: Done computing "TvNonsense" mutation rates for edge_6104.
fdr estimate @ 886.23s: Computing mutation rates for decoy context "TvNonsyn" (11 / 11)...
fdr estimate @ 907.39s: Done computing "TvNonsyn" mutation rates for edge_6104.
fdr estimate @ 907.41s: Done computing mutation rates for all decoy context(s) for edge_6104.
--------
fdr estimate @ 907.41s: Computing mutation rates and FDR estimates for the 467 target contig(s)...
fdr estimate @ 941.67s: Done.
--------
fdr estimate @ 941.72s: Done.

fdr-CP2Nonsense.tsv    fdr-CP2Tv.tsv	 fdr-TvNonsense.tsv
fdr-CP2Nonsyn.tsv      fdr-Full.tsv	 fdr-TvNonsyn.tsv
fdr-CP2.tsv	       fdr-Nonsense.tsv  fdr-Tv.tsv
fdr-CP2TvNonsense.tsv  fdr-Nonsyn.tsv	 num-mutations-per-mb.tsv

Usage: strainFlye fdr fix [OPTIONS]

  Fix contigs' mutation calls' estimated FDRs to an upper limit.

  This takes as input the estimated FDRs from "strainFlye fdr estimate" (if
  you used multiple decoy contexts, then you will need to choose which set of
  FDR estimates to use here) to guide us on how to fix the FDR for each
  contig. Note that mutations that passed the "high" p or r threshold
  specified for "strainFlye fdr estimate", and thus were not used for FDR
  estimation, will all be included in the output BCF file from this command;
  these mutations are considered "indisputable."

  We include indisputable mutations from the decoy contig and from all target
  contigs in our output BCF file. We will only consider including non-
  indisputable mutations from the target contigs: the decision of which non-
  indisputable mutations will be included is based on the lowest p or r
  parameter for a target contig that yields an estimated FDR ≤ the fixed FDR
  given here.

Options:
  -b, --bcf PATH            Indexed BCF file describing naïvely called p- or
                            r-mutations, produced by "strainFlye call".
                            [required]
  -fi, --fdr-info FILE      One of the estimated FDR TSV files produced by
                            "strainFlye fdr estimate".  [required]
  --fdr FLOAT RANGE         False discovery rate at which identified mutations
                            will be fixed. This is interpreted as a
                            percentage: the default of 1 corresponds to an FDR
                            of 1%. We do not restrict this to an upper limit,
                            since estimated FDRs can technically exceed 100%.
                            [default: 1; x>0]
  -o, --output-bcf FILE     Filepath to which an output indexed BCF file
                            (describing the called mutations at the optimal p
                            or r value for each contig, along with any
                            "indisputable" mutations relative to the --high-p
                            or --high-r parameter used for strainFlye fdr
                            estimate) will be written.  [required]
  --verbose / --no-verbose  Display extra details while writing the filtered
                            BCF.  [default: no-verbose]
  -h, --help                Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
fdr fix @ 0.00s: Starting...
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/call-p15/naive-calls.bcf
Input FDR estimate file: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr-info/fdr-CP2.tsv
FDR at which to fix mutation calls: 1.00%
Verbose?: No
Output BCF file with mutation calls at the fixed FDR: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr1pct.bcf
--------
fdr fix @ 0.00s: Loading and checking BCF and TSV files...
fdr fix @ 17.34s: Looks good so far; decoy contig seems to be edge_6104.
fdr fix @ 17.34s: Looks like the cutoff for "indisputable" mutations was p = 500.
fdr fix @ 17.34s: All mutations passing this cutoff will be included in the output BCF file.
--------
fdr fix @ 17.34s: Based on the FDR information, finding optimal values of p for each contig...
fdr fix @ 17.36s: Done.
fdr fix @ 17.36s: For 239 / 467 contigs, there exist values of p (at least, considering the range from p = 15 to p = 499) that yield estimated FDRs ≤ 1.0%.
fdr fix @ 17.36s: These values range from p = 37 (edge_1300) to p = 475 (edge_1016).
fdr fix @ 17.36s: The mean of these values is p = 255.96.
--------
fdr fix @ 17.36s: Writing a filtered BCF file (to a temporary location, for now) including both (1) indisputable mutations from all contigs and (2) non-indisputable mutations from the target contigs that result in a FDR ≤ 1.00%...
fdr fix @ 32.55s: Done.
--------
fdr fix @ 32.55s: Updating the filtered BCF file's header to indicate that FDR fixing was done...
fdr fix @ 36.15s: Done.
--------
fdr fix @ 36.15s: Indexing the BCF file...
fdr fix @ 37.04s: Done indexing the BCF file.
--------
fdr fix @ 37.06s: Done.

Usage: strainFlye spot [OPTIONS] COMMAND [ARGS]...

  [+] Identify putative mutational hotspots or coldspots.

  Many methods exist for identifying these sorts of hotspots or coldspots;
  strainFlye's implementations of these methods are intended mostly as a quick
  proof-of-concept for replicating the results shown in our paper, and are not
  extremely "feature-rich" quite yet.

Options:
  -h, --help  Show this message and exit.

Commands:
  hot-features  Identify hotspot features (for example, genes).
  cold-gaps     Identify long coldspot "gaps" without any mutations.

Usage: strainFlye spot hot-features [OPTIONS]

  Identify hotspot features (for example, genes).

  By "feature", we refer to a single continuous region within a contig, as
  described in the file given for --features. These regions could describe
  anything: predicted protein-coding genes, introns or exons, intergenic
  regions of interest, etc. For now, we treat each feature independently (e.g.
  we don't lump together exons from the same "Parent" gene; each feature is
  considered separately as a potential "hotspot").

  You can configure whether or not we classify a feature as a hotspot by
  adjusting the --min-num-mutations and --min-perc-mutations parameters; at
  least one of these parameters must be specified. If both parameters are
  specified, then both checks (number of mutations in a feature, and
  percentage of mutations in a feature) will need to pass in order for us to
  label a feature as a hotspot.

Options:
  -b, --bcf PATH                  Indexed BCF file describing single-
                                  nucleotide mutations in a set of contigs.
                                  [required]
  -f, --features PATH             Generic Feature Format version 3 (GFF3) file
                                  describing "features" (e.g. predicted
                                  protein-coding genes) in the contigs. We
                                  will ignore features located on contigs that
                                  are not described in the BCF file's header.
                                  [required]
  -mn, --min-num-mutations INTEGER RANGE
                                  Label a feature as a "hotspot" if it
                                  contains at least this many mutations.
                                  [default: (no check); x>0]
  -mp, --min-perc-mutations FLOAT RANGE
                                  Label a feature as a "hotspot" if its
                                  percentage of mutations ((# mutations /
                                  feature length) × 100) is at least this
                                  value.  [default: (no check); 0<x<=100]
  -o, --output-hotspots FILE      Filepath to which an output tab-separated
                                  values (TSV) file describing hotspot
                                  features across all contigs will be written.
                                  [required]
  -h, --help                      Show this message and exit.

-------------------------------------
PRODIGAL v2.6.3 [February, 2016]         
Univ of Tenn / Oak Ridge National Lab
Doug Hyatt, Loren Hauser, et al.     
-------------------------------------
Request:  Metagenomic, Phase:  Training
Initializing training files...done!
-------------------------------------
Request:  Metagenomic, Phase:  Gene Finding
Finding genes in sequence #1 (2153394 bp)...done!
Finding genes in sequence #2 (2806161 bp)...done!
Finding genes in sequence #3 (1289244 bp)...done!

Using strainFlye version "0.1.0-dev".
--------
spot hot-features @ 0.00s: Starting...
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr1pct.bcf
Input feature file: /Poppy/mfedarko/sftests/tutorial-output/prodigal_anonymous_selected3.gff
Minimum number of mutations needed to call a feature a hotspot: 5
Minimum % of mutations needed to call a feature a hotspot: 2.0
Output file describing hotspot features: /Poppy/mfedarko/sftests/tutorial-output/hotspot-features-n5p2.tsv
--------
spot hot-features @ 0.00s: Loading and checking the BCF file...
spot hot-features @ 12.41s: Looks good so far.
--------
spot hot-features @ 12.41s: Going through features in the GFF3 file and identifying hotspot features...
spot hot-features @ 12.64s: Identified 165 hotspot feature(s) across all 468 contigs in the BCF file.
--------
spot hot-features @ 12.64s: Writing out this information to a TSV file...
spot hot-features @ 12.64s: Done.
--------
spot hot-features @ 12.65s: Done.

Usage: strainFlye spot cold-gaps [OPTIONS]

  Identify long coldspot "gaps" without any mutations.

  To clarify, we define a "gap" of length L on a contig as a collection of
  continuous positions [N, N + 1, ... N + L - 2, N + L - 1] in which no
  positions are mutations (based on the input BCF file).

  If the --circular flag is specified, then we can loop around the contig from
  right to left; otherwise, the left and right sides of the contig are hard
  boundaries. To give an example of this, consider a 9-nucleotide contig that
  has mutations at positions 4 and 6:

                             Mut.    Mut.
                  1   2   3   4   5   6   7   8   9

  If --circular is specified, then this contig has two gaps: one gap of length
  1 (covering just position 5, between the two mutations), and another gap of
  length 6 (starting at position 7 and looping around to position 3: [7, 8, 9,
  1, 2, 3]).

  If --circular is not specified, then this contig has three gaps: [1, 2, 3],
  [5], and [7, 8, 9].

Options:
  -b, --bcf PATH                  Indexed BCF file describing single-
                                  nucleotide mutations in a set of contigs.
                                  [required]
  -l, --min-length INTEGER RANGE  Label a gap between mutations in a contig as
                                  a "coldspot" if the gap is at least this
                                  long (in bp).  [default: 5000; x>0;
                                  required]
  --circular / --no-circular      If --circular is specified, we'll assume
                                  that all contigs are circular: we'll
                                  consider the gap "looping around" from the
                                  rightmost mutation in each contig to the
                                  leftmost mutation in the contig as a
                                  potential coldspot. Otherwise, we will
                                  assume all contigs are linear.  [default:
                                  no-circular]
  --exact-pvals / --no-exact-pvals
                                  If --exact-pvals is specified, we'll use the
                                  method for computing exact longest-gap
                                  p-values given in equation (3.1) of (Naus
                                  1982). We use Python's decimal library with
                                  default parameters to try to deal with large
                                  numbers, but we can't guarantee that this
                                  won't cause the program to fail in certain
                                  cases. If --exact-pvals is not specified,
                                  we'll use equation (3.3) of (Naus 1982),
                                  which gives an approximation of the p-value.
                                  [default: exact-pvals]
  -o, --output-coldspots FILE     Filepath to which an output tab-separated
                                  values (TSV) file describing coldspots will
                                  be written. The longest gap in each contig
                                  will also be given a p-value, indicating the
                                  probability of the longest gap being at
                                  least this long (given the length of and
                                  number of mutated positions in this contig,
                                  and assuming that mutations occur with a
                                  constant probability at any position in the
                                  contig). See (Naus 1982) for details.
                                  [required]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
spot cold-gaps @ 0.00s: Starting...
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr1pct.bcf
Minimum coldspot gap length: 5,000 bp
Check for circular coldspot gaps?: No
Compute exact longest-gap p-values?: Yes
Output file describing coldspot gaps: /Poppy/mfedarko/sftests/tutorial-output/coldspot-gaps-minlen5000-nocircular.tsv
--------
spot cold-gaps @ 0.00s: Loading and checking the BCF file...
spot cold-gaps @ 12.16s: Looks good so far.
--------
spot cold-gaps @ 12.16s: Going through contigs and identifying coldspot gaps...
spot cold-gaps @ 15.74s: Identified 16,154 coldspot gap(s) across all 468 contigs in the BCF file.
--------
spot cold-gaps @ 15.74s: Writing out this information to a TSV file...
spot cold-gaps @ 15.75s: Done.
--------
spot cold-gaps @ 15.76s: Done.

Usage: strainFlye smooth [OPTIONS] COMMAND [ARGS]...

  [+] Create and assemble smoothed and virtual reads.

Options:
  -h, --help  Show this message and exit.

Commands:
  create    Create smoothed and virtual reads for each contig.
  assemble  Assemble contigs' smoothed and virtual reads using LJA.

Usage: strainFlye smooth create [OPTIONS]

  Create smoothed and virtual reads for each contig.

Options:
  -c, --contigs PATH              FASTA file of contigs for which we will
                                  create smoothed and virtual reads. All
                                  contigs in this FASTA file should also be
                                  contained in the BAM and BCF files; it's ok
                                  if the BAM or BCF files contain contigs not
                                  in this FASTA file (we'll ignore them).
                                  [required]
  --bam PATH                      Sorted and indexed BAM file representing an
                                  alignment of reads to contigs.  [required]
  --bcf PATH                      Indexed BCF file describing single-
                                  nucleotide mutations in a set of contigs.
                                  [required]
  -di, --diversity-indices PATH   TSV file describing the diversity indices of
                                  a set of contigs, produced by one of the
                                  "strainFlye call" commands. Only used if
                                  --virtual-reads is specified. Along with
                                  diversity indices, these files list contigs'
                                  average coverages. If --virtual-reads is
                                  specified, we will need to know contigs'
                                  average coverages. So, if a diversity
                                  indices file is provided here, then we can
                                  avoid re-computing average coverages.
                                  (Otherwise, if --virtual-reads is specified
                                  but no diversity indices file is provided,
                                  we'll need to compute average coverages;
                                  this will take some extra time.)  [default:
                                  (nothing)]
  --virtual-reads / --no-virtual-reads
                                  If --virtual-reads is specified, we'll
                                  create virtual reads covering "low-coverage"
                                  regions in contigs.  [default: virtual-
                                  reads]
  -vrp, --virtual-read-well-covered-perc FLOAT RANGE
                                  Only used if --virtual-reads is specified.
                                  In a contig with average coverage (based on
                                  the BAM file, and only considering match +
                                  mismatch counts) C, we will define a
                                  position in this contig (with coverage P,
                                  based only on the smoothed reads) as low-
                                  coverage if ((P / C) × 100) is less than
                                  this percentage.  [default: 95; 0<=x<=100;
                                  required]
  -vrf, --virtual-read-flank INTEGER RANGE
                                  Only used if --virtual-reads is specified.
                                  When we add virtual reads spanning a single
                                  continuous low-coverage region, these reads
                                  will start and end this many positions
                                  before and after the region. (For example,
                                  the default of 100 means that the virtual
                                  reads created for a low-coverage region of
                                  5,000 bp will all be 5,200 bp.)  [default:
                                  100; x>=0; required]
  -o, --output-dir DIRECTORY      Directory to which smoothed reads (and
                                  virtual reads, if --virtual-reads is
                                  specified) will be written. Each contig's
                                  reads will be written to a gzipped FASTA
                                  file in this directory named
                                  [contig].fasta.gz.  [required]
  --verbose / --no-verbose        Display extra details while running.
                                  [default: no-verbose]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
smooth create @ 0.00s: Starting...
Input contig file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs_selected3.fasta
Input BAM file: /Poppy/mfedarko/sftests/tutorial-output/alignment/final.bam
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr1pct.bcf
Input diversity indices file: /Poppy/mfedarko/sftests/tutorial-output/call-p15/diversity-indices.tsv
Create virtual reads?: Yes
"Well-covered" percentage used for virtual read creation: 95.00%
Create "flanks" of this length on either side of each virtual read: 100
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/smooth/reads
--------
smooth create @ 0.00s: Loading and checking FASTA, BAM, and BCF files...
smooth create @ 0.10s: The FASTA file describes 3 contig(s).
smooth create @ 0.19s: All FASTA contig(s) are included in the BAM file (this BAM file has 78,793 reference(s)).
smooth create @ 12.33s: All FASTA contig(s) are included in the BCF file (the header of this BCF file describes 468 contig(s)).
smooth create @ 12.34s: The lengths of all contig(s) in the FASTA file match the corresponding lengths in the BAM and BCF files.
smooth create @ 12.34s: So far, these files seem good.
--------
smooth create @ 12.34s: All contigs must be > (2 × --virtual-read-flank) = 200 bp long. Checking this...
smooth create @ 12.34s: All contigs meet this minimum length.
--------
smooth create @ 12.34s: Going through contigs and creating smoothed and virtual reads...
smooth create @ 12.34s: On contig edge_1671 (2,153,394 bp) (1 / 3 contigs = 33.33%).
smooth create @ 12.39s: Contig edge_1671 has 23,821 mutated position(s).
smooth create @ 5,086.59s: From the 262,757 linear alignment(s) to contig edge_1671: created 228,902 smoothed read(s) and ignored 33,855 linear alignment(s).
smooth create @ 5,087.21s: Contig edge_1671 (average coverage 1,415.55x, based on the BAM file) has 538 run(s) of consecutive low-coverage (≤ 1,344.77x) positions (based on smoothed read coverages). Creating virtual reads...
smooth create @ 5,542.74s: Created 43,596 virtual read(s) total for contig edge_1671.
smooth create @ 5,542.77s: On contig edge_2358 (2,806,161 bp) (2 / 3 contigs = 66.67%).
smooth create @ 5,542.78s: Contig edge_2358 has 17 mutated position(s).
smooth create @ 8,102.10s: From the 742,784 linear alignment(s) to contig edge_2358: created 719,654 smoothed read(s) and ignored 23,130 linear alignment(s).
smooth create @ 8,102.85s: Contig edge_2358 (average coverage 2,993.59x, based on the BAM file) has 612 run(s) of consecutive low-coverage (≤ 2,843.91x) positions (based on smoothed read coverages). Creating virtual reads...
smooth create @ 8,518.15s: Created 103,347 virtual read(s) total for contig edge_2358.
smooth create @ 8,518.17s: On contig edge_6104 (1,289,244 bp) (3 / 3 contigs = 100.00%).
smooth create @ 8,518.18s: Contig edge_6104 has 35 mutated position(s).
smooth create @ 10,215.96s: From the 476,373 linear alignment(s) to contig edge_6104: created 476,232 smoothed read(s) and ignored 141 linear alignment(s).
smooth create @ 10,216.38s: Contig edge_6104 (average coverage 4,158.79x, based on the BAM file) has 361 run(s) of consecutive low-coverage (≤ 3,950.85x) positions (based on smoothed read coverages). Creating virtual reads...
smooth create @ 10,398.82s: Created 80,090 virtual read(s) total for contig edge_6104.
smooth create @ 10,398.83s: Done.
--------
smooth create @ 10,398.93s: Done.

Usage: strainFlye smooth assemble [OPTIONS]

  Assemble contigs' smoothed and virtual reads using LJA.

  Please note that this command relies on the "simple_ec" branch of LJA being
  installed on your system. See strainFlye's README (and/or LJA's manual) for
  details on installing LJA.

Options:
  -r, --reads-dir DIRECTORY   Directory produced by "strainFlye smooth create"
                              containing smoothed (and optionally virtual)
                              reads. We will use LJA to assemble each
                              *.fasta.gz file in this directory (representing
                              reads from different contigs) independently.
                              [required]
  -p, --lja-params TEXT       Additional parameters to pass to LJA, besides
                              the --reads and --output-dir parameters. To
                              explain our defaults: the --simpleec flag is
                              currently only available on the simple_ec branch
                              of LJA. This flag tells LJA to perform "simple"
                              error correction by removing all edges in the de
                              Bruijn graph with k-mer coverage less than
                              --Cov-threshold (10), as well as reads passing
                              through these edges. This "simple" error
                              correction is done instead of running the
                              default LJA error correction module, mowerDBG.
                              Please note that we do not perform any
                              validation on this string before passing it to
                              LJA (so if you are allowing users to run
                              strainFlye through a web server, be careful
                              about shell injection).  [default: --simpleec
                              --Cov-threshold 10]
  -b, --lja-bin FILE          Location of LJA's "lja" binary file, which can
                              be used to run LJA. This should be located in
                              the bin/ directory constructed when you compiled
                              LJA. If this is not provided, we will check to
                              see if LJA is available in your $PATH.
                              [default: (Look for LJA's bin in the $PATH
                              environment variable)]
  -o, --output-dir DIRECTORY  Directory to which output LJA assemblies (one
                              subdirectory per *.fasta.gz file in the input
                              --reads-dir) will be written.  [required]
  --verbose / --no-verbose    Display extra details while running.  [default:
                              no-verbose]
  -h, --help                  Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
smooth assemble @ 0.00s: Starting...
Input reads directory: /Poppy/mfedarko/sftests/tutorial-output/smooth/reads
LJA parameters: --simpleec --Cov-threshold 10
LJA binary (if None, we'll look in $PATH): /home/mfedarko/software/LJA-branch/bin/lja
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies
--------
smooth assemble @ 0.00s: Assembling each *.fasta.gz file in /Poppy/mfedarko/sftests/tutorial-output/smooth/reads...
smooth assemble @ 0.00s: Found file edge_1671.fasta.gz, presumably for contig edge_1671. Assembling.
smooth assemble @ 0.00s: Running this command: /home/mfedarko/software/LJA-branch/bin/lja --reads /Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_1671.fasta.gz --simpleec --Cov-threshold 10 --output-dir /Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671
00:00:00 2Mb  INFO: Hello! You are running La Jolla Assembler (LJA), a tool for genome assembly from PacBio HiFi reads
00:00:10 5Mb  INFO: e213d10b906cc633b78a27eceb52d44073478d6d
00:00:10 5Mb  INFO: LJA pipeline started
00:00:10 5Mb  INFO: Performing initial correction with k = 5001
00:00:10 0Mb  INFO: Reading reads
00:00:10 0Mb  INFO: Extracting minimizers
00:00:22 3.1Gb  INFO: Finished read processing
00:00:22 3.1Gb  INFO: 3804650 hashs collected. Starting sorting.
00:00:22 3.1Gb  INFO: Finished sorting. Total distinct minimizers: 26485
00:00:22 3.1Gb  INFO: Starting construction of sparse de Bruijn graph
00:00:22 3.1Gb  INFO: Vertex map constructed.
00:00:22 3.1Gb  INFO: Filling edge sequences.
00:00:36 3.2Gb  INFO: Finished sparse de Bruijn graph construction.
00:00:36 3.2Gb  INFO:  Collecting tips 
00:00:36 3.2Gb  INFO: Added 1672 artificial minimizers from tips.
00:00:36 3.2Gb  INFO: Collected 52384 old edges.
00:00:36 3.2Gb  INFO: New minimizers added to sparse graph.
00:00:36 3.2Gb  INFO: Refilling graph with old edges.
00:00:37 3.2Gb  INFO: Filling graph with new edges.
00:00:37 3.2Gb  INFO: Finished fixing sparse de Bruijn graph.
00:00:37 3.2Gb  INFO: Starting to extract disjointigs.
00:00:37 3.2Gb  INFO: Finished extracting 2434 disjointigs of total size 18962614
00:00:37 0Mb  INFO: Loading disjointigs from file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/k5001/disjointigs.fasta"
00:00:38 35Mb  INFO: Filling bloom filter with k+1-mers.
00:00:39 35Mb  INFO: Filled 31425472 bits out of 217285760
00:00:39 35Mb  INFO: Finished filling bloom filter. Selecting junctions.
00:00:39 36Mb  INFO: Collected 5299 junctions.
00:00:39 36Mb  INFO: Starting DBG construction.
00:00:39 36Mb  INFO: Vertices created.
00:00:39 36Mb  INFO: Filled dbg edges. Adding hanging vertices 
00:00:39 36Mb  INFO: Added 1 hanging vertices
00:00:39 36Mb  INFO: Merging unbranching paths
00:00:39 36Mb  INFO: Ended merging edges. Resulting size 2748
00:00:40 36Mb  INFO: Cleaning edge coverages
00:00:40 36Mb  INFO: Collecting alignments of sequences to the graph
00:00:51 3.2Gb  INFO: Alignment collection finished. Total length of alignments is 677345
00:00:52 3.2Gb  INFO: Could not correct 4910 reads. They will be removed.
00:00:52 3.2Gb  INFO: Uncorrected reads were removed.
00:00:52 3.2Gb  INFO: Applying changes to the graph
00:00:55 3.2Gb  INFO: Printing reads to fasta file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/k5001/corrected_reads.fasta"
00:01:14 5Mb  INFO: Performing repeat resolution by transforming de Bruijn graph into Multiplex de Bruijn graph
00:01:14 0Mb  INFO: Loading graph from fasta
00:01:14 13Mb  INFO: Finished loading graph
00:01:15 126Mb  INFO: Looking for unique edges
00:01:15 126Mb  INFO: Marked 76 long edges as unique
00:01:15 126Mb  INFO: Marking extra edges as unique based on read paths
00:01:15 126Mb  INFO: Marked 81 edges as unique
00:01:15 126Mb  INFO: Splitting graph with unique edges
00:01:15 126Mb  INFO: Processing 100 components
00:01:15 126Mb  INFO: Finished unique edges search. Found 89 unique edges
00:01:15 126Mb  INFO: Analysing repeats of multiplicity 2 and looking for additional unique edges
00:01:15 126Mb  INFO: Finished processing of repeats of multiplicity 2. Found 0 erroneous edges.
00:01:15 126Mb  INFO: Resolving repeats
00:01:15 126Mb  INFO: Constructing paths
00:01:17 169Mb  INFO: Building graph
00:01:17 169Mb  INFO: Increasing k
00:01:17 170Mb  INFO: Finished increasing k
00:01:17 170Mb  INFO: Exporting remaining active transitions
00:01:17 170Mb  INFO: Export to Dot
00:01:17 170Mb  INFO: Export to GFA and compressed contigs
00:01:17 178Mb  INFO: Finished repeat resolution
00:01:17 5Mb  INFO: Performing polishing and homopolymer uncompression
00:01:18 6Mb  INFO: Aligning reads back to assembly
00:01:27 2.2Gb  INFO: Finished alignment.
00:01:27 2.2Gb  INFO: Printing alignments to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/uncompressing/alignments.txt"
00:01:27 2.2Gb  INFO: Reading and processing initial reads from ["/Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_1671.fasta.gz"]
00:02:45 3.4Gb  INFO: Uncompressing homopolymers in contigs
00:02:46 3.4Gb  INFO: Total zero covered nucleotides 0
00:02:46 3.4Gb  INFO: Calculating overlaps between adjacent uncompressed edges
00:02:46 3.4Gb  INFO: Printing final gfa file to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/mdbg.gfa"
00:02:46 3.4Gb  INFO: Printing final assembly to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/assembly.fasta"
00:02:46 5Mb  INFO: Final homopolymer compressed and corrected reads can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/k5001/corrected_reads.fasta"
00:02:46 5Mb  INFO: Final graph with homopolymer compressed edges can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/mdbg/mdbg.hpc.gfa"
00:02:46 5Mb  INFO: Final graph can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/mdbg.gfa"
00:02:46 5Mb  INFO: Final assembly can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_1671/assembly.fasta"
00:02:46 5Mb  INFO: LJA pipeline finished
smooth assemble @ 166.78s: Finished running LJA on file edge_1671.fasta.gz.
smooth assemble @ 166.78s: Found file edge_2358.fasta.gz, presumably for contig edge_2358. Assembling.
smooth assemble @ 166.78s: Running this command: /home/mfedarko/software/LJA-branch/bin/lja --reads /Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_2358.fasta.gz --simpleec --Cov-threshold 10 --output-dir /Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358
00:00:00 5Mb  INFO: Hello! You are running La Jolla Assembler (LJA), a tool for genome assembly from PacBio HiFi reads
00:00:09 5Mb  INFO: e213d10b906cc633b78a27eceb52d44073478d6d
00:00:09 5Mb  INFO: LJA pipeline started
00:00:09 5Mb  INFO: Performing initial correction with k = 5001
00:00:09 0Mb  INFO: Reading reads
00:00:09 0Mb  INFO: Extracting minimizers
00:00:45 5.4Gb  INFO: Finished read processing
00:00:45 5.4Gb  INFO: 10414993 hashs collected. Starting sorting.
00:00:46 5.4Gb  INFO: Finished sorting. Total distinct minimizers: 8259
00:00:46 5.4Gb  INFO: Starting construction of sparse de Bruijn graph
00:00:46 5.4Gb  INFO: Vertex map constructed.
00:00:46 5.4Gb  INFO: Filling edge sequences.
00:01:21 7.8Gb  INFO: Finished sparse de Bruijn graph construction.
00:01:21 7.8Gb  INFO:  Collecting tips 
00:01:21 7.8Gb  INFO: Added 4 artificial minimizers from tips.
00:01:21 7.8Gb  INFO: Collected 16526 old edges.
00:01:21 7.8Gb  INFO: New minimizers added to sparse graph.
00:01:21 7.8Gb  INFO: Refilling graph with old edges.
00:01:21 7.8Gb  INFO: Filling graph with new edges.
00:01:21 7.8Gb  INFO: Finished fixing sparse de Bruijn graph.
00:01:21 7.8Gb  INFO: Starting to extract disjointigs.
00:01:21 7.8Gb  INFO: Finished extracting 20 disjointigs of total size 2169195
00:01:21 0Mb  INFO: Loading disjointigs from file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/k5001/disjointigs.fasta"
00:01:21 13Mb  INFO: Filling bloom filter with k+1-mers.
00:01:22 14Mb  INFO: Filled 9578413 bits out of 66213600
00:01:22 14Mb  INFO: Finished filling bloom filter. Selecting junctions.

00:01:22 14Mb  INFO: Collected 779 junctions.
00:01:22 14Mb  INFO: Starting DBG construction.
00:01:22 14Mb  INFO: Vertices created.
00:01:22 14Mb  INFO: Filled dbg edges. Adding hanging vertices 
00:01:22 14Mb  INFO: Added 0 hanging vertices
00:01:22 14Mb  INFO: Merging unbranching paths
00:01:23 14Mb  INFO: Ended merging edges. Resulting size 16
00:01:23 14Mb  INFO: Cleaning edge coverages
00:01:23 14Mb  INFO: Collecting alignments of sequences to the graph
00:01:54 6.7Gb  INFO: Alignment collection finished. Total length of alignments is 662331
00:01:55 6.7Gb  INFO: Could not correct 3 reads. They will be removed.
00:01:55 6.7Gb  INFO: Uncorrected reads were removed.
00:01:55 6.7Gb  INFO: Applying changes to the graph
00:02:00 6.7Gb  INFO: Printing reads to fasta file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/k5001/corrected_reads.fasta"
00:03:04 5Mb  INFO: Performing repeat resolution by transforming de Bruijn graph into Multiplex de Bruijn graph
00:03:04 0Mb  INFO: Loading graph from fasta
00:03:04 11Mb  INFO: Finished loading graph
00:03:06 265Mb  INFO: Looking for unique edges
00:03:06 265Mb  INFO: Marked 8 long edges as unique
00:03:06 265Mb  INFO: Marking extra edges as unique based on read paths
00:03:06 265Mb  INFO: Marked 8 edges as unique
00:03:06 265Mb  INFO: Splitting graph with unique edges
00:03:06 265Mb  INFO: Processing 5 components
00:03:06 265Mb  INFO: Finished unique edges search. Found 22 unique edges
00:03:06 265Mb  INFO: Analysing repeats of multiplicity 2 and looking for additional unique edges
00:03:06 265Mb  INFO: Finished processing of repeats of multiplicity 2. Found 0 erroneous edges.
00:03:06 265Mb  INFO: Resolving repeats
00:03:06 265Mb  INFO: Constructing paths
00:03:11 496Mb  INFO: Building graph
00:03:11 496Mb  INFO: Increasing k
00:03:11 496Mb  INFO: Finished increasing k
00:03:11 496Mb  INFO: Exporting remaining active transitions
00:03:11 496Mb  INFO: Export to Dot
00:03:11 496Mb  INFO: Export to GFA and compressed contigs
00:03:11 0.5Gb  INFO: Finished repeat resolution
00:03:12 5Mb  INFO: Performing polishing and homopolymer uncompression
00:03:12 9Mb  INFO: Aligning reads back to assembly
00:03:38 5.8Gb  INFO: Finished alignment.
00:03:38 5.8Gb  INFO: Printing alignments to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/uncompressing/alignments.txt"
00:03:38 5.8Gb  INFO: Reading and processing initial reads from ["/Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_2358.fasta.gz"]
00:06:12 7Gb  INFO: Uncompressing homopolymers in contigs
00:06:12 7Gb  INFO: Total zero covered nucleotides 0
00:06:12 7Gb  INFO: Calculating overlaps between adjacent uncompressed edges
00:06:12 7Gb  INFO: Printing final gfa file to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/mdbg.gfa"
00:06:12 7Gb  INFO: Printing final assembly to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/assembly.fasta"
00:06:13 5Mb  INFO: Final homopolymer compressed and corrected reads can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/k5001/corrected_reads.fasta"
00:06:13 5Mb  INFO: Final graph with homopolymer compressed edges can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/mdbg/mdbg.hpc.gfa"
00:06:13 5Mb  INFO: Final graph can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/mdbg.gfa"
00:06:13 5Mb  INFO: Final assembly can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_2358/assembly.fasta"
00:06:13 5Mb  INFO: LJA pipeline finished
smooth assemble @ 539.88s: Finished running LJA on file edge_2358.fasta.gz.
smooth assemble @ 539.88s: Found file edge_6104.fasta.gz, presumably for contig edge_6104. Assembling.
smooth assemble @ 539.89s: Running this command: /home/mfedarko/software/LJA-branch/bin/lja --reads /Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_6104.fasta.gz --simpleec --Cov-threshold 10 --output-dir /Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104
00:00:00 2Mb  INFO: Hello! You are running La Jolla Assembler (LJA), a tool for genome assembly from PacBio HiFi reads
00:00:09 5Mb  INFO: e213d10b906cc633b78a27eceb52d44073478d6d
00:00:09 5Mb  INFO: LJA pipeline started
00:00:09 5Mb  INFO: Performing initial correction with k = 5001
00:00:09 0Mb  INFO: Reading reads
00:00:09 0Mb  INFO: Extracting minimizers
00:00:30 5.6Gb  INFO: Finished read processing
00:00:30 5.6Gb  INFO: 5736312 hashs collected. Starting sorting.
00:00:30 5.6Gb  INFO: Finished sorting. Total distinct minimizers: 3721
00:00:30 5.6Gb  INFO: Starting construction of sparse de Bruijn graph
00:00:30 5.6Gb  INFO: Vertex map constructed.
00:00:30 5.6Gb  INFO: Filling edge sequences.
00:00:50 5.7Gb  INFO: Finished sparse de Bruijn graph construction.
00:00:50 5.7Gb  INFO:  Collecting tips 
00:00:50 5.7Gb  INFO: Added 14 artificial minimizers from tips.
00:00:50 5.7Gb  INFO: Collected 7448 old edges.
00:00:50 5.7Gb  INFO: New minimizers added to sparse graph.
00:00:50 5.7Gb  INFO: Refilling graph with old edges.
00:00:50 5.7Gb  INFO: Filling graph with new edges.
00:00:50 5.7Gb  INFO: Finished fixing sparse de Bruijn graph.
00:00:50 5.7Gb  INFO: Starting to extract disjointigs.
00:00:50 5.7Gb  INFO: Finished extracting 26 disjointigs of total size 1074351
00:00:50 0Mb  INFO: Loading disjointigs from file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/k5001/disjointigs.fasta"
00:00:50 9Mb  INFO: Filling bloom filter with k+1-mers.
00:00:51 9Mb  INFO: Filled 4370802 bits out of 30218400
00:00:51 9Mb  INFO: Finished filling bloom filter. Selecting junctions.
00:00:51 9Mb  INFO: Collected 392 junctions.
00:00:51 9Mb  INFO: Starting DBG construction.
00:00:51 9Mb  INFO: Vertices created.
00:00:51 9Mb  INFO: Filled dbg edges. Adding hanging vertices 
00:00:51 9Mb  INFO: Added 0 hanging vertices
00:00:51 9Mb  INFO: Merging unbranching paths
00:00:51 9Mb  INFO: Ended merging edges. Resulting size 34
00:00:51 9Mb  INFO: Cleaning edge coverages
00:00:51 9Mb  INFO: Collecting alignments of sequences to the graph
00:01:11 5.7Gb  INFO: Alignment collection finished. Total length of alignments is 435418
00:01:11 5.7Gb  INFO: Could not correct 24 reads. They will be removed.
00:01:11 5.7Gb  INFO: Uncorrected reads were removed.
00:01:11 5.7Gb  INFO: Applying changes to the graph
00:01:14 5.7Gb  INFO: Printing reads to fasta file "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/k5001/corrected_reads.fasta"
00:01:50 5Mb  INFO: Performing repeat resolution by transforming de Bruijn graph into Multiplex de Bruijn graph
00:01:50 0Mb  INFO: Loading graph from fasta
00:01:50 7Mb  INFO: Finished loading graph
00:01:51 243Mb  INFO: Looking for unique edges
00:01:51 243Mb  INFO: Marked 6 long edges as unique
00:01:51 243Mb  INFO: Marking extra edges as unique based on read paths
00:01:51 243Mb  INFO: Marked 6 edges as unique
00:01:51 243Mb  INFO: Splitting graph with unique edges
00:01:51 243Mb  INFO: Processing 4 components
00:01:51 243Mb  INFO: Finished unique edges search. Found 10 unique edges
00:01:51 243Mb  INFO: Analysing repeats of multiplicity 2 and looking for additional unique edges
00:01:51 243Mb  INFO: Finished processing of repeats of multiplicity 2. Found 0 erroneous edges.
00:01:51 243Mb  INFO: Resolving repeats
00:01:51 243Mb  INFO: Constructing paths
00:01:54 326Mb  INFO: Building graph
00:01:54 326Mb  INFO: Increasing k
00:01:54 326Mb  INFO: Finished increasing k
00:01:54 326Mb  INFO: Exporting remaining active transitions
00:01:54 326Mb  INFO: Export to Dot
00:01:54 326Mb  INFO: Export to GFA and compressed contigs
00:01:54 337Mb  INFO: Finished repeat resolution
00:01:54 5Mb  INFO: Performing polishing and homopolymer uncompression
00:01:54 5Mb  INFO: Aligning reads back to assembly
00:02:10 3.5Gb  INFO: Finished alignment.
00:02:10 3.5Gb  INFO: Printing alignments to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/uncompressing/alignments.txt"
00:02:10 3.5Gb  INFO: Reading and processing initial reads from ["/Poppy/mfedarko/sftests/tutorial-output/smooth/reads/edge_6104.fasta.gz"]

00:03:52 3.6Gb  INFO: Uncompressing homopolymers in contigs
00:03:52 3.6Gb  INFO: Total zero covered nucleotides 0
00:03:52 3.6Gb  INFO: Calculating overlaps between adjacent uncompressed edges
00:03:52 3.6Gb  INFO: Printing final gfa file to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/mdbg.gfa"
00:03:52 3.6Gb  INFO: Printing final assembly to "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/assembly.fasta"
00:03:53 5Mb  INFO: Final homopolymer compressed and corrected reads can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/k5001/corrected_reads.fasta"
00:03:53 5Mb  INFO: Final graph with homopolymer compressed edges can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/mdbg/mdbg.hpc.gfa"
00:03:53 5Mb  INFO: Final graph can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/mdbg.gfa"
00:03:53 5Mb  INFO: Final assembly can be found here: "/Poppy/mfedarko/sftests/tutorial-output/smooth/assemblies/edge_6104/assembly.fasta"
00:03:53 5Mb  INFO: LJA pipeline finished
smooth assemble @ 772.97s: Finished running LJA on file edge_6104.fasta.gz.
smooth assemble @ 772.97s: Done.
--------
smooth assemble @ 772.97s: Done.

Usage: strainFlye link [OPTIONS] COMMAND [ARGS]...

  [+] Create link graphs showing co-occurring alleles.

  The "nt" command should be run first; this generates information about the
  counts and co-occurrences of nucleotides at mutated positions in a contig.

  The "graph" command takes as input the information produced by "nt" and
  creates link graphs (one link graph per contig) from it. There are many
  parameters that impact the graph creation, so we have split this into two
  steps in order to make it easy to rerun the "graph" step with different
  parameter settings if needed (the "nt" command will likely take longer to
  run).

Options:
  -h, --help  Show this message and exit.

Commands:
  nt     Compute (co-)occurrence information for mutations' nucleotides.
  graph  Convert (co-)occurrence information into a link graph structure.

Usage: strainFlye link nt [OPTIONS]

  Compute (co-)occurrence information for mutations' nucleotides.

Options:
  -c, --contigs PATH          FASTA file of contigs for which we will create
                              link graphs. All contigs in this FASTA file
                              should also be contained in the BAM and BCF
                              files; it's ok if the BAM or BCF files contain
                              contigs not in this FASTA file (we'll ignore
                              them).  [required]
  --bam PATH                  Sorted and indexed BAM file representing an
                              alignment of reads to contigs.  [required]
  --bcf PATH                  Indexed BCF file describing single-nucleotide
                              mutations in a set of contigs.  [required]
  -o, --output-dir DIRECTORY  Directory to which information about nucleotide
                              frequencies at the mutated positions in each
                              contig, as well as co-occurrence information
                              between pairs of mutated positions, will be
                              written. Each contig will be represented by two
                              "pickle" files within this directory, one
                              describing nucleotide frequencies at each
                              position and one describing nucleotide pair
                              frequencies at pairs of positions; each file's
                              name will be prefixed with the corresponding
                              contig name.  [required]
  --verbose / --no-verbose    Display extra details while running.  [default:
                              no-verbose]
  -h, --help                  Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
link nt @ 0.00s: Starting...
Input contig file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs_selected3.fasta
Input BAM file: /Poppy/mfedarko/sftests/tutorial-output/alignment/final.bam
Input BCF file: /Poppy/mfedarko/sftests/tutorial-output/p15-fdr1pct.bcf
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/link/nt
--------
link nt @ 0.00s: Loading and checking FASTA, BAM, and BCF files...
link nt @ 0.12s: The FASTA file describes 3 contig(s).
link nt @ 0.22s: All FASTA contig(s) are included in the BAM file (this BAM file has 78,793 reference(s)).
link nt @ 11.96s: All FASTA contig(s) are included in the BCF file (the header of this BCF file describes 468 contig(s)).
link nt @ 11.98s: The lengths of all contig(s) in the FASTA file match the corresponding lengths in the BAM and BCF files.
link nt @ 11.98s: So far, these files seem good.
--------
link nt @ 11.98s: Going through contigs and computing nucleotide (co-)occurrence information...
link nt @ 11.98s: On contig edge_1671 (2,153,394 bp) (1 / 3 contigs = 33.33%).
link nt @ 11.99s: Contig edge_1671 has 23,821 mutated position(s). Going through alignments to this contig...
link nt @ 899.59s: Found 256,509 read(s) spanning ≥ 1 mutated position in contig edge_1671. Computing (co-)occurrence info...
link nt @ 3,115.82s: Wrote out (co-)occurrence info for contig edge_1671.
link nt @ 3,115.82s: On contig edge_2358 (2,806,161 bp) (2 / 3 contigs = 66.67%).
link nt @ 3,115.82s: Contig edge_2358 has 17 mutated position(s). Going through alignments to this contig...
link nt @ 4,448.62s: Found 22,384 read(s) spanning ≥ 1 mutated position in contig edge_2358. Computing (co-)occurrence info...
link nt @ 4,451.40s: Wrote out (co-)occurrence info for contig edge_2358.
link nt @ 4,451.40s: On contig edge_6104 (1,289,244 bp) (3 / 3 contigs = 100.00%).
link nt @ 4,451.40s: Contig edge_6104 has 35 mutated position(s). Going through alignments to this contig...
link nt @ 5,292.62s: Found 8,670 read(s) spanning ≥ 1 mutated position in contig edge_6104. Computing (co-)occurrence info...
link nt @ 5,293.61s: Wrote out (co-)occurrence info for contig edge_6104.
link nt @ 5,293.61s: Done.
--------
link nt @ 5,293.66s: Done.

Usage: strainFlye link graph [OPTIONS]

  Convert (co-)occurrence information into a link graph structure.

  We create one link graph per contig. These graphs are undirected. A node (i,
  Ni) in the graph represents an allele: in other words, the occurrence of
  nucleotide Ni (one of {A, C, G, T}) at position i (1-indexed).

  What do edges between nodes represent? To answer that, let's define reads(i,
  Ni) as the number of reads at which Ni is aligned to i. Let's also define
  reads(i, j, Ni, Nj) as the number of reads at which Ni is aligned to i, and
  Nj is aligned to j.

  Given these definitions, we define the link weight between two alleles (i,
  Ni) and (j, Nj) as

                               reads(i, j, Ni, Nj)
   link(i, j, Ni, Nj) = -----------------------------------
                          max(reads(i, Ni), reads(j, Nj))

  We create an edge between nodes (i, Ni) and (j, Nj) if both of the following
  conditions hold:

  1. The number of reads spanning both i and j with any (non-degenerate)
  nucleotides aligned to either position is ≥ the --min-span parameter.

  2. link(i, j, Ni, Nj) > the --low-link parameter.

Options:
  -n, --nt-dir DIRECTORY          Directory produced by "strainFlye link nt"
                                  containing "pickle" files describing
                                  nucleotide (co-)occurrence information for
                                  certain contigs.  [required]
  --min-nt-count INTEGER RANGE    We will only create a node representing a
                                  given nucleotide at a mutated position if
                                  this nucleotide is supported by at least
                                  this many reads at this position.  [default:
                                  2; x>=1]
  --min-span INTEGER RANGE        One of the prerequisites for creating an
                                  edge between two nodes (representing
                                  nucleotides Ni, Nj at mutated positions i
                                  and j) is that the total number of reads
                                  spanning both i and j, across all
                                  combinations of (non-degenerate) nucleotides
                                  at both i and j, is at least this parameter.
                                  [default: 501; x>=1]
  --low-link FLOAT RANGE          The other prerequisite for creating an edge
                                  between two nodes (representing nucleotides
                                  Ni, Nj at mutated positions i and j) is that
                                  link(i, j, Ni, Nj) > this parameter. Note
                                  this check is not inclusive (i.e. the
                                  default of 0 means that, if reads(i, j, Ni,
                                  Nj) = 0, we will never create an edge
                                  between (i, Ni) and (j, Nj).  [default: 0;
                                  0<=x<1]
  -f, --output-format [dot|nx]    Output format used for each contig's link
                                  graph. "dot" will write out graphs in
                                  Graphviz' DOT language; "nx" will write out
                                  graphs to "pickle" files, in which each link
                                  graph is a NetworkX object.  [default: dot]
  --min-penwidth-clamp FLOAT RANGE
                                  Only used if --output-format is "dot". Edges
                                  will be assigned a "penwidth" attribute by
                                  interpolating their link weight from the
                                  range [0, 1] to the range [0, --max-
                                  penwidth]; penwidths less than --min-
                                  penwidth-clamp will then be clamped to
                                  --min-penwdith-clamp.  [default: 0.1; x>=0]
  --max-penwidth FLOAT RANGE      Only used if --output-format is "dot". See
                                  the description for --min-penwidth-clamp.
                                  [default: 5; x>=0]
  -o, --output-dir DIRECTORY      Directory to which graphs will be written.
                                  We'll write one graph file per contig; each
                                  file's name will be prefixed with the
                                  corresponding contig name.  [required]
  --verbose / --no-verbose        Display extra details while running.
                                  [default: no-verbose]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
link graph @ 0.00s: Starting...
Input nucleotide information directory: /Poppy/mfedarko/sftests/tutorial-output/link/nt
To make a node (i, Ni): reads(i, Ni) must be ≥ 2
To make an edge btwn (i, Ni) & (j, Nj): spanCount(i, j) must be ≥ 501
To make an edge btwn (i, Ni) & (j, Nj): link(i, j, Ni, Nj) must be > 0.00
Output graph file format: dot
In dot output, edges' penwidths will be scaled from link weights to [0, 5.0] and then clamped to a minimum of 0.1
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/link/graph
--------
link graph @ 0.00s: Going through (co-)occurrence information and creating link graphs...
link graph @ 25.65s: Creating a link graph for contig edge_1671...
link graph @ 25.73s: The link graph for contig edge_1671 has 49,167 node(s).
link graph @ 54.96s: The link graph for contig edge_1671 has 8,801,972 edge(s).
link graph @ 84.48s: Wrote out a link graph (format: "dot") for contig edge_1671.
link graph @ 86.67s: Creating a link graph for contig edge_2358...
link graph @ 86.67s: The link graph for contig edge_2358 has 41 node(s).
link graph @ 86.67s: The link graph for contig edge_2358 has 13 edge(s).
link graph @ 93.57s: Wrote out a link graph (format: "dot") for contig edge_2358.
link graph @ 93.57s: Creating a link graph for contig edge_6104...
link graph @ 93.57s: The link graph for contig edge_6104 has 80 node(s).
link graph @ 93.57s: The link graph for contig edge_6104 has 2,190 edge(s).
link graph @ 93.58s: Wrote out a link graph (format: "dot") for contig edge_6104.
link graph @ 93.58s: Done.
--------
link graph @ 93.58s: Done.

/home/mfedarko/miniconda3/envs/strainflye/bin/sfdp

Usage: strainFlye matrix [OPTIONS] COMMAND [ARGS]...

  [+] Create codon and amino acid mutation matrices.

Options:
  -h, --help  Show this message and exit.

Commands:
  count  Count 3-mers aligned to contigs' codons.
  fill   Call codon mutations from 3-mer counts and fill in matrices.

Usage: strainFlye matrix count [OPTIONS]

  Count 3-mers aligned to contigs' codons.

Options:
  -c, --contigs PATH          FASTA file of contigs in which we'll count
                              3-mers aligned to codons in predicted genes. All
                              contigs in this FASTA file should also be
                              contained in the BAM file; it's ok if the BAM
                              file contains contigs not in this FASTA file
                              (we'll ignore them). We'll also ignore contigs
                              for which no gene coordinates are given.
                              [required]
  -b, --bam PATH              Sorted and indexed BAM file representing an
                              alignment of reads to contigs.  [required]
  -g, --genes PATH            Generic Feature Format version 3 (GFF3) file
                              describing gene coordinates within contigs.
                              We'll only consider features with a "type" of
                              one of {"CDS", "SO:0000316"}.  [required]
  -o, --output-dir DIRECTORY  Directory to which 3-mer count information for
                              each contig will be written. Each contig will be
                              represented by one "pickle" file in this
                              directory named [contig]_3mers.pickle.
                              [required]
  --verbose / --no-verbose    Display extra details while running.  [default:
                              no-verbose]
  -h, --help                  Show this message and exit.

Usage: strainFlye matrix fill [OPTIONS]

  Call codon mutations from 3-mer counts and fill in matrices.

  You can call either p-mutations (based on codon frequencies) or r-mutations
  (based just on codon counts) using -p or -r. These methods work analogously
  to how they are described in "strainFlye call", with the addition that we do
  not call codon mutations for "unreasonable" codons where the most common
  aligned 3-mer (or, in the case of a tie, any of the most common aligned
  3-mers) at this codon is not exactly the same as the "reference" codon's
  3-mer.

  It's worth noting that the amino acid mutation matrix is derived directly
  from the codon mutation matrix: that is, as soon as we identify a mutation
  from one codon (C1) into another codon (C2), we also record a mutation in
  the amino acid matrix from translate(C1) into translate(C2). We could
  arguably do this in a different way (count the translations of each 3-mer
  aligned to a codon, then call amino acid mutations separately), but for now
  we just do things the one way.

Options:
  -c, --count-dir DIRECTORY       Directory produced by "strainFlye matrix
                                  count" containing "pickle" files describing
                                  the counts of 3-mers aligned to contigs'
                                  codons.  [required]
  -p INTEGER RANGE                If specified, we'll call codon p-mutations
                                  at this p threshold. This is scaled up by
                                  100 (i.e. "-p 50" corresponds to 50 / 100 =
                                  0.5%) in order to bypass floating-point
                                  precision issues. Mutually exclusive with
                                  -r.  [default: (nothing); 0<x<=5000]
  --min-alt-pos INTEGER RANGE     Only used if -p is specified. In order for
                                  us to call a p-mutation at a position, this
                                  position's alternate nucleotide must be
                                  supported by at least this many reads.
                                  [default: 2; x>=1]
  -r INTEGER RANGE                If specified, we'll call codon r-mutations
                                  at this r threshold. Mutually exclusive with
                                  -p.  [default: (nothing); x>=1]
  -f, --output-format [tsv|json]  Output format for matrix information.
                                  [default: tsv]
  -o, --output-dir DIRECTORY      Directory to which matrix information will
                                  be written. Within this directory, we'll
                                  create one subdirectory for each contig
                                  represented in --count-dir. Each
                                  subdirectory will contain a codon mutation
                                  matrix, amino acid mutation matrix, and
                                  counts of codons and amino acids in the
                                  "reference" contig sequence.  [required]
  --verbose / --no-verbose        Display extra details while running.
                                  [default: no-verbose]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
matrix fill @ 0.00s: Starting...
Input 3-mer count directory: /Poppy/mfedarko/sftests/tutorial-output/matrix/count
p: 50
Minimum number of supporting reads needed to call a p-mutation: 2
r: None
Output matrix information file format: tsv
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/matrix/fill
--------
matrix fill @ 0.00s: Performing codon p-mutation calling (p = 0.50%).
--------
matrix fill @ 0.00s: Going through aligned 3-mer counts and creating matrices...
matrix fill @ 1.26s: Creating matrices for contig edge_1671...
matrix fill @ 66.95s: Created an output directory for contig edge_1671.
matrix fill @ 68.38s: Creating matrices for contig edge_2358...
matrix fill @ 143.64s: Created an output directory for contig edge_2358.
matrix fill @ 144.56s: Creating matrices for contig edge_6104...
matrix fill @ 183.52s: Created an output directory for contig edge_6104.
matrix fill @ 183.52s: Done.
--------
matrix fill @ 183.64s: Done.

/home/mfedarko/miniconda3/envs/strainflye/lib/python3.7/site-packages/ipykernel_launcher.py:48: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.

Usage: strainFlye dynam covskew [OPTIONS]

  Compare coverage and GC skew within contigs.

  Bins
  ----
  We split each contig into bin(s) of a fixed number of positions (determined
  by the --bin-length parameter), starting from the leftmost end of the
  contig (which we define as position 1). We continue creating bins until
  we run out of positions: the rightmost bin will contain fewer positions
  than --bin-length if the contig length is not divisible by --bin-length.
  (And if the contig's length is ≤ --bin-length, we'll just create one
  bin for all positions in this contig.)

  Normalized coverages
  --------------------
  For each bin, we compute the median coverage of all positions in this bin.
  We then compute M, the median of these medians across all bins in a contig.
  We finally compute "normalized" coverages by dividing each bin's median
  coverage by M, and then clamping these normalized coverages to the range
  [1 - ϵ, 1 + ϵ]. (If M = 0, then we cannot do this normalization:
  in this case, we define each bin's normalized coverage as "NA".) The
  --norm-coverage-epsilon parameter can be used to adjust ϵ.

  Cumulative GC skews
  -------------------
  For each bin, we compute the skew (G - C) / (G + C), where G is the number
  of G nucleotides in this bin and C is the number of C nucleotides in this
  bin.  (If a bin has zero G or C nucleotides, then we define its skew as 0.)
  We make these skews cumulative by, starting from the second-from-the-left
  bin, defining this bin's cumulative skew as the sum of its skew and the
  cumulative skew of the bin exactly to the left of it. This approach is
  derived from (Grigoriev, 1998).

Options:
  -c, --contigs PATH              FASTA file of contigs for which we will
                                  compute coverage and skew information. All
                                  contigs in this FASTA file should also be
                                  contained in the BAM file; it's ok if the
                                  BAM file contains contigs not in this FASTA
                                  file (we'll ignore them).  [required]
  -b, --bam PATH                  Sorted and indexed BAM file representing an
                                  alignment of reads to contigs.  [required]
  -l, --bin-length INTEGER RANGE  Bin length (in bp).  [default: 10000; x>=1]
  -e, --norm-coverage-epsilon FLOAT RANGE
                                  We clamp all bins' normalized coverages to
                                  within this distance of 1.0. For example,
                                  the default of ϵ = 0.3 means we clamp to the
                                  range [0.7, 1.3].  [default: 0.3; x>0]
  -o, --output-dir DIRECTORY      Directory to which TSV files describing
                                  contigs' bins, coverages, and skews will be
                                  written. We'll write one TSV file per contig
                                  (named [contig]_covskew.tsv); each bin
                                  within a contig will be described in one row
                                  in its TSV file.  [required]
  --verbose / --no-verbose        Display extra details while running.
                                  [default: no-verbose]
  -h, --help                      Show this message and exit.

Using strainFlye version "0.1.0-dev".
--------
dynam covskew @ 0.00s: Starting...
Input contig file: /Poppy/mfedarko/sftests/tutorial-output/sheepgut_contigs_selected3.fasta
Input BAM file: /Poppy/mfedarko/sftests/tutorial-output/alignment/final.bam
Bin length: 10,000 bp
Normalized coverage epsilon: 0.3
Verbose?: Yes
Output directory: /Poppy/mfedarko/sftests/tutorial-output/covskew
--------
dynam covskew @ 0.00s: Loading and checking FASTA and BAM files...
dynam covskew @ 0.35s: The FASTA file describes 3 contig(s).
dynam covskew @ 0.35s: All of these are included in the BAM file (which has 78,793 reference(s)), with the same lengths.
--------
dynam covskew @ 0.35s: Going through contigs and computing coverage/skew information...
dynam covskew @ 0.35s: On contig edge_1671 (2,153,394 bp) (1 / 3 contigs = 33.33%).
dynam covskew @ 0.35s: Creating 215 bins of length 10,000 bp and 1 smaller bin of length 3,394 bp for contig edge_1671...
dynam covskew @ 231.35s: For contig edge_1671's bins' covs, pre-norm.: min = 4.0x; max = 1,606.0x; median (M) = 1,435.5x.
dynam covskew @ 231.35s: For contig edge_1671's cumulative bin skews: min = 0.0963; max = 4.9053.
dynam covskew @ 231.39s: On contig edge_2358 (2,806,161 bp) (2 / 3 contigs = 66.67%).
dynam covskew @ 231.39s: Creating 280 bins of length 10,000 bp and 1 smaller bin of length 6,161 bp for contig edge_2358...
dynam covskew @ 1,172.61s: For contig edge_2358's bins' covs, pre-norm.: min = 2,427.0x; max = 4,194.0x; median (M) = 2,940.0x.
dynam covskew @ 1,172.61s: For contig edge_2358's cumulative bin skews: min = -2.1523; max = 2.1512.
dynam covskew @ 1,172.67s: On contig edge_6104 (1,289,244 bp) (3 / 3 contigs = 100.00%).
dynam covskew @ 1,172.67s: Creating 128 bins of length 10,000 bp and 1 smaller bin of length 9,244 bp for contig edge_6104...
dynam covskew @ 1,888.12s: For contig edge_6104's bins' covs, pre-norm.: min = 1,455.0x; max = 5,104.5x; median (M) = 4,108.0x.
dynam covskew @ 1,888.12s: For contig edge_6104's cumulative bin skews: min = 0.0137; max = 7.6426.
dynam covskew @ 1,888.15s: Done.
--------
dynam covskew @ 1,888.19s: Done.

The normalized coverage in the bin with the minimum skew is 1.1831.
The normalized coverage in the bin with the maximum skew is 0.9282.

The PTR for CAMP is 1.1831 / 0.9282 = 1.2746.