The Minigraph-Cactus Pangenome Pipeline

July 27, 2026 · View on GitHub

Minigraph-Cactus is included in the Cactus Software Package and is suitable for aligning similar samples, such as those from the same species. See Progressive Cactus for aligning different species.

Please cite the Minigraph-Cactus paper when using Minigraph-Cactus.

Table of Contents

Pangenome Data

Some pangenomes constructed with Minigraph-Cactus, along with all material to reproduce, can be found here.

Quick-Start

The Mingraph-Cactus pipeline is run with cactus-pangenome. Its input is a text file with two columns, separated by whitespace, denoting sample name and fasta URL, respectively. All output is written to a directory specified with --outDir. All output filenames are prefixed with the value of --outName.

Begin by making a pangenome for tiny simulated primates:

cactus-pangenome ./js ./examples/evolverPrimates.txt --outDir primates-pg --outName primates-pg --reference simChimp --vcf --giraffe --gfa --gbz

If it worked properly, the input sequences should show up as contigs in the GFA:

zcat primates-pg/primates-pg.gfa.gz | grep '^W' | awk '{print \$1 "\t" \$2 "\t" \$3 "\t" \$4 "\t" \$5 "\t" \$6 }' | more
W	simChimp	0	simChimp.chr6	0	596350
W	simGorilla	0	simGorilla.chr6	0	599081
W	simHuman	0	simHuman.chr6	0	597871
W	simOrang	0	simOrang.chr6	0	591073

The reference path, which will be treated differently by some vg tools (for efficiency), will be identified in the header:

zcat primates-pg/primates-pg.gfa.gz | head -1
H	VN:Z:1.1	RS:Z:simChimp

The VCF will be based on the reference path (simChimp) and have a sample for each haplotype :

gzip -dc primates-pg/primates-pg.vcf.gz | grep CHROM -A 1
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	simGorilla	simHuman	simOrang
simChimp.chr6	15	>1>4	T	G	60.0	.	AC=1;AF=0.333333;AN=3;AT=>1>2>4,>1>3>4;NS=3;LV=0	GT	0	0	1

The three input files needed for giraffe are produced:

ls -hs primates-pg/primates-pg.d2.*
4.8M primates-pg/primates-pg.d2.dist  2.5M primates-pg/primates-pg.d2.gbz  6.1M primates-pg/primates-pg.d2.min

By default, --giraffe will produce frequency filtered indexes, with a default minimum coverage of 2 (hence the .d2). This means only nodes covered by two haplotypes will appear in the index. This helps vg giraffe performance considerably (though a version of Giraffe that no longer needs it is under development). The dataset here is too small for this to be useful. To index the clipped but unfiltered graph, use --giraffe clip or use --giraffe full to index the full, unclipped graph. See more detailed explanations below.

Introduction

Minigraph-Cactus uses minigraph to construct a pangenome graph of structural variation in a set of input assemblies. The assemblies are then mapped back to this graph using minigraph. These mappings are used as input to Cactus to construct a new graph that contains variants of all sizes, allowing the input assemblies to be encoded as embedded paths in the graph. The graph is output in pangenome graph formats such as vg and GFA, in addition to the usual HAL.

Pangenomes from Minigraph-Cactus are indexable for and ready for read mapping with vg Giraffe.

Unlike Progressive Cactus, Minigraph-Cactus does depend on a predetermined reference genome. This genome is guaranteed to be acyclic and unclipped in the final graph. Other genomes in the graph can still be used as reference coordinate systems, however. For example, we achieved great variant calling performance when projecting mapped reads on a CHM13-referenced graph onto GRCh38.

Interface

To run the Minigraph-Cactus pangenome pipeline, the basic format is:

cactus-pangenome <jobStorePath> <seqFile> --outDir <output directory> --outName <output file prefix> --reference <reference sample name>

The jobStorePath is where intermediate files, as well as job metadata, will be kept by Toil. It must be accessible to all worker systems.

The seqFile is a two-column mapping sample names to fasta paths (gzipped fastas are supported). The seqfile is the same as Progressive Cactus, except a tree is not specified.

Unlike for Progressive Cactus, Minigraph-Cactus ignores softmasking of the input sequences. You therefore do not need to mask the input. And like Progressive Cactus, hardmasking is not recommended.

Running on Clusters or Cloud

Please see the Progressive Cactus documentation for how to run Cactus on SLURM and AWS.

Sample Names

A naming convention must be followed for sample names where "." characters have special meaning: The "." character is used to specify haplotype, and should be avoided in sample names unless it is being used that way. For haploid samples, just don't use a ".". For diploid or polyploid samples, use the form SAMPLE.HAPLOTYPE where HAPLOTYPE is 1 or 2 for a diploid sample etc:

# Diploid sample:
HG002.1  ./HG002.paternal.fa
HG002.2  ./HG002.maternal.fa

# Haploid sample:
CHM13  ./chm13.fa

Contig Names

Ideally, the sequence/contig names in the input FASTA files should be as simple as possible. Ex:

>chr1
AAAA....
>chr2
CCCC....
etc

They do not need to be globally unique: the same contig name can be present in more than one input FASTA file. But the same name can only appear once per file: you cannot have a genome with two chr1s (you will get an error almost right away in this case). Some other things to know about contig names:

  • If you have a name of the form HG002#0#chr1, then cactus-pangenome will automatically ignore everything up to and including the last # character (leaving just chr1). The sample name and haplotype will be derived from the Sample Name in the seqfile (as described above).
  • vg can sometimes output paths with names like HG002#0#chr1#0. This will be interpreted as 0 and almost surely lead to trouble. You must manually correct this type of input before running cactus-pangenome.
  • Path subranges denoted by suffixes of the form :start-end (1-based, inclusive) are now supported (v2.8.1). The subrange information will be properly handled in the final output (VCF,GFA,GBZ etc).
  • : characters in contig names that are not part of subranges will be automatically replaced with _ in the output, as they can otherwise confuse the pipeline.
  • Anything after the first space in a contig name is ignored.

Reference Sample

The --reference option must be used to select a "reference" sample. This sample will:

  • Never be clipped.
  • Never be self-aligned ie its path in the output graph will be acyclic
  • Only visit nodes in their forward orientations
  • Be a "reference-sense" path in vg/gbz and will therefore be indexably for fast coordinate lookup
  • Be the basis for the output VCF and therefore won't appear as a sample in the VCF
  • Be used to divide the graph into chromosomes

It is therefore extremely important that the reference sample's assembly be chromosome scale. If there are many small contigs in the addition to chromosomes in the reference assembly, then please consider specifying the chromosomes with --refContigs. If you still want to keep the other contigs, add --otherContig chrOther (see explanation below).

Multiple Alignment Scoring Parameters

The --lastTrain option, added as an experimental prototype in v2.9.3, can be used to train substitution, match and gap scores from the input data using last-train. The model is inferred from an alignment between the (first) referene sample and the most diverged input genome from it in the minigraph construct phase, and these scores are then used in the cactus alignment phase to compute the base-level multiple alignment.

The default scoring parameters, found in the config XML file, are derived from the HOXD70 matrix used in lastz and are most appropriate for very diverged genomes. But for pangenomes these scores can lead to, for example, long runs of transitions that really should be gaps. In this case, using the data-derived scores should be beneficial.

You can also use the --scoresFile to input any .train file from last-train and use the scores contained in it to override cactus's defaults. Be careful: only fully symmetric scoring matrices are accepted.

Self-Alignment and the Collapse Option

The --collapse option, added as an experimental prototype in v2.9.1, can be used to incorporate self-alignments into the pangenome, including the reference sample. This will produce a more compact graph with, for example, tandem duplications being represented as cycles (like PGGB) rather than insertions (like minigraph). It also drops the invariant (see above) that the reference path by acyclic -- with this option, the reference path can contain cycles. The --collapse option is implemented by running minimap2 using options found in the minimapCollapseOptions parameter in the configuration XML to align each input contig with itself. These self alignments are fed into Cactus alongide the usual sequence-to-minigraph alignments.

Note: Haplotype indexes (--haplo) cannot always be built for collapsed graphs. If haplotype indexing fails when using --collapse, a warning will be printed in the log and the failing index will not appear in the output (but the program will not crash).

Multiple Reference Samples

The --reference option can accept multiple samples (separated by space). If multiple samples are specified beyond the first, they will be clipped as usual, but end up as "reference-sense" paths in the vg/gbz output. They can also be used as basis for VCF, and VCF files can be created based on them with the --vcfReference sample.

For example, for human data one might consider using --reference CHM13 GRCh38 GRCh37 --vcfReference CHM13 GRC3h8. This will make a graph referenced on CHM13, but will promote GRCh38 and GRCh38 to reference-sense paths so that they could be used, for example, to project BAMs on in vg giraffe. Two VCFs will be output, one based on CHM13 and one based on GRCh38.

Note: if you are using the step by step interface instead of cactus-pangenome, make sure to pass the same --reference options to each tool!

Running Minigraph on Chromosomes Independently

By default, minigraph construction and mapping are performed at the whole-genome level. The results are then split by (reference) chromosome so that Cactus (and some subsequent processing) can be run one chromosome at a time. The cactus-pangenome --mgSplit option was introduced to run the minigraph construction and mapping independently at the chromosome level instead. The intention is that by allowing the pipeline to choose a different ordering for minigraph construction for each chromosome, the accuracy of the output can be improved. This comes at the cost of extra computation (minigraph construction and mapping run twice). The step-by-step interface is likewise supported to support this and, using the small included example, would look like

cactus-minigraph js examples/evolverPrimates.txt ep.sv.gfa.gz --refOnly --reference simChimp 
cactus-graphmap js examples/evolverPrimates.txt ep.sv.gfa.gz ep.paf --reference simChimp --mgSplit --outputFasta ep.sv.fa.gz
cactus-graphmap-split js examples/evolverPrimates.txt ep.sv.gfa.gz ep.paf --reference simChimp --mgSplit --outDir out-split
cactus-minigraph js ./out-split/chromfile.txt out-construct --reference simChimp --batch
cactus-graphmap js ./out-construct/chromfile.mg.txt out-map --reference simChimp --batch
cactus-align js ./out-map/chromfile.gm.txt out-align --reference simChimp --outVG --pangenome --batch
cactus-graphmap-join js --vg out-align/*.vg --hal out-align/*.hal --sv-gfa out-construct/*.gfa.gz --reference simChimp --outDir out-join --outName ep --gb

Pipeline

The Minigraph-Cactus pipeline is run via the cactus-pangenome command. It consists of five stages which can also be run individually (below). cactus-pangenome writes output files into --outDir at the end of each stage. So different stages can be rerun with if necessary using the lower-level commands.

Before running large jobs, it is important to consider the following options:

  • --mgCores the number of cores for minigraph construction (default: all available)
  • --mgMemory the amount of memory for minigraph construction. The default estimate can be quite conservative (ie high), so if it is too high for your system, you can lower it with this option (default: estimate based on input size).
  • --mapCores the number of cores for each minigraph mapping job (default: up to 6)
  • --consCores the number of cores for each cactus-consolidated job (default: all available)
  • --consMemory the amount of memory for each cactus-consolidated job. By default, it is estimated from the data. When running on SLURM, Cactus now automatically clamps this (and every other) memory request down to the largest available node, so an over-large estimate will no longer leave the job unschedulable. You can still set --consMemory explicitly if an estimate comes out too low and the job runs out of memory, though --doubleMem usually handles that (this seems to be more of an issue for non-human data).
  • --indexCores the number of cores for each vg indexing job (default: all available - 1)
  • --indexMemory Like --consMemory above, you may want to set this when running on a cluster.
  • The various output options: --gbz, --gfa, --giraffe, --vcf, --odgi, etc. which are explained in detail below. If you forget to add one of these and are missing the corresponding output, you will need to rerun cactus-graphmap-join (or use vg to manually make the file yourself).

PLEASE NOTE While many Minigraph-Cactus parameters' default values were tuned on high-quality human assemblies from the HPRC where ample benchmarking data was available, we believe they will be suitable for other datasets and species, so long as the contigs can be mapped with minigraph. By default, small contigs are filtered out during chromosome assignment using more stringent thresholds. This might lead to a surprisingly low sensitivity on small, fragmented, diverse assemblies or difficult-to-assemble regions. Users wishing to keep these contigs in their graph can use the following option:

  • --permissiveContigFilter

The application and impact of this option is demonstrated in the explanation of the Yeast pangenome example below.

You can disable reference chromosome splitting using the cactus-pangenome --noSplit option. This will effectively remove the cactus-graphmap-split step from the pipeline, meaning that all chromosomes will be aligned in a single job. Contigs that map to mulitple chromosomes will not be filtered, potentially resulting in fewer components in the graph. This all comes at a potentially substantial increase in memory usage for all pipeline stpes beginning with Cactus.

Important The reference genome assembly must be chromosome scale. If your reference assembly also consists of many small fragments (ex GRCh38) then you can use the --refContigs option to specify the chromosomes. Ex for GRCh38 --refContigs $(for i in $(seq 22); do printf "chr$i "; done ; echo "chrX chrY chrM"). If you want to include the remaining reference contig fragments in your graph, add the --otherContig chrOther option. If you do not specify --refContigs, they will be determined automatically and all small contigs will be included. For example, if your reference has the following contigs and lengths

chr1 1000000
chr2 900000
chr3 500000
chrM 1000
random_1 400000
random_2 5000
random_3 5000

Then the pipeline will automatically determine, using the names and sizes (see the section of the configuration XML for details) that the reference contigs are chr1, chr2, chr3, chrM and random_1. It will make a graph file for each of these, and a single chrOther graph file for random_2 and random_3. All contigs will end up in the final whole-genome indexes. This logic only affects the chromosome-scale output and workload distribution. If you do not want the random contigs in the graph, you would need to specify --refContigs chr1 chr2 chr3 chrM. By default, refContigs are limited to 128 in number, contigs that are at least 10% as long as the longest contig or contigs whose name begins with chr followed by up to 3 characters (any case).

Also Important We do not yet automatically support the alternate loci from GRCh38, ex the various HLA contigs. They must be excluded from the input fasta file to get sane results. They can be included in the graph by providing a separate sample / fasta pair in the input for each contig. Please here for an example of how to do so.

The individual parts of the pipeline can be run independently using the following commands. See also the "step-by-step" yeast example below.

  1. cactus-minigraph <jobStore> <seqFile> <outputGFA> --reference: Construct a minigraph in GFA format (may be gzipped) from a set of FASTA files (may also be gzipped). This is a very thin wrapper over minigraph -cxggs. The reference is added first and the remainder of samples are added in decreasing order of decreasing mash distance to the reference (see the minigraphSortInput parameter in the XML config to change or disable this). Use the --mapCores option to specify the number of cores.

  2. cactus-graphmap <jobStore> <seqFile> <inputGFA> <outputPAF> --reference: Map each input assembly back to the graph using minigraph. The number of cores for each mapping job can be set with --mapCores.

  3. (Optional) cactus-graphmap-split <jobStore> <seqFile> <inputGFA> <inputPAF> --reference --outDir: Split the input assemblies and PAF into chromosomes using the rGFA tags in the GFA. Doing so reduces the memory requirements in the following steps. It assigns each contig to a single chromosome according to the alignment in the input PAF, so all inter-chromosomal events will be filtered out. Contigs that can't be assigned to a chromosome are deemed "ambiguous" and not considered in later steps.

  4. cactus-align <jobStore> <seqFile> <inputPAF> <outHal> --reference --pangenome --outVG --maxLen: Compute the Cactus multiple genome alignment from the assembly-to-graph minigraph mappings. The --maxLen parameter specifies the maximum gap between minigraph mappings that Cactus will attempt to fill at once, and is recommended to be set to 10000. If cactus-graphmap-split was used, the cactus-align --batch interface should be used instead (on the chromfile.txt output rather than seqfile.txt) (see examples below).

  5. cactus-graphmap-join <jobStore> --vg --outDir --outName --reference: Produce the final graph and indexes. This should be run whether or not cactus-graphmap-split was used.

Clipping, Filtering and Indexing

cactus-pangenome (options also available in cactus-graphmap-join) normalizes, clips and filters the graph in addition to producing some useful indexes. It can produce up to three graphs (now in a single invocation) in addition to the direct minigraph output, and a variety of indexes for any combination of them. The different graphs can be distinguished by their filenames. Suppose the tool as run with --outName yeast, then you may have these files in the output:

  • yeast.sv.gfa.gz: This graph is output by minigraph. It contains SVs only, and doesn't have embedded paths for the input sequences.
  • yeast.full.gfa.gz: This is the full minigraph-cactus graph. It is normalized, but no sequence is removed. It and its indexes will have .full in their filenames.
  • yeast.gfa.gz. This is the default or clip graph. Stretches of sequence >10kb that were not aligned to the underlying SV/minigraph are removed. "Dangling" nodes (ie that don't have an edge on each side) that aren't on the reference path are also removed, so that each chromosome only has two tips in the graph.
  • yeast.d2.gfa.gz: This filter graph is made by removing nodes covered by fewer than 2 haplotypes (this value can be changed using the --filter option) from the clip graph. Note in newer versions of vg, you can usually get away without allele frequency filtering by way of haplotype sampling (using the --haplo option (without --giraffe)) to make an index for this).

The clip graph is a subgraph of the full graph and the filter graph is a subgraph of the clip graph. Put another way, any node in the filter graph exists with the exact same ID and sequence in the clip graph, etc.

The different graphs have different uses. For instance, the current version of vg giraffe performs best on the filtered graph (this will hopefully be soon remedied in an update to vg). For the HPRC v1.0 graph and paper, we used d9. When you use --giraffe, it will make the giraffe indexes on the filtered graph by default. But you can override this behaviour to produces the indexes for any of the graphs by passing in any combination of [full, clip and filter] to the --giraffe options. For example:

--giraffe: Make the giraffe indexes for the filtered graph (default choice).

--giraffe clip: Make the giraffe indexes for the clipped graph.

--giraffe clip filter: Make the giraffe indexes for both the clipped and filtered graph.

Mapping long reads with giraffe requires a different type of minimizer index, which can be created with the --lrGiraffe option. --giraffe and --lrGiraffe can be used together or separately. They create the same .dist index which can be used for both short and long read mapping.

The same type of interface applies to all the output specification options: --vcf, --gbz, --gfa, --xg, --giraffe, --chrom-vg, --odgi, --chrom-og, --viz, --draw. They can all be used without arguments to apply to the default graph (generally the clip graph for everything except --giraffe which defaults to the filter graph, and anything odgi-related which defaults to full), or with any combination of full, clip and filter to be applied to different graphs.

Note that by default, only GFA is output, so the above options need to be used to toggle on any other output types.

Different clipping and filtering thresholds can be specified using the --clip and --filter options, respectively. For larger graphs, you probably want to use --filter N where N represents about 10% of the haplotypes. It is indeed a shame to remove rarer variants before mapping, but is a necessity to get the best performance out of (the current version) of vg giraffe.

Re-indexing with Pre-processed Graphs

If you have previously run cactus-pangenome or cactus-graphmap-join with --chrom-vg and want to regenerate indexes without repeating the clipping and filtering steps, you can use the bypass input options: --vgFull, --vgClip, and/or --vgFilter (available only in cactus-graphmap-join). Each accepts the per-chromosome VG files from a previous --chrom-vg full, --chrom-vg clip, or --chrom-vg filter output respectively, and skips directly to index building. These are incompatible with --vg but can be combined with each other. For example:

cactus-graphmap-join <jobstore> --vgClip yeast.chroms/chr*.vg --outDir reindex --outName yeast --reference S288C --gbz --giraffe

This will produce GBZ and Giraffe indexes from the pre-clipped chromosome VGs without re-running any graph processing.

For --vgFilter, the filter threshold is inferred from the .dX.vg filename pattern (e.g., chr1.d2.vg implies a filter of 2). All filter VG files must use the same threshold.

Note: per-chromosome output options (--chrom-vg, --chrom-og, --viz, --draw) cannot be used with bypass options, as you already have those files from the previous run. Also, bypass options are not compatible with graphs that were originally built with --collapse.

VCF Output

The --vcf option runs vg deconstruct to represent the graph as sites of variation along a reference. A single run can write several VCFs, all prefixed with --outName:

SuffixEnabled byWhat it is
.raw.vcf.gz--vcfThe vg deconstruct output, with sample ploidies padded so they agree across chromosomes. Nested sites are kept, marked with the LV (nesting level) and PS (parent snarl) tags.
.vcf.gz--vcf, --vcfbub non-zeroThe raw VCF put through vcfbub, which drops sites nested inside a larger site and pops open any site whose reference allele is longer than 100kb. Use this one unless you want to handle nesting yourself.
.wave.vcf.gz--vcfwavevcfwave realigns each alt allele against the reference and breaks it into primitive SNPs and indels. Built from the raw VCF, not from .vcf.gz: it runs its own vcfbub filtering on the longest allele rather than the reference allele, so it pops more sites and drops INFO/AT.
.L<NN>.vcf.gz, .L<NN>.raw.vcf.gz--vcfLAllele-clustered copies of the two VCFs above, written in addition to them (see below).
.gref.*--grefThe same set over again, deconstructed against the synthetic graph reference paths instead of the reference sample. Experimental.

Those suffixes attach to a base name that records which graph and which reference the VCF came from:

  • the clip graph contributes nothing: <outName>.vcf.gz. A bare --vcf selects it, except under --clip 0, where the default becomes full
  • --vcf full contributes full: <outName>.full.vcf.gz
  • --vcf filter contributes d<N>, where N is the --filter threshold: <outName>.d2.vcf.gz
  • a --vcfReference SAMPLE that is not the first --reference contributes its sample name: <outName>.SAMPLE.vcf.gz. VCFs built on the first reference carry no sample component

These combine, so --reference CHM13 GRCh38 --vcfReference GRCh38 --vcf filter --vcfwave writes <outName>.GRCh38.d2.wave.vcf.gz.

Switch off vcfbub with --vcfbub 0 or specify a different threshold with --vcfbub N.

The bcftoolsNorm config attribute runs bcftools norm -f over the vcfbub VCF, left-aligning and normalizing indels. It ships off (bcftoolsNorm="0"), since left-shifting can leave overlapping variants behind; set it to "1" in the config XML to turn it on. The wave VCF is normalized by default instead, under the separate vcfwaveNorm attribute.

Also new in v2.8.2, you can use the --vcfwave option to create a version of the VCF(s) that has been normalized with vcfwave. vcfwave can take a while to run on larger graphs, but it can do a good job of smoothing out small variants in complex regions. vcfwave is included in the Cactus docker images but not the the Cactus binary release. So you must either run Cactus from inside Docker, or run outside docker with the --binariesMode docker option. If you cannot run docker and still want to use --vcfwave, you must build it yourself and make sure its on your PATH before running cactus. You can set the number of cores used for each vcfwave job with --vcfwaveCores, and the memory with --vcfwaveMemory.

The experimental --vcfL option passes -L to vg deconstruct, merging alt alleles whose traversals are at least that similar (length-weighted Jaccard) into a single allele, which collapses near-identical alleles in complex bubbles. It is lossy: a sample merged into another allele is genotyped as that allele — possibly the reference allele — and the difference survives only in the TS/TL FORMAT fields. So it never replaces anything. The ordinary VCFs are still written, and clustered copies appear beside them tagged with L<NN>, e.g. --vcfL 0.95 gives <outName>.L95.vcf.gz. No clustered .wave.vcf.gz is produced, since vcfwave exists to undo precisely the allele merging that -L does.

Whenever normalization runs — so the wave VCF always, and the vcfbub VCF only if you turn bcftoolsNorm on — the VCF then goes through merge_duplicates.py (from the collapse-bubble submodule), which combines records left at the same position by left-alignment, and finally vcffixup to recompute AC/AF/AN. Set mergeDuplicatesOptions to "0" in the config to skip the merge step.

Graph Reference Paths (--gref)

Experimental. Non-reference sequence reaches a reference-based VCF only as ALT alleles, so variation inside an insertion or a non-reference SV allele has no reference position of its own: it appears in .raw.vcf.gz as nested LV>0 records, which vcfbub then drops from .vcf.gz. --gref gives that sequence coordinates so it can be addressed directly.

vg paths -u computes a graph reference path cover: it finds the parts of the graph the reference does not walk, and promotes fragments of the haplotype paths covering them into synthetic reference paths. These form a new sample, gref_<reference>, holding a copy of the reference paths plus fragments suffixed _<N>_alt (gref_GRCh38#0#chr1_1_alt). Deconstructing against that sample puts the non-reference material in the VCF.

cactus-pangenome ./js ./seqfile.txt --outDir pg --outName pg --reference GRCh38 --vcf --gref

--gref takes the graph to build from: full, clip or filter, defaulting to clip. --minGrefLen sets the minimum fragment length (default 50, or the minGrefLen config attribute). Only the first --reference sample is used, as its paths are the ones guaranteed acyclic. --gref clip requires clipping to be enabled, --gref filter requires filtering to be enabled, and neither can be combined with --collapse.

FileWhat it is
<outName>.gref.gbz, .gref.gfa.gzThe source graph with the gref_<reference> sample added.
<outName>.gref.vcf.gz etc.The VCF set from the table above (.raw, .wave, .L<NN>), deconstructed against gref_<reference> rather than the reference sample. Produced even without --vcf; .gref.wave.vcf.gz still requires --vcfwave, which itself requires --vcf.
<outName>.gref.gref-segs.tsv.gzWhere every synthetic fragment came from (see below).

The gref graph adds paths only, no nodes or edges, so it is topologically identical to the base graph and gets no topology indexes of its own. Use the base graph's <outName>.snarls, .dist and .hapl with gref.gbz. .gref.snarls and .gref.hapl are not produced.

The segment table maps a gref coordinate back onto an assembly one. It is headerless, tab-separated and BED-like, one row per fragment: the source haplotype path and interval, the new gref path name, then the reference path the fragment hangs off and the span of the enclosing snarl there. That span includes the flanking nodes, so the fragment sits strictly inside it rather than filling it.

SAMP2#0#chr1    200    700    gref_REF#0#chr1_1_alt    REF#0#chr1    0    901

gref_REF#0#chr1_1_alt is SAMP2#0#chr1:200-700, placed within REF#0#chr1:0-901. Fragments that could not be placed against the reference get . and 0 0 in the last three columns.

The gref VCF has records in two coordinate systems: on the reference contigs and on the gref contigs. The same sample sequence appears in both, as an ALT allele at a reference-contig site and as the reference for a gref contig's own records. RC/RS/RD join the two: on a gref-contig record they give the reference-contig site containing it. RC=chr22 RS=15470047 on a chr22_2_alt record means that fragment is the insertion called at chr22:15470047. RC/RS/RD survive vcfbub and vcfwave; INFO/AT, which would give the specific ALT allele, is stripped by the wave path.

Do not pool reference-contig and gref-contig records when computing allele frequencies or site counts, or that sequence is counted twice. In .wave.vcf.gz, match RS with a tolerance: realignment can shift the reference-contig record a few bases from where RS points.

vg counts LV within each reference contig rather than as absolute depth in the snarl tree, so gref-contig sites reach LV=0 and the whole VCF goes through a single vcfbub --max-level 0 pass.

Haplotype Sampling Instead of Filtering (NEW)

The .dX graphs created with --filter were necessary for read mapping, but now vg supports dynamic haplotype subsampling (ie personalized pangenomes) and, in most cases, filtering is no longer necessary. In order to use haplotype sampling, run cactus-pangenome / cactus-graphmap-join with the --haplo option (and do not use --giraffe). This will create the giraffe indexes for the special .hapl haplotype index which (with the .gbz`` is all you need to run vg giraffe` using the current best practices.

While this process will give better mapping performance than using the --filter graphs there is one downside:

  • Read mapping will now require an invocation of kmc to compute a kmer index (see links below). While this adds complexity, it does not seriously affect runtime (the time used making the kmer index and doing the subsampling is balanced out by faster mappings times).

Further reading:

Output

  • hal: Cactus's native alignment format can be used to convert to MAF, build assembly hubs, run liftover and comparative annotation.
  • gfa: A standard text-based graph format. Minigraph-Cactus uses GFA 1.1 as it represents haplotypes as Walks. You can use vg convert -gfW to convert from GFA 1.1 to 1.0 and vg convert -gf to convert from 1.0 to 1.1.
  • vcf: A standard text-based format that represents a pangenome graph as sites of variation along a reference. VCFs exported from the graph are nested, and by default vcfbub is used to flatten them.
  • vg: vg's native packed-graph format, can be read and written by vg but does not scale well with the number of paths.
  • gbz: A read-only format that scales extremely efficiently with the number of paths. Readable by vg tools and required for giraffe.
  • xg : A read-only format that is indexed on all paths (not just reference paths like .gbz). As a consequence, it requires more memory than .gbz. Readable by vg tools, and ideal for vg find and vg chunk.
  • snarls: The start and end nodes of the bubbles in the graph, as well as their nesting relationships. Used by some vg tools like call and deconstruct.
  • dist: Snarl distance index required for vg giraffe.
  • min: Minimizer index required for vg giraffe.
  • hapl : Haplotype sampling index. Created with --haplo (new alternative to --filter) and used for new best practice vg giraffe pipeline.
  • stats/: Tables describing the graph and a complete accounting of input sequence that did not make it in. See Statistics and Clipping Reports.
  • og: odgi's native format, can be read and written by odgi. Very useful for visualization.
  • snarl-stats.tsv.gz: Table with one row per snarl (bubble) in the graph, sorted in decreasing order of the distance they span on the (first) reference sample. See the table header for a description of each column, and note that snarls can be nested in eachother. This table is produced with the --snarlStats option.

Statistics and Clipping Reports

Every run writes a <outName>.stats/ directory describing what went into the graph and what did not, plus two files at the top level:

<outName>.WARNING                    only exists if something looks wrong -- see below
<outName>.input-contig-sizes.tsv.gz  input contig lengths, as cactus saw them
<outName>.stats/
    clipped-by-genome.tsv               one row per input genome and reason
    clipped-by-reference-contig.tsv     one row per genome per reference contig
    clipped-by-input-contig.tsv.gz      one row per input contig
    clipped.beds.tar.gz                 the intervals themselves, for the clipped graph
    clipped.full.beds.tar.gz            ...for the full graph
    clipped.d2.beds.tar.gz              ...for the frequency-filtered graph
    graph-stats.tsv                     one row per reference contig: nodes, edges, length
    path-stats.tsv.gz                   one row per path in the clipped graph
    refgaps.bed.gz                      reference regions no other assembly aligns to
What "clipped" means here

Any input base that is not in the output graph, whatever removed it. That is broader than the --clip phase alone, and the reason column says which stage was responsible:

reasonmeaning
ambiguouscontig binned to _AMBIGUOUS_ during chromosome splitting
unassignedcontig that splitting made no decision about at all
no_chromosome_graphcontig assigned to a reference contig no graph was built for (e.g. one left out of --refContigs)
unalignedcontig reached its chromosome graph, but some of its sequence is not in the full graph
clipremoved by the clip phase (--clip)
filterremoved by the allele-frequency filter (--filter)

The first four mean the sequence is absent from every graph. Only phases that actually ran are reported: with the default options there is no frequency-filtered graph, so there are no filter rows and no .d2 files. The report never causes a phase to run.

The three tables

They are the same accounting at three zoom levels, and all reconcile with each other. clipped-by-genome.tsv is the headline:

genome        reason      intervals  bp        pct_of_input
UWOPS034614   ambiguous   4          3255663   27.714
UWOPS034614   clip        15         389645    3.317
UWOPS034614   TOTAL       25390      4067475   34.625

Every reason the run could measure gets a row for every genome, so a 0 means measured-and-nothing-lost, while a reason with no rows at all means that phase never ran.

The other two say how much of the input survived into each graph, rather than what was lost, so the columns decrease left to right and the loss at any stage is the drop between adjacent _bp columns:

#ref_chrom  genome  contigs  input_bp  full_bp  full_frags  clip_bp  clip_frags  filter_bp  filter_frags
chrI        SK1     1        228861    228861   1           214802   1           191745     869

_frags is how many contiguous pieces the sequence is broken into in that graph, so the example above shows clipping removing 14,059 bp while leaving one piece, then the frequency filter removing another 23,057 and shattering it into 869. A ref_chrom of _NONE_ means the contig reached no chromosome graph at all.

clipped-by-input-contig.tsv.gz is the same thing one level down, one row per input contig. It is most useful for fragmented assemblies, where a genome contributes many contigs to one reference contig; where each contig maps to its own reference contig it adds little over the rollup.

The BED archives

Each clipped*.beds.tar.gz unpacks to a directory of one plain 4-column BED per input genome, in input contig coordinates, usable with bedtools as-is:

chrVI      0        290867    unassigned
chrI       0        17071     clip

They are cumulative, since each graph is built from the previous one: everything missing from the full graph is also missing from the clipped graph. So clipped.beds.tar.gz is what is not in the clipped graph -- the default output, and what most people want.

The accounting is measured, not estimated: path coverage is read back out of each graph that was built and subtracted from the input contig lengths, so for every genome bases in the graph + bases in that graph's BED == input bases, exactly.

Graph statistics

graph-stats.tsv and path-stats.tsv.gz describe the clipped graph rather than what was clipped from it. graph-stats.tsv gives the size and complexity of each chromosome graph, where length is the total sequence in the graph -- much less than the sum of the inputs, because shared sequence is stored once. path-stats.tsv.gz has one row per path; a [start-end] suffix on a path name means it is a fragment of a longer input contig.

Reference gaps

refgaps.bed.gz is different in kind: the reference is never clipped, so it loses nothing. These are instead regions of the reference that no other assembly aligns to (runs of at least refGapMinLength, 10kb by default). This sequence is present in the graph and is never counted as clipped.

<outName>.WARNING

A run can finish successfully and still have gone badly wrong. If the accounting shows one of the shapes that usually means something is off, cactus writes <outName>.WARNING -- the file only exists when there is something to say, so its absence is the all-clear. It looks for a genome with a large fraction of itself in no graph at all, whole input contigs absent from every graph, one genome losing far more than the rest of the panel, a reference little of the panel aligns to, and an accounting that fails to add up.

It deliberately stays quiet on losses that are uniform across the panel, since those reflect how the run was configured rather than a problem: dropping whole chromosomes with --refContigs and removing rare alleles with --filter both do that. Each note is a heuristic, so check it against your data -- a genuinely rearranged assembly will trip the first one, correctly.

Running cactus-graphmap-join on its own

Pass --inputContigSizes <outName>.input-contig-sizes.tsv.gz from the cactus-pangenome run that produced the graphs. It is the only record of how much input sequence there was, so without it no clipping report is written -- only refgaps.bed.gz and the graph statistics, which do not depend on it. Measuring the graphs against each other instead would produce the same filenames and column names meaning something different, which is worse than producing nothing.

Node Chopping

As of v2.9.1, all output graphs will have node IDs of at most 1024bp. This is because the gbz and dist indexes require this (they use 10bits for node offsets), and as a result so do an increasing number of vg tools. There is also a major benefit from the simplicity of having all output files sharing the same ID space.

If you must have longer node lengths, perhaps for viewing very large simple graphs in BandageNG, then you can use the --unchopped-gfa option to produce additional output graph(s), with .unchopped.gfa.gz suffix. The node IDs in these graphs will be incompatible with all other output, including VCF.

If you want to revert to the previous logic of only the vg giraffe-specific files having chopped IDs, you can do so by setting <graphmap_join maxNodeLength="-1"> in the configuration XML. If you do this, then you can still use vg gbwt -Z graph.gbz --translation mapping.tsv to recover the translation between the IDs from the GBZ, or -O in vg call/deconstruct to toggle which IDs to use from the GBZ.

Patching Assemblies (cactus-panpatch)

panpatch uses a pangenome graph to patch an assembly into telomere-to-telomere chromosomes — filling N-gaps, scaffolding contigs, and (with --requireTelomeres) completing telomeres. cactus-panpatch is a Toil workflow that builds the required graphs and runs panpatch all in a single script and, with the --batch option, can process many samples at once.

cactus-panpatch takes a seqfile listing assemblies of the same sample. The target (to be patched) assembly comes first, with donor assemblies following in priority order. For example:

PAN028-verkko.1   PAN028.hap1.verkko.fa      # <- the assembly to patch (chromosome-scale)
PAN028-verkko.2   PAN028.hap2.verkko.fa
PAN028-hifiasm.1  PAN028.hifiasm.hap1.fa     # <- donor
PAN028-hifiasm.2  PAN028.hifiasm.hap2.fa

cactus-panpatch ./js seqfile.txt --outDir patched

By default patching is reference-free: because Minigraph-Cactus references must be haploid, each target haplotype is patched in its own pangenome (so a diploid sample costs N times a single run), and its output is named for it: patched/PAN028-verkko.hap1.fa.gz, patched/PAN028-verkko.hap2.fa.gz.

If the target assembly is not chromosome scale (you want the patching to scaffold across its contigs), you can include a chromosome-scale reference genome (ex hs1) as the first line of the seqfile, and then flag it as a reference using the --reference option (it must be haploid).

Patching many samples

--batch patches many samples in one workflow; the argument is a chromfile with one <name> <seqfile> line per sample (as with cactus-align --batch), and works with or without --reference:

PAN028   PAN028.seqfile.txt
PAN027   PAN027.seqfile.txt

cactus-panpatch ./js chromfile.txt --outDir patched --batch

Cluster example

cactus-panpatch ./js chromfile.txt --outDir patched --batch --requireTelomeres \
    --batchSystem slurm --slurmTime 10:00:00 --doubleMem true --maxMemory 1.5T --retryCount 10 \
    --consCores 16 --mgCores 32 --mapCores 16

Please cite the panpatch paper when using cactus-panpatch.

Advanced Configuration

See here

Visualization

Sequence Tube Maps

Sequence Tube Maps let you use your web browser to display pangenome graphs and their alignments. You should be able to load the gbz output from minigraph-cactus directly into a tube map. Any alignments (in .gam) format from vg giraffe can also be viewed, but they will need to be sorted first (a prepare script is provided). Note that while you should be able to load your entire graph into the tube map, you can only use it to display relatively small regions at a time.

This is the tool that was used to make most of Figure 1 in the Minigraph-Cactus paper.

Please cite Sequence Tube Maps for images you create with it!

Bandage-NG

Bandage-NG is an interactive viewer for graphs in GFA format. It is not included in Cactus, but it is very easy to install on Mac and Linux.

Depending on the size of the data set, you may need to extract a subregion from the graph in order to display it with Bandage-NG, which generally scales to about 10 or so mb of minigraph-cactus graphs. You can do this with vg chunk on the output .gbz file specifying -O gfa to set the output format to GFA.

Please see the GRCh38 Alts Graph section below for how to use vg chunk and Bandage-NG to view regions of interest in a human pangenome.

Please cite Bandage-NG for images your create with it!

ODGI

odgi is a pangenome analysis toolkit and format. As of version v2.6.1, Minigraph-Cactus includes the odgi binary and supports output in odgi format via the following cactus-pangenome / cactus-graphmap-join options:

  • --odgi : Output the graph to odgi (.og) format. Valid options are full and/or clip, with the default being full if none are specified (there is no reason to output the filter graphs to odgi, and the tiny path fragments can potentially make conversion very slow).
  • --chrom-og: Output each graph chromosome in odgi (.og) format. This is recommended if you want to run odgi yourself to do any visualization, since you will generally want to deal with one chromosome at a time. As above, valid options are full and clip, defaulting to full if none specified.

Unlike some of the related options such as --gbz, --xg, --chrom-vg, --gfa, etc, the odgi options above default to working on the full graph as opposed the the clip graphs. As such, the output (by default) will contain the unaligned chromosomes. The rationale is that these unaligned sequences do not seem to hinder visualization, whereas the path fragments that arise from clipping can bog odgi down a little bit (ex: chr1 from the clipped 10-chicken pangenome takes several hours to convert to odgi, but the full graph is fine). Also, odgi's coordinate system does not support path fragments, so its extraction tools etc. will only work properly on the full graphs. This is still a work in progress, so you can use, say, --chrom-og clip full to experiment with both.

ODGI Viz

odgi's 1-dimensional visualization is a very efficient way to view pangenome graphs at the chromosome scale. You can use odgi on the --chrom-og output described following the instructions from odgi's manual or use the --viz option in cactus-pangenome / cactus-graphmap-join to automatically produce a .png file for each chromosome in the graph.

Here is an example of using --viz in the yeast example below.

cactus-pangenome ./js ./examples/yeastPangenome.txt --reference S288C --outDir yeast-pg --outName yeast-pg --gbz --viz

The output is found here:

ls -h yeast-pg/yeast-pg.viz
chrI.full.viz.png    chrIV.full.viz.png  chrVI.full.viz.png    chrX.full.viz.png    chrXIII.full.viz.png  chrXVI.full.viz.png
chrII.full.viz.png   chrIX.full.viz.png  chrVII.full.viz.png   chrXI.full.viz.png   chrXIV.full.viz.png
chrIII.full.viz.png  chrV.full.viz.png   chrVIII.full.viz.png  chrXII.full.viz.png  chrXV.full.viz.png

And the first two chromosomes look like

ODGI Draw

odgi's 2-dimensional "draw" visualization resembles Bandage-NGs, but can be done at a larger scale. You can produce chromosome images using this functionality with the --draw option in cactus-pangenome / cactus-graphmap-join.

WARNING This option is very experimental -- I do not yet know the best set of parameters (and have not had much time to look) for drawing minigraph-cactus graphs. Please use with caution, especially on larger graphs, since any major issue with drawing may prevent you from obtaining the various other indexes. For large graphs, consider using the --chrom-og option to produce the odgi graphs for each chromosome, and then running draw yourself following the ODGI documentation. Please report back any discoveries here!

For example, we can add the 2D layouts to those created above using the following command (note the substantially longer running time):

cactus-graphmap-join ./js --vg ./yeast-pg/chrom-alignments/*.vg  --reference S288C --outDir yeast-pg --outName yeast-pg --draw

The output is written to the same place as --viz (see above):

ls -h yeast-pg/yeast-pg.viz/*.draw.png
yeast-pg/yeast-pg.viz/chrI.full.draw.png    yeast-pg/yeast-pg.viz/chrVI.full.draw.png    yeast-pg/yeast-pg.viz/chrXIII.full.draw.png
yeast-pg/yeast-pg.viz/chrII.full.draw.png   yeast-pg/yeast-pg.viz/chrVII.full.draw.png   yeast-pg/yeast-pg.viz/chrXIV.full.draw.png
yeast-pg/yeast-pg.viz/chrIII.full.draw.png  yeast-pg/yeast-pg.viz/chrVIII.full.draw.png  yeast-pg/yeast-pg.viz/chrXV.full.draw.png
yeast-pg/yeast-pg.viz/chrIV.full.draw.png   yeast-pg/yeast-pg.viz/chrX.full.draw.png     yeast-pg/yeast-pg.viz/chrXVI.full.draw.png
yeast-pg/yeast-pg.viz/chrIX.full.draw.png   yeast-pg/yeast-pg.viz/chrXI.full.draw.png
yeast-pg/yeast-pg.viz/chrV.full.draw.png    yeast-pg/yeast-pg.viz/chrXII.full.draw.png

And the first chromosome looks like

Please cite ODGI for any images and analysis you use it for.

VG

vg has two visualization tools, vg view -d and vg viz.

vg view -d outputs graphs in GraphViz "dot" format, and can be converted to images with dot -Tpng. This command only works on very small subgraphs (though you can increase its scope by leaving out bases with -S). As such, it is best used with vg chunk to extract subgraphs.

An example for yeast, may look like

vg chunk -x yeast-pg/yeast-pg.gbz -S yeast-pg/yeast-pg.snarls -p S288C#0#chrI:100400-100500 | vg view -pd - | dot -Tpng > yeast-pg-chunk-view.png

And the result would be

We can pull out a larger region with vg viz (note that vg viz can also display coverage output from vg pack which is not demonstrated here)

vg chunk -x yeast-pg/yeast-pg.gbz -S yeast-pg/yeast-pg.snarls -p S288C#0#chrI:100400-101400 | vg viz -x - -o yeast-pg-chunk-viz.png

Please cite vg when using these visualizations.

Pangenome Statistics (panacus)

panacus computes pangenome coverage and growth statistics: how much sequence is shared across the samples, how the pangenome grows as samples are added, and how big its core is. Passing --panacus to cactus-graphmap-join (or cactus-pangenome) runs panacus on the resulting whole-genome graph and writes both the source table data and a self-contained interactive HTML report of the plots.

cactus-pangenome ./js ./seqfile.txt --outDir ./out --outName mygraph --reference GRCh38 --panacus

--panacus takes an optional list of graph types (clip, full, filter); with no argument it defaults to the clipped graph (full if clipping is disabled), matching the other graph-type options like --gfa and --vcf. The results are written into the <outName>.stats/ subdirectory (along with the other run statistics), for example:

  • <outName>.stats/<outName>.panacus.report.html: a single self-contained interactive report with a section for every requested graph type (clip, full, filter) and count type (bp, node) — coverage histogram, growth and core-size curves, path-similarity heatmap and node distribution — that opens in any web browser
  • <outName>.stats/<outName>.panacus.histgrowth.<count>.tsv: the source coverage/growth tables, one per count type (the full and filter graphs' tables are named <outName>.full.panacus.* and <outName>.d<N>.panacus.*)

panacus renders its plots directly from the graph, so no extra plotting dependencies are needed. The count type(s) reported (bp, node, edge) can be changed with the panacus element in the configuration XML.

Please cite panacus when using these statistics.

Yeast Graph

This is a small test case whose input data is included in cactus that illustrates how to split by chromosome.

cactus-pangenome ./js ./examples/yeastPangenome.txt --reference S288C --outDir yeast-pg --outName yeast-pg --vcf --giraffe

All the output will be written to ./yeast-pg as it is computed and is described below in the context of the equivalent lower-level commands that could be used to compute the same pangenome.

Yeast: Getting Started

Below is an example of creating a yeast pangenome chromosome by chromosome, referenced on S288C.

# make the seqfile
mkdir -p yeast-pg
cp ./examples/yeastPangenome.txt yeast-pg/

# make the minigraph
cactus-minigraph ./jobstore  ./yeast-pg/yeastPangenome.txt ./yeast-pg/yeast.sv.gfa  --reference S288C

# map back to the minigraph
cactus-graphmap ./jobstore ./yeast-pg/yeastPangenome.txt ./yeast-pg/yeast.sv.gfa ./yeast-pg/yeast.paf \
--outputFasta ./yeast-pg/yeast.sv.gfa.fa  --reference S288C

Yeast: Splitting By Chromosome

Now the PAF and GFA minigraph output can be used to partition the graph and mappings based on the reference genome's (S288C's) chromosomes:

cactus-graphmap-split ./jobstore ./yeast-pg/yeastPangenome.txt ./yeast-pg/yeast.sv.gfa \
./yeast-pg/yeast.paf --outDir yeast-pg/chroms  --reference S288C

This command makes a cactus subproblem for each reference chromosome. By default, it uses all contigs in the reference. A subset can be specified using the --refContigs option.

In this example, for instance, the chrI data can be found as follows. This is everything required to run cactus-align on it as described previously.

ls -hs yeast-pg/chroms/chrI/* yeast-pg/chroms/seqfiles/chrI.seqfile 
264K yeast-pg/chroms/chrI/chrI.gfa   44K yeast-pg/chroms/chrI/chrI.paf  4.0K yeast-pg/chroms/seqfiles/chrI.seqfile

yeast-pg/chroms/chrI/fasta:
total 656K
 68K DBVPG6044.0_chrI.fa.gz   68K S288C_chrI.fa.gz   68K UWOPS034614.0_chrI.fa.gz   72K YPS128.0_chrI.fa.gz
244K _MINIGRAPH__chrI.fa      72K SK1.0_chrI.fa.gz   64K Y12.0_chrI.fa.gz

Some contigs cannot be assigned to a reference chromosome. These end up in the _AMBIGUOUS_ directory:

ls -hs yeast-pg/chroms/_AMBIGUOUS_/*
188K yeast-pg/chroms/_AMBIGUOUS_/_AMBIGUOUS_.paf

yeast-pg/chroms/_AMBIGUOUS_/fasta:
total 1.2M
   0 DBVPG6044.0__AMBIGUOUS_.fa.gz     0 S288C__AMBIGUOUS_.fa.gz  1.2M UWOPS034614.0__AMBIGUOUS_.fa.gz     0 YPS128.0__AMBIGUOUS_.fa.gz
   0 _MINIGRAPH___AMBIGUOUS_.fa        0 SK1.0__AMBIGUOUS_.fa.gz     0 Y12.0__AMBIGUOUS_.fa.gz

Here we can see that a few contigs from UWOPS034614 were left unplaced (and would be left out of any future cactus jobs).

zcat yeast-pg/chroms/_AMBIGUOUS_/fasta/UWOPS034614__AMBIGUOUS_.fa.gz | grep '>'
>id=UWOPS034614|chrXI
>id=UWOPS034614|chrX
>id=UWOPS034614|chrVII
>id=UWOPS034614|chrVIII

The reason why these contigs are unassigned to a chromosome can normally be found in minigraph.split.log:

Query contig is ambiguous: id=UWOPS034614|chrXI  len=792116 cov=0.573045 (vs 0.5) uf=1.47861 (vs 2)
 Reference contig mappings:
  chrVII: 306989
  chrXI: 453918
--
Query contig is ambiguous: id=UWOPS034614|chrVIII  len=738767 cov=0.481758 (vs 0.5) uf=1.06071 (vs 2)
 Reference contig mappings:
  chrVII: 355907
  chrVIII: 335536
--
Query contig is ambiguous: id=UWOPS034614|chrVII  len=632616 cov=0.407576 (vs 0.5) uf= infinity (vs 2)
Assigned contig to chrXI: id=DBVPG6044|chrXI  len=695907 cov=0.972054 (vs 0.5) uf= infinity (vs 2)
Query contig is ambiguous: id=UWOPS034614|chrX  len=1092164 cov=0.49082 (vs 0.25) uf=1.04957 (vs 2)
 Reference contig mappings:
  chrX: 510741
  chrXIII: 536056

This log shows that these contigs were omitted because they could not be confidently mapped (presumably due to an inter-chromosomal event) to a single reference chromosome. These filters can be turned down with the --permissiveContigFilter option (other either cactus-graphmap-split or cactus-pangenome). Using that here will assign all UWOPS034614 contigs to the reference chromosome they best align to. For example, UWOPS034614|chrXI would end up in the graph for chrXI and the approximately 40% of the contig that better aligns to chrVII would most likely end up unaligned and removed in postprocessing. The exact thresholds are explained, and can be adjusted in, the <graphmap-split> section of the cactus config.

Yeast: Batch Aligning the Chromosomes

cactus-align can be run individually on each chromosome using the seqFiles created above by way of the --batch option. When this is used, the input is expected to be the "chromfile" created by cactus-graphmap-split rather than the usual seqfile. Also, its output parameter is interpreted as an output directory, rather than an output HAL file.

cactus-align ./jobstore ./yeast-pg/chroms/chromfile.txt yeast-pg/chrom-alignments \
--batch --pangenome --reference S288C --outVG

The results are a HAL and VG file for each chromosome. The .raw.vg graphs are the direct hal2vg output: they have not been clipped, filtered or otherwise normalized, which is what distinguishes them from the <outName>.chroms/ graphs that cactus-graphmap-join makes below.

ls -hs yeast-pg/chrom-alignments/
total 72M
920K chrI.hal      1.3M chrIII.hal     1.6M chrIX.hal     1.1M chrVI.hal      1.6M chrVIII.hal     2.0M chrXI.hal      2.9M chrXIII.hal     3.6M chrXV.hal
1.3M chrI.raw.vg   1.3M chrIII.raw.vg  1.7M chrIX.raw.vg  1.2M chrVI.raw.vg   1.3M chrVIII.raw.vg  1.4M chrXI.raw.vg   2.2M chrXIII.raw.vg  3.1M chrXV.raw.vg
2.7M chrII.hal     5.1M chrIV.hal      2.0M chrV.hal      3.1M chrVII.hal     2.1M chrX.hal        3.5M chrXII.hal     2.7M chrXIV.hal      3.2M chrXVI.hal
2.1M chrII.raw.vg  4.4M chrIV.raw.vg   1.6M chrV.raw.vg   2.4M chrVII.raw.vg  1.7M chrX.raw.vg     3.1M chrXII.raw.vg  2.4M chrXIV.raw.vg   2.5M chrXVI.raw.vg

Yeast: Joining the Chromosome Alignments

As in the primates example, cactus-graphmap-join is used to make the final indexes. Its use is identical, except multiple graphs are passed as input. We also pass in the HAL files so it can merge them too.

cactus-graphmap-join ./jobstore --vg yeast-pg/chrom-alignments/*.vg --hal yeast-pg/chrom-alignments/*.hal \
--outDir ./yeast-pg --outName yeast-pg --reference S288C --vcf --giraffe clip

The GFA, VCF and all vg giraffe indexes will now be in yeast-pg:

ls -hs yeast-pg/yeast-pg*
 21M yeast-pg/yeast-pg.dist      103M yeast-pg/yeast-pg.min             4.3M yeast-pg/yeast-pg.vcf.gz
 38M yeast-pg/yeast-pg.full.hal  4.6M yeast-pg/yeast-pg.raw.vcf.gz      8.0K yeast-pg/yeast-pg.vcf.gz.tbi
 16M yeast-pg/yeast-pg.gbz       8.0K yeast-pg/yeast-pg.raw.vcf.gz.tbi
 14M yeast-pg/yeast-pg.gfa.gz     12K yeast-pg/yeast-pg.input-contig-sizes.tsv.gz
4.0K yeast-pg/yeast-pg.stats

zcat yeast-pg/yeast-pg.gfa.gz | grep '^W' | awk '{print \$1 "\t" \$2 "\t" \$3 "\t" \$4 "\t" \$5 "\t" \$6 }' | grep S288C
W	S288C	0	chrIII	0	341580
W	S288C	0	chrII	0	813597
W	S288C	0	chrI	0	219929
W	S288C	0	chrIV	0	1566853
W	S288C	0	chrIX	0	440036
W	S288C	0	chrVIII	0	581049
W	S288C	0	chrVII	0	1091538
W	S288C	0	chrVI	0	271539
W	S288C	0	chrV	0	583092
W	S288C	0	chrXIII	0	930506
W	S288C	0	chrXII	0	1075542
W	S288C	0	chrXI	0	666862
W	S288C	0	chrXIV	0	777615
W	S288C	0	chrX	0	751611
W	S288C	0	chrXVI	0	954457
W	S288C	0	chrXV	0	1091343

halStats yeast-pg/yeast-pg.full.hal --sequenceStats S288C
SequenceName, Length, NumTopSegments, NumBottomSegments
chrI, 219929, 5766, 0
chrII, 813597, 15896, 0
chrIII, 341580, 7778, 0
chrIV, 1566853, 32488, 0
chrIX, 440036, 10368, 0
chrV, 583092, 12226, 0
chrVI, 271539, 6587, 0
chrVII, 1091538, 19912, 0
chrVIII, 581049, 9996, 0
chrX, 751611, 12635, 0
chrXI, 666862, 11861, 0
chrXII, 1075542, 21991, 0
chrXIII, 930506, 16183, 0
chrXIV, 777615, 17052, 0
chrXV, 1091343, 23011, 0
chrXVI, 954457, 19805, 0

Yeast: Making a UCSC Genome Browser Assembly Hub

The HAL file can be used to produce an assembly hub has follows. Note that PYTHONPATH must set as described in Cactus's installation instructions.

hal2assemblyHub.py ./jobstore ./yeast-pg/yeast-pg.full.hal yeast-pg/hub --shortLabel yeast --longLabel "yeast pangenome"

Move yeast-pg/hub to somewhere web-accessible, and pass the full URL of yeast-pg/hub/hub.txt to the Genome Browser in the "My Data -> Track Hubs" menu. Select S288C as the reference and display the hub. Right-click on the display and select "Configure yeast track set" to toggle on all the assemblies (and toggle off Anc0 and MINIGRAPH).

MHC Graph

This is an example of how to build a MHC Pangenome using Minigraph example data of 61 haplotypes from the HPRC year 1 release. The first step is to convert the AGC archive into a seqfile:

# download the agc binary (if you don't already have it)
wget -q https://github.com/refresh-bio/agc/releases/download/v3.0/agc-3.0_x64-linux-noavx.tar.gz
tar zxf agc-3.0_x64-linux-noavx.tar.gz
export PATH=$(pwd)/agc-3.0_x64-linux-noavx/:${PATH}

# download the sequences
wget -q https://zenodo.org/record/6617246/files/MHC-61.agc

# make the seqfile
mkdir -p mhc-fa ; rm -f mhc-seqfile.txt
for s in $(agc listset MHC-61.agc); do printf "${s}\tmhc-fa/${s}.fa\n" >> mhc-seqfile.txt; agc getset MHC-61.agc $s > mhc-fa/${s}.fa; done

# (Optional) clean up the sample names. if you don't do this, adjust --reference accordingly below
sed -i mhc-seqfile.txt -e 's/^MHC-00GRCh38/MHC-GRCh38/g' -e 's/^MHC-CHM13.0/MHC-CHM13/g'

# (Optional) clean up the contig names (replacing numbers with "MHC")
for f in mhc-fa/*.fa; do sed -i ${f} -e 's/#0/#MHC/g' -e 's/#1/#MHC/g' -e 's/#2/#MCH/g'; done

The indexed pangenome can now be constructed as follows. You can set the reference to MHC-CHM13 instead of MHC-GRCh38. --mapCores 1 is used to increase parallelism. The pipeline should take about 15 minutes on a system with 8 cores.

cactus-pangenome ./js ./mhc-seqfile.txt --outDir mhc-pg --outName mhc --reference MHC-GRCh38 --gbz --giraffe --vcf --mapCores 1

GRCh38 Alts Graph

It's not often utilized as such, but the GRCh38 reference genome is actually a pangenome in and of itself, due to the various alternic loci scaffolds it contains. In general, including these sequences in the input to Cactus (both progressive and pangenome) along with the rest of GRCh38 will cause errors, so it is always important to filter them out. But, we can use them by adding them back as separate samples. Proof of concept instructions on doing so, resulting in a pangenome graph constructed from GRCh38 and its alt loci (excluding patches), are below.

The results can be found with the example pangenomes here.

We first download GRCh38 and extract each alt contig into its own fasta file:

wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
gzip -d hg38.fa.gz
samtools faidx hg38.fa
mkdir -p fa
rm -f fa/hg38_no_alts.fa
for contig in $(grep -v _alt hg38.fa.fai | awk '{print \$1}'); do samtools faidx hg38.fa $contig >> fa/hg38_no_alts.fa; done
bgzip fa/hg38_no_alts.fa --threads 16
printf "GRCh38\tfa/hg38_no_alts.fa.gz\n" > grch38_alts.seqfile
for contig in $(grep  _alt hg38.fa.fai | awk '{print \$1}'); \
do samtools faidx hg38.fa $contig | bgzip > fa/${contig}.fa.gz; \
printf "${contig}\tfa/${contig}.fa.gz\n" >> grch38_alts.seqfile ; done

Make the pangenome, setting a very permissive overlap threshold to make sure all contigs get into a chromosome as some of them are very tiny and diverse. You will need about 100Gb of RAM to run this and, if you have fewer than 32 cores, need to adjust some of the settings below.

mkdir -p work
cactus-pangenome ./js ./grch38_alts.seqfile --reference GRCh38 --gbz clip --giraffe clip --vcf --outName grch38-alts-apr13 --outDir grch38-alts --logFile grch38-alts-apr13.log --indexCores 32 --mapCores 8 --consCores 8 --refContigs $(for i in `seq 22`; do printf "chr$i "; done ; echo "chrX chrY chrM") --otherContig chrOther --permissiveContigFilter 0.05 --workDir work

The GRCh38 alt loci are concentrated in a handful of regions. Here's an example of how to pull out MHC and LRC-KIR (with a little bit of padding):

vg chunk -x grch38-alts-apr13.gbz -S grch38-alts-apr13.snarls -p GRCh38#0#chr6:28500120-33490577 -O gfa > mhc.gfa
vg chunk -x grch38-alts-apr13.gbz -S grch38-alts-apr13.snarls -p GRCh38#0#chr19:54015634-55094318 -O gfa > lrc_kir.gfa

which can then be visualized with Bandage-NG.

MHC (partial view) LRC-KIR >

HPRC Graph

The Human Pangenome Reference Consortium is producing an ever-growing number of high quality phased assemblies. This section will demonstrate how to use the Minigraph-Cactus Pangenome Pipeline to construct a Pangenome from them. Note the instructions here are slightly different than were used to create the v1.0 Minigraph-Cactus pangenome that's been released by the HPRC, as they are based on a more recent and improved version of the pipeline.

WDL / cactus-prepare support is in progress!

HPRC v1.1 Graph

Please see this page for a description of all commands to produce the v1.1 HPRC Minigraph-Cactus release. The below instructions are now out of date.

HPRC Graph: Setup and Name Munging

Important The Cactus-Minigraph Pipeline does not support alt contigs in the reference. If you really want them in your graph, then you will need to pull them out into separate samples (ie one alt contig per region per sample). Otherwise they will end up as separate reference contigs and not align together. As such we advice using the GRCh38 fasta file referenced in the hprc-${VERSION}-mc.seqfile generated below for any graph using GRCh38 as a reference.

The fasta sequences for the Year-1 HPRC assemblies are available here. We begin by using them to create an input seqfile for Cactus:

export VERSION=may4
wget -q https://raw.githubusercontent.com/human-pangenomics/HPP_Year1_Assemblies/main/assembly_index/Year1_assemblies_v2_genbank.index
grep GRCh38 Year1_assemblies_v2_genbank.index | sed -e 's/_no_alt_analysis_set\t/\t/g' | awk '{print \$1 "\t" \$2}' > hprc-${VERSION}-mc.seqfile
printf "CHM13v2\thttps://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0_maskedY_rCRS.fa.gz\n" >> hprc-${VERSION}-mc.seqfile
tail -n +2 Year1_assemblies_v2_genbank.index | awk '{print \$1 ".1\t" \$2}' | grep -v CHM13 | grep -v GRCh38 >> hprc-${VERSION}-mc.seqfile
tail -n +2 Year1_assemblies_v2_genbank.index | awk '{print \$1 ".2\t" \$3}' | grep -v CHM13 | grep -v GRCh38 >> hprc-${VERSION}-mc.seqfile
sort -k1 hprc-${VERSION}-mc.seqfile > hprc-${VERSION}-mc.seqfile.sort ; mv hprc-${VERSION}-mc.seqfile.sort hprc-${VERSION}-mc.seqfile
sed hprc-${VERSION}-mc.seqfile -i -e 's%s3://human-pangenomics/working/%https://s3-us-west-2.amazonaws.com/human-pangenomics/working/%g'

We have been holding out three samples for evaluation. This is, of course, optional, but can be done here with:

grep -v 'HG002\|HG005\|NA19240' hprc-${VERSION}-mc.seqfile > t && mv t hprc-${VERSION}-mc.seqfile

Also, a misjoin in HG02080#1#JAHEOW010000073.1 was manually corrected by using samtools faidx to break it into HG02080#1#JAHEOW010000073.1_sub_0_7238466 and HG02080#1#JAHEOW010000073.1_sub_7238466_12869124. The sub_X_Y (0-based, open-ended like BED) coordinates are understood by the pipeline, and the offsets will be preserved in the GFA W-lines at the end. If we don't apply this change, then path names with ":"'s will end up in the HAL which will prevent it from working with assembly hubs.

wget -q $(grep HG02080\.1 hprc-${VERSION}-mc.seqfile | tail -1 | awk '{print \$2}') -O HG02080.1.fa.gz
gzip -d HG02080.1.fa.gz
samtools faidx HG02080.1.fa
keep_contigs=$(awk '{print \$1}' HG02080.1.fa.fai | grep -v JAHEOW010000073\.1)
samtools faidx HG02080.1.fa ${keep_contigs} > HG02080.1.fix.fa
samtools faidx HG02080.1.fa "HG02080#1#JAHEOW010000073.1:1-7238466" | sed -e 's/\([^:]*\):\([0-9]*\)-\([0-9]*\)/echo "\1_sub_$((\2-1))_\3"/e' >> HG02080.1.fix.fa
samtools faidx HG02080.1.fa "HG02080#1#JAHEOW010000073.1:7238467-12869124" | sed -e 's/\([^:]*\):\([0-9]*\)-\([0-9]*\)/echo "\1_sub_$((\2-1))_\3"/e' >> HG02080.1.fix.fa
bgzip HG02080.1.fix.fa --threads 8

If on aws, but the revised sequence in your bucket:

aws s3 cp HG02080.1.fix.fa.gz ${MYBUCKET}/fasta/
grep -v HG02080\.1 hprc-${VERSION}-mc.seqfile > t && mv t hprc-${VERSION}-mc.seqfile
printf "HG02080.1\t${MYBUCKET}/fasta/HG02080.1.fix.fa.gz\n" >> hprc-${VERSION}-mc.seqfile

otherwise, just leave it local:

grep -v HG02080\.1 hprc-${VERSION}-mc.seqfile > t && mv t hprc-${VERSION}-mc.seqfile
printf "HG02080.1\t./HG02080.1.fix.fa.gz\n" >> hprc-${VERSION}-mc.seqfile

head -4 hprc-${VERSION}-mc.seqfile
CHM13v2 https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0_maskedY_rCRS.fa.gz
GRCh38  https://s3-us-west-2.amazonaws.com/human-pangenomics/working/HPRC_PLUS/GRCh38/assemblies/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
HG00438.1       https://s3-us-west-2.amazonaws.com/human-pangenomics/working/HPRC/HG00438/assemblies/year1_f1_assembly_v2_genbank/HG00438.paternal.f1_assembly_v2_genbank.fa.gz
HG00438.2       https://s3-us-west-2.amazonaws.com/human-pangenomics/working/HPRC/HG00438/assemblies/year1_f1_assembly_v2_genbank/HG00438.maternal.f1_assembly_v2_genbank.fa.gz

HPRC Graph: Running all at once

The new (as of v2.5.0) cactus-pangenome interface can create the pangenome from the seqfile created above in one command. To run locally (recommended 64 cores, 512Gb RAM, 3Tb disk in current directory -- the below command runs on a SLURM cluster):

cactus-pangenome ./js hprc-${VERSION}-mc.seqfile --outDir ./hprc-${VERSION} --outName hprc-${VERSION}-mc-grch38 --gbz --giraffe --vcf --chrom-vg --maxCores 64 --indexCores 32 --mapCores 8 --alignCores 16 --batchSystem slurm

To run on a Mesos cluster on AWS (initialized as described here, the command is the same as above except buckets (defined in variables above) are specified for the jobstore and outDir, and a few Toil options are provided to select instance types:

cactus-pangenome ${MYJOBSTORE} hprc-${VERSION}-mc.seqfile --outDir ${MYBUCKET} --outName hprc-${VERSION}-mc-grch38 --gbz --giraffe --vcf --chrom-vg --maxCores 64 --indexCores 64 --mapCores 8 --alignCores 16 --batchSystem mesos --provisioner aws --defaultPreemptable --nodeType r5.8xlarge:1.25,r5.16xlarge --nodeStorage 500

HPRC Graph: Running step by step

(This was the only way to run the pipeline prior to v2.5.0)

The steps below are run on AWS/S3, and assume everything is written to s3://MYBUCKET. All jobs are run on r5.8xlarge (32 cores / 256G RAM) nodes. In theory, the entire pipeline could therefore be run on a single machine (ideally with 64 cores). It would take several days though. They can be run on other batch systems, at least in theory. Most of the compute-heavy tasks spawn relatively few jobs, and may be amenable to SLURM environments.

The following environment variables must be defined: MYBUCKET and MYJOBSTORE. All output will be placed in MYBUCKET, and MYJOBSTORE will be used by TOIL for temporary storage. For example

export VERSION=may4
export MYBUCKET=s3://vg-k8s/vgamb/wg/cactus/GRCh38-f1g-90/${VERSION}
export MYJOBSTORE=aws:us-west-2:cactus-hprc-jobstore
export MINIGRAPH=https://zenodo.org/record/6499594/files/GRCh38-90c.r518.gfa.gz

HPRC Graph (step by step): Mapping to the Graph

Note: since there is already a minigraph available for this data, we just use it instead of constructing it ourselves. See the previous examples for how to construct a minigraph with cactus-minigraph.

Now that the sequences are ready, we run cactus-graphmap as before. There is a new option:

--delFilter N : Filter out mappings that would induce a deletion bubble of >N bases w.r.t. a path in the reference. If this option is used, the unfiltered paf will also be output (with a .unfiltered suffix) as well as a log detailing what was filtered and why (.filter.log suffix). This option is very important as minigraph will produce a small number of split-mappings that can cause chromosome-scale bubbles. By default, it is set to 1000000.

cactus-graphmap ${MYJOBSTORE} hprc-${VERSION}-mc.seqfile ${MINIGRAPH} ${MYBUCKET}/hprc-${VERSION}-mc-grch38.paf \
--outputGAFDir ${MYBUCKET}/gaf-hprc-${VERSION}-mc-grch38 --outputFasta ${MYBUCKET}/fasta/minigraph.grch38.sv.gfa.fa.gz \
--reference GRCh38 --mapCores 16 --delFilter 10000000  --batchSystem mesos --provisioner aws --defaultPreemptable \
--nodeType r5.8xlarge:1.5 --nodeStorage 650 --maxNodes 25 --betaInertia 0 --targetTime 1 \
--logFile hprc-${VERSION}-mc-grch38.paf.log

Note: The --betaInertia 0 --targetTime 1 options force Toil to create AWS instances as soon as they are needed.

This command uses the spot market by specifying :1.35 after the node type to bid $1.35/hr (on-demand pricing at time of writing is about $2.00).

HPRC Graph (step-by-step): Splitting by Chromosome

There are too many reference contigs to make a graph for each because of all the unplaced contigs in GRCh38. Ideally, we would drop them but it simplifies some downstream pipelines that use tools that expect them to be in BAM headers etc. to just include them in the graph. To do this, we use the --otherContig option to lump them all into a single job, and --refContigs to spell out all the contigs we want to treat separately. Note that the final output will be the same whether or not --otherContig is used. This option serves only to reduce the number of output files (and therefore alignment jobs).

cactus-graphmap-split ${MYJOBSTORE}  hprc-${VERSION}-mc.seqfile ${MINIGRAPH} \
${MYBUCKET}/hprc-${VERSION}-mc-grch38.paf --outDir ${MYBUCKET}/chroms-hprc-${VERSION}-mc-grch38 \
--otherContig chrOther --refContigs $(for i in `seq 22`; do echo chr$i; done ; echo "chrX chrY chrM") \
--reference GRCh38  --batchSystem mesos --provisioner aws --defaultPreemptable --nodeType r5.8xlarge \
--nodeStorage 1000 --maxNodes 5 --betaInertia 0 --targetTime 1 --logFile hprc-${VERSION}-mc-grch38.split.log

HPRC Graph (step-by-step): Batch Alignment

The rest of the pipeline is proceeds as in the yeast example. We need to manually download the chromfile though. We also use a new option

--maxLen N : Do not attempt to align more than N bases with the Cactus base aligner (activated with --base). This will save aligning too far into anchorless regions, which cannot be properly resolved with base alignment alone. It is 1000000 by default.

This command will create a vg and hal file for each chromosome in ${MYBUCKET}/align-batch-grch38/


aws s3 cp ${MYBUCKET}/chroms-hprc-${VERSION}-mc-grch38/chromfile.txt .
cactus-align --batch ${MYJOBSTORE} ./chromfile.txt ${MYBUCKET}/align-hprc-${VERSION}-mc-grch38 \
--alignCores 16  --pangenome --maxLen 10000 --reference GRCh38 --outVG \
--batchSystem mesos --provisioner aws --defaultPreemptable --nodeType r5.8xlarge:1.5 --nodeStorage 1000 \
--maxNodes 20 --betaInertia 0 --targetTime 1 --logFile hprc-${VERSION}-mc-grch38.align.log

HPRC Graph: Creating the Whole-Genome Graph

Important:

  • we use --filter 9 and --giraffe to make giraffe indexes on the subgraph covered by at least 9/90 haplotypes (filter does not apply to the reference).
  • we use --reference GRCh38 CHM13v2 to specify an additional reference. In this case CHM13v2 will still be clipped (but not filtered), and it will be treated as a reference path in vg (and therefore easier to query). You can add --vcfReference GRCh38 CHM13v2 to make a VCF based on CHM13 too.
  • we use --indexCores 63 to allow indexing and hal merging to be done in parallel, which can save quite a bit of time.
cactus-graphmap-join ${MYJOBSTORE} --vg $(for j in $(for i in `seq 22`; do echo chr$i; done ; \
echo "chrX chrY chrM chrOther"); do echo ${MYBUCKET}/align-hprc-${VERSION}-mc-grch38/${j}.vg; done) \
--hal $(for j in $(for i in `seq 22`; do echo chr$i; done ; echo "chrX chrY chrM chrOther"); \
do echo ${MYBUCKET}/align-hprc-${VERSION}-mc-grch38/${j}.hal; done) --outDir ${MYBUCKET}/ \
--outName hprc-${VERSION}-mc-grch38 --reference GRCh38 --filter 9 --giraffe --vcf --gbz --gfa \
--vg-chroms --batchSystem mesos --provisioner aws --defaultPreemptable --nodeType r5.16xlarge \
--nodeStorage 1000 --maxNodes 1 --indexCores 63  --logFile hprc-${VERSION}-mc-grch38.join.log

All sequences clipped out by cactus-graphmap-join will be saved in a BED file in the ".stats.gz" file in its output directory.

HPRC Graph: Changing the Reference

The selection of the reference genome is very important, as it will be used as the backbone for the graph. It is the only genome that is guaranteed to not have any cycles nor to ever be clipped, and therefore provides a coordinate system in the graph. Any input genome can be used as a reference, provided it's consistently passed as the --reference option to all the commands. It also must not have a "." in its genome name. In practice, there are usually two possible references for the HPRC graphs: GRCh38 and CHM13.

It is advisable to also pass --reference CHM13v2 GRCh38 --vcfReference CHM13v2 GRCh38 to cactus-graphmap-join to tell it to make a second VCF based on the GRCh38 reference.

Note: some contig names like chrY (if it is not included) and options like --otherContig will not be necessary for CHM13

HPRC Graph: Other Approaches for Masking or Clipping out Complex Regions

The Pangenome Pipeline supports options to for special handling of masked regions at pretty much every step. These were added to address various issues during initial development and testing. The approach of just aligning everything and filtering based on the minigraph described above is much simpler and seems at least as effective.

cactus-preprocess

Most satellite sequence can be detected with dna-brnn, which can be run with via the --maskAlpha --minLength 100000 --brnnCores 8 options in cactus-preprocess. The entire pipeline supports sub-sequence fragments via naming conventions, so the masked sequence can be clipped out instead of masked by using --clipAlpha instead of --maskAlpha

The minigraph mappings themselves can also be used to derive regions to mask, by finding gaps in the alignments. This can be done by passing a PAF file (output from cactus-graphmap) back into cactus-preprocess via the --maskFile option. This option can also accept BED files to mask any user-specified regions. When using this option, the --maskAction option can be used to specify whether masked sequence is clipped out or not.

cactus-graphmap

Softmasked input sequence can be ignored by using the --maskFilter 100000 option. This will force such sequence to remain unaligned.

cactus-graphmap-split

Softmasked input sequence can (and should) be ignored when computing coverage in order to assign contigs to reference chromosomes. This is done with --maskFilter 100000

cactus-align

Softmasked input can be ignored (and forced to stay unaligned) with the --barMaskFilter 100000 option to cactus-align.

HPRC Version 1.0 Graphs

These graphs were created with the cactus-pangenome.sh script using Cactus commit 6cd9a42cdf40ad61843664ed82c9d5bc26445570. The seqfile input was constructed as above, except chrY was only added to CHM13 for the CHM13-based graph (and chrEBV was never added). Instead, a decoy graph consisting of chrEBV and all the hs38d1 contigs was added to both graphs in the cactus-graphmap-join step.

The other main differences between this pipeline and ${VERSION} are

  • Input fasta files were softmasked with dna-brnn regions >100kb
  • After mapping to the graph, minimizer gaps >100kb were masked using a second call to cactus-preprocess
  • The two sets of masked regions were merged together and clipped out the input sequences.
  • The clipped sequences were remapped to the graph once again and the pipeline continued from there
  • The --base option was never used to perform sequence-to-graph base alignment (it didn't exist)
  • The --delFilter option didn't exist either, so several large spurious bubbles made it into the graphs
  • Much more stringent options were used to assign contigs to chromosomes with cactus-graphmap-join. This was possible to some extent because the contigs were clipped, but also caused more sequence to be classified as ambiguous.
  • cactus-graphmap-join clipped out sequence that was unaligned to anything else (including minigraph), rather than unaligned to minigraph. (this is less stringent).
  • A few small bugs were fixed in Cactus between the two versions, notably one that caused erroneous tiny duplications and inversions.

GRCh38 graph command line

./cactus-pangenome.sh -j aws:us-west-2:glennhickey-jobstore7 -s ./hprc-year1-f1g.fix.HG02080.1.brnn.leaveout.seqfile \
-m ftp://ftp.dfci.harvard.edu/pub/hli/minigraph/HPRC-f1g/GRCh38-f1g-90.gfa.gz  \
-o s3://vg-k8s/vgamb/wg/cactus/GRCh38-f1g-90/aug11 -n GRCh38-f1g-90-mc-aug11 -r GRCh38 \
-d s3://vg-k8s/vgamb/wg/fasta/hs38d1.decoys.only.vg  -g  -F  -C -M 100000 -K 10000  2>> stderr.aug11.2.log > /dev/null

CHM13 graph command line

 ./cactus-pangenome.sh -j aws:us-west-2:glennhickey-jobstore-hprc4 -s \
 ./hprc-year1-f1g.chmy.fix.HG02080.1.brnn.leaveout.seqfile \
 -m ftp://ftp.dfci.harvard.edu/pub/hli/minigraph/HPRC-f1g/CHM13-f1g-90.gfa.gz  \
 -o s3://vg-k8s/vgamb/wg/cactus/CHM13-f1g-90/aug11 -n CHM13-f1g-90-mc-aug11  \
 -r CHM13 -v GRCh38 -d s3://vg-k8s/vgamb/wg/fasta/hs38d1.decoys.only.vg  -g  \
 -F -C -M 100000 -K 10000 -y  2>> stderr.aug11.chm13.3.log > /dev/null

Tutorials

Frequently Asked Questions

Q: Do my input FASTA files need to be masked?

A: No. While Progressive Cactus requires input repeats to be softmasked, Minigraph-Cactus does not use masking information. Whether or not your input is softmasked will have zero effect on the output. It is strongly advised not to hardmask your input!

Q": cactus-graphmap-join keeps crashing with a segfault and/or running out of memory.

A: This is usually because vg index -j crashes while computing the distance index for --giraffe. Make sure you did not disable clipping or filtering. Try to use --filter N where N is about 10% of your input samples. If that fails, try on a system with more memory.

Q: Why are the node id's different in my GFA and GBZ

A: See Node Chopping section above. Note: this shouldn't happen when running v2.9.1 or later.

Q: My tools can't read GFA 1.1 and W-lines.

A: As mentioned above, you can use vg convert -gfW to convert from GFA 1.1 to GFA 1.0.

Q: I hate the idea of clipping and filtering sequence from my graph. Why do I have to do it?!

A: So current toolchains can work with your graphs. But clipping and filtering is optional in cactus-pangenome (and cactus-graphmap-join), just specify full for the various outputs (keeping in mind that much satellite sequence won't be aligned to anything).

Q: I get an error to the effect of ERROR: No matching distribution found for toil[aws]==xxxx when trying to install Toil.

A: This is probably happening because you are using Python 3.6. Toil and Cactus require Python >= 3.7. Use python3 --version to check your Python version.

Q: cactus-align-batch spawns too many cactus-align jobs and runs out of memory. How do I fix this?

A: You can control the number of jobs with --alignCores and --maxCores which set the cores per align job and total cores, respectively. So to only do two align jobs at a time using 8 cores total, you can set --alignCores 4 --maxCores 8. Update: cactus-align-batch is now deprecated -- use cactus-align --batch instead and control the number of jobs with --consCores, just like in progressive cactus.

Q: The node IDs referred to in the output VCF don't match the GBZ!

A: Indeed they do not (prior to v2.9.1). They refer to the (unchopped) GFA IDs. Please see the Node Chopping section above.

Q: Some contigs or even entire samples are getting mysteriously dropped from my output. What happened?

A: This is probably due to the reference contig assignment thresholds. The defaults (found in minQueryCoverages and minQueryCoverageThresholds in the cactus_progressive_config.xml configuration file) are quite stringent for tiny contigs. For example, with the current defaults (circa v2.6.7), a contig shorter than 100kb would need to map with at least 75% of its bases to a reference graph component in order to be included. For diverse inputs and / or very fragmented assemblies, this may be too strict. You can inspect which of your contigs were filtered and why by looking at chrom-subproblems/minigraph-split-log in your output directory. In most cases, you can resolve this by using the --permissiveContigFilter option, which by default, applies a 25% threshold to all contig sizes (you can further lower it by passing in a value, ex --permissiveContigFilter 0.1. This option (available in cactus-pangenome and cactus-graphmap-split) is also mentioned in the Yeast example above. **

Q: Can I adjust the minigraph parameters?

A: Yes, they can be tuned via the <graphmap> element in the configuration XML. The construction parameters are found in the minigraphConstructOptions attribuate, and the mapping parameters in minigraphMapOptions.