VCF output has chromosome symbol "A" but contig ID "1" in header

[pmckenz1]

Hi, I'm writing to report what appears (?) to be a bug in VCF output from SLiM 5.2.

In the VCF produced by the script below, the contig declaration uses the chromosome ID "1":

##contig=<ID=1,URL=https://github.com/MesserLab/SLiM>

However, the CHROM field in the variant records uses the implicit chromosome's symbol, "A":

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ...
  A      3   .  A   T   1000 PASS   ...

I expected the declared contig ID to match the value used in the CHROM field.

The same header/record mismatch also appears in the example VCF output on page 199 of the updated SLiM 5.2 manual.

You can reproduce this with the script below.

initialize() {
    setSeed(12345);
    initializeMutationRate(1e-3);
    initializeMutationType("m1", 0.5, "f", 0.0);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 999);
    initializeRecombinationRate(1e-4);
}

1 early() {
    sim.addSubpop("p1", 20);
}

20 late() {
    p1.outputVCFSample(
        4,
        replace=F,
        outputMultiallelics=F,
        filePath="slim_vcf_header_bug.vcf"
    );
    community.simulationFinished();
}

Thanks!
Patrick McKenzie

[bhaller]

Thanks for the report! Yes, this sounds like it is probably a bug.

[pmckenz1]

Thanks! Very easy to work around for now.

[bhaller]

marking "priority" because I intend to fix this in the SLiM 6 beta.

[bhaller]

Hi @pmckenz1! OK, looking into this, it looks like the ##contig ID is simply 1, always; it is hard-coded into SLiM's VCF output. Obviously this was overlooking when multiple-chromosome support was added to SLiM. :-> Specifically, the line emitted is always:

##contig=<ID=1,URL=https://github.com/MesserLab/SLiM>

So, I'd like a bit of clarification as to what the best behavior would be for multi-chromosome SLiM now. I take it you want the ID to be the symbol property of the chromosome. And I guess in a multi-chromosome model that uses outputIndividualsToVCF() to output data from multiple chromosomes to a single VCF, you'd like to see a separate ##contig line emitted for each chromosome emitted, is that correct? And what about the URL field? I am tempted to add a contigURL property to the Chromosome class that allows the user to configure this; by default it would be set to https://github.com/MesserLab/SLiM but the user could change it, and then see the resulting URL in their VCF output. I guess that would be nice to allow.

Does all of this sound correct? Anything else that should be on my radar as far as multichromosome VCF output is concerned?

[bhaller]

Fixed; closing.