-
Notifications
You must be signed in to change notification settings - Fork 219
Graph Coordinates (gRef)
Pangenome Graph Reference (GRef) coordinates give canonical positions to graph locations that do not lie on a linear reference. Their primary application is representing variation within large inserted sequences. They also serve as a bridge for using graph coordinates within standard linear-reference formats such as VCF and BAM. In vg, GRef coordinates are stored as a covering set of reference-sense paths. They include a "base" reference path alongside smaller paths that denote various inserted sequences. GRef coordinates are intended to be usable as a drop-in replacement for existing reference samples, and are identified with a gref_ prefix. For example, you could toggle between the base and GRef reference by using vg surject -n CHM13 or vg surject -n gref_CHM13.
GRef coordinates can be added to a graph in .vg format using vg paths -u. For example, chr1.vg may have reference paths
CHM13#0#chr1
HG002#1#chr1
HG002#2#chr1
HG003#1#chr1
HG003#2#chr1
Running vg paths -x chr1.vg -u -Q CHM13 --gref-segs chr1.gref.segs.tsv > chr1.gref.vg could produce a graph with
CHM13#0#chr1
HG002#1#chr1
HG002#2#chr1
HG003#1#chr1
HG003#2#chr1
gref_CHM13#0#chr1
gref_CHM13#0#chr1_1_alt
gref_CHM13#0#chr1_2_alt
gref_CHM13#0#chr1_3_alt
(etc)
The GRef paths are computed greedily, with the objective of maximizing their lengths. By default, only paths with length >= 50 are stored. The threshold can be changed with vg paths -l. All input paths are preserved, and gref_CHM13#0#chr1 is the same as CHM13#0#chr1 in everything but name. In other words, gref_CHM13 is a superset of CHM13. Each off-reference path (with an _alt suffix) corresponds to a contiguous interval of one of the other input paths. In this sense, they are similar to rGFA, and to the alternate contigs of GRCh38 (hence the _alt suffix). The mapping between the GRef path and its source interval is stored in chr1.gref.segs.tsv.
GRef paths are acyclic, and their computation is currently limited to graphs with acyclic reference paths such as those from minigraph-cactus. There is a higher-level interface in Cactus for GRef computation.
Everything below assumes you are working with a graph that has GRef coordinates embedded. So vg paths -u is a prerequisite to all GRef functionality.
vg paths --gref-segs produces a BED-like table that provides a mapping between the graph and the GRef coordinates. It is tab-separated; the columns are padded here to line up:
#source_path source_start source_end gref_contig level strand ref_contig ref_start ref_end top_level_snarl
HG03130#2#CM087266.1 2081065 2416307 gref_CHM13#0#chr1_1_alt 1 + CHM13#0#chr1 2095656 2253401 >909923>1299637
NA20282#1#JBIRDT010000003.1 125491785 125756875 gref_CHM13#0#chr1_2_alt 1 + CHM13#0#chr1 127205277 127333597 >14349205>15944529
HG03688#1#CM086800.1 230925228 231178279 gref_CHM13#0#chr1_3_alt 1 + CHM13#0#chr1 227747005 228024746 >8915034>9040880
The first four columns map the original path interval to the GRef path. The next column stores the nesting level, where 1 means the interval branches directly off the base reference. The strand column is the strand of the GRef path relative to the original interval. It will be positive in most cases but the format supports both directions (though only positive is rGFA compatible). The last 4 columns map the interval to its position at the 0th level (ie the linear reference path), alongside the name of the top level snarl.
GRef coordinates can be used to create VCF files with nested variation inside insertions. No changes to the VCF format are required, though some INFO tags are included to help place off-reference contigs onto the base reference (more below). Nested output is activated with -a; without it you get a flat VCF over the GRef contigs, with nothing nested inside them:
vg deconstruct chr1.gref.vg -P gref_CHM13 -C -a > chr1.gref.vcf
-P gref_CHM13 selects the GRef paths as the references, -a reports nested snarls as well as top-level ones, and -C (optional) shortens CHROM from the full path name to just the contig.
Records whose CHROM is a base contig describe variation against the linear reference, exactly as they would without GRef. Records whose CHROM is an _alt contig describe variation inside an inserted sequence, in that insertion's own coordinates. These are variants that are effectively invisible without GRef.
#CHROM POS ID REF ALT INFO
chr1 1 >1>6 C CAAG,CATG ...;LV=0;CH=0;RC=chr1;RS=1;RD=2
chr1_1_alt 2 >2>5 A T ...;LV=0;CH=1;PS=>1>6;RC=chr1;RS=1;RD=2
Nesting is described by six INFO tags:
-
LV— level in the snarl tree, counting only ancestors whose record is on this record's own CHROM. 0 means top level for this contig, so a record at the start of an insertion hasLV=0even though it is nested. -
CH— how many times CHROM changes on the way up, i.e. how many insertions deep the record is. A record directly inside an insertion off the base reference hasCH=1. This is the tag to filter on. -
PS— the ID of the parent snarl's record. -
RC,RS,RD— the contig, start and end of the topmost enclosing site. For a nested record these give a linear-reference anchor for the whole nest.
So bcftools view -i 'INFO/CH==1' selects everything exactly one insertion deep, and INFO/CH==0 is the plain linear-reference VCF.
RC/RS/RD are what make a region query on the linear reference reach inside insertions. A plain bcftools view -r chr1:2000000-2500000 returns only records literally on chr1, since a nested record lives on an _alt contig with its own coordinates. Filtering on the anchor instead picks up both:
bcftools view -i 'INFO/RC="chr1" && INFO/RS>=2000000 && INFO/RD<=2500000' chr1.gref.vcf.gz
This returns the on-reference records in the interval and every record nested inside an insertion anchored there. Add && INFO/CH>0 for just the nested ones. (Records get RC even when they are top level, where it is their own CHROM; the tags are only absent for a fragment with no enclosing site at all.)
Nested VCFs carry redundant information in that the off-reference variants are buried inside insertion alleles they are contained in. This can be simplified using the -L flag, which merges similar alleles (of at least 50bp, change with --cluster-min-len). For example, L 0.95 will merge alleles if they are more than 95% similar. This results in a much smaller, cleaner VCF where resolution is lost at large SV sites, but represented on the nested contigs.
vg call can output nested GRef VCFs as well, using the -A option. They have the same tags as described above for vg deconstruct. All sites are genotyped independently with -A, but --top-down can be used to genotype top-level sites then propagate consistent genotypes to their children. vg call supports the same -L allele clustering described above.
vg pack -x chr1.gref.vg -g aln.gam -o aln.pack
vg call chr1.gref.vg -k aln.pack -S gref_CHM13 -A > chr1.gref.calls.vcf
#CHROM POS ID REF ALT INFO FORMAT SAMPLE
chr1 1 >1>6 C CATG,CAAG AT=...;DP=149;LV=0;CH=0;RC=chr1;RS=1;RD=2 GT 1/2
chr1_1_alt 2 >2>5 A T AT=...;DP=300;LV=0;CH=1;PS=>1>6;RC=chr1;RS=1;RD=2 GT 1/0
Note that vg call shortens CHROM to the contig name on its own, so its output matches vg deconstruct -C rather than plain vg deconstruct. Surjected BAMs keep the full path name, so all three do not agree by default -- worth checking before joining them.
GRef coordinates can be used to make BAM / FASTA pairs. In this context, they resemble alt contigs from GRCh38, with three differences: they never overlap each other, they cover the whole genome, and they are explicitly linked to the base reference through the graph. In theory, you can pass a GRef BAM to any tool that operates on BAM but in practice the high fragmentation will often cause issues. For this reason, you may want to filter down to a minimum contig length, either upstream with vg paths -l or on the BAM / FASTA itself with custom scripts. I used a 1kb minimum for FreeBayes with the HPRC graphs for example (and this still wasn't sufficient for DeepVariant to work properly). Besides that, there is no interface change required, just use the gref_ version of the reference instead of the base version. Ex:
vg surject -x chr1.gref.vg -n gref_CHM13 -N HG002 -b aln.gam > aln.bam
A read that aligns inside a large insertion surjects onto the _alt contig at its own offset, instead of being left unmapped or pushed onto the nearest linear position. The _alt contigs appear in the BAM header next to the base ones, so downstream tools see them as ordinary extra contigs:
@SQ SN:gref_CHM13#0#chr1 LN:248387328
@SQ SN:gref_CHM13#0#chr1_1_alt LN:335242
For tools that need a FASTA of the reference, extract the same path set:
vg paths -x chr1.gref.vg -S gref_CHM13 -F | bgzip > chr1.gref.fa.gz
samtools faidx chr1.gref.fa.gz