How to build an index for an mpmap without using VCF and FASTA files

[xiaoxinhan]

Hello,
I want to use vg mpmap for RNA data alignment, but the reference genome I'm using is nonlinear. I do this with codes below:

vg autoindex --thread 40  --workflow mpmap --prefix graph_gfa/ath --gfa data/NNJ.combine.gfa -T graph_gfa/ --hap-tx-gff data/Col_exon_graph.gff3 

And I got an error report like this:

vg [warning]: System's vm.overcommit_memory setting is 2 (never overcommit). vg does not work well under these conditions; you may appear to run out of memory with p
[vg autoindex] Executing command: vg autoindex --thread 40 --workflow mpmap --prefix graph_gfa/ath --gfa data/NNJ.combine.gfa -T graph_gfa/ --hap-tx-gff data/Col_exo
[IndexRegistry]: Checking for haplotype lines in GFA.
error:[vg autoindex] Input is not sufficient to create indexes
Inputs
        Haplotype GTF/GFF
        Reference GFA w/ Haplotypes
are insufficient to create target index Spliced Distance Index

Besides,I ues vg index to built an index,likes:

vg index -p -t 20 -x graph/ath.xg graph/ath_spliced.pg
vg index -p -t 20 -g graph/ath.gcsa graph/ath_spliced.pg
vg index -p -t 20 -j graph/ath.dist graph/ath_spliced.pg

But there is only one .xg file as a result, without the gcsa and dist indexes.

[adamnovak]

It looks like to make a "Spliced Distance Index" you need a "Spliced XG", which you get from a "Spliced VG w/ Transcript Paths", which in turn comes from one of these recipes:

vg/src/index_registry.cpp

Lines 3256 to 3272 in 3a1c22a

	auto vg_rna_gbz_graph_only =
	registry.register_recipe({"Spliced MaxNodeID", "Spliced VG w/ Transcript Paths"}, {"GBZ", "GTF/GFF", "Haplotype GTF/GFF"},
	[gbz_vg_rna](const vector<const IndexFile*>& inputs,
	const IndexingPlan* plan,
	AliasGraph& alias_graph,
	const IndexGroup& constructing) {
	return gbz_vg_rna(inputs, plan, alias_graph, constructing);
	});
	auto vg_rna_gbz_liftover_graph_only =
	registry.register_recipe({"Spliced MaxNodeID", "Spliced VG w/ Transcript Paths"}, {"GBZ", "GTF/GFF"},
	[gbz_vg_rna](const vector<const IndexFile*>& inputs,
	const IndexingPlan* plan,
	AliasGraph& alias_graph,
	const IndexGroup& constructing) {
	return gbz_vg_rna(inputs, plan, alias_graph, constructing);
	});

The "GBZ" can be made from the "Reference GFA w/ Haplotypes", which you have. But both these recipes need a plain "GTF/GFF", not (just) a "Haplotype GTF/GFF".

So it looks like you need to use --tx-gff and use one reference's gene models instead of (or in addition to) gene models for each haplotype with --hap-tx-gff. (@faithokamoto, would it make sense for us to have a recipe that only used the per-haplotype transcripts to splice the graph?)

[faithokamoto]

We could do that, yeah.

[xiaoxinhan]

@adamnovak
I used the vg file generated by PGGB together with vg rna to construct a spliced vg file, and then used vg index to build the xg, gcsa, and dist indexes.
But,the sizes of the generated indexes differed substantially, and the construction of the gcsa index resulted in memory exhaustion.likes,

vg mod -X 256 ath_spliced_new.vg >ath_new.vg

vg index -p -t 40 -x ath_new.xg ath_new.vg
vg index -p -t 10 -g ath_new.gcsa ath_new.vg
vg index -p -t 40 -j ath_new.dist ath_new.vg

its results is

Is this result normal?

[adamnovak]

It's fairly common for complex graphs to have large distance indexes, and to not be able to make GCSA indexes for them without aggressive pruning.

vg autoindex automates some of the backing-off and applying progressively more aggressive pruning. If you want to make the GCSA with vg index, you will probably need to follow https://github.com/vgteam/vg/wiki/Index-Construction#complex-graph to prune the graph you are indexing first yourself. If that isn't enough to make the GCSA indexing problem fit in your memory and disk space, you can start increasing the -k and maybe decreasing the -e option values for vg prune to make it remove more stuff.