--ref-allele force changes phasing information in VCF file

[JMR JMR]

I am using the `--ref-allele force` command to swap alleles in a VCF file. The VCF file has both PHASED genotypes (it's actually 1KG) and UNPHASED genotypes. I noticed that after using this command some unphased genotypes appear as phased.

For example, keeping just one sample from my VCF file and estimating the frequency of the 0/0 0/1 and 0/1 genotypes before the --ref-alleles gives me:

   0/0    0/1    1/1
167361   1494  19030

And after:

   0/0    0/1    0|0    1/1    1|1
167317   1494   1744  17120    210

I noticed that only homozygous genotypes have phased information now, which is fine. I just want to know if this is the default/intended behavior of this command. What will happen if all my samples are phased? I assume the original phase will be retained? I saw on this post that there were some issues with this command before but appear to be solved using the latest version, which I am using.

Thank you.

JMR

Log file info:

PLINK v2.00a2LM 64-bit Intel (5 Jul 2019)      www.cog-genomics.org/plink/2.0/
(C) 2005-2019 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to test.log.
Options in effect:
  --export vcf-iid
  --out test
  --ref-allele force test_chr22_alleles_2swap.txt 2 1
  --vcf test_chr22_4AAswap.vcf.gz

Start time: Mon Jul  8 14:28:09 2019
Note: --export 'vcf-iid' modifier is deprecated.  Use 'vcf' + 'id-paste=iid'.
64453 MiB RAM detected; reserving 32226 MiB for main workspace.
Using up to 12 threads (change this with --threads).
--vcf: 187885 variants scanned.
--vcf: test-temporary.pgen + test-temporary.pvar + test-temporary.psam written.
312 samples (0 females, 0 males, 312 ambiguous; 312 founders) loaded from
test-temporary.psam.
187885 variants loaded from test-temporary.pvar.
Note: No phenotype data present.
--ref-allele: 21409 sets of allele codes rotated.
--export vcf to test.vcf ... done.
End time: Mon Jul  8 14:28:11 2019

[JMR JMR]

This is the post on biostars where I noticed the --ref-allele force command was having some issues before treating heterozygous:

https://www.biostars.org/p/271765/#388641

[Christopher Chang]

No, because it's irrelevant.  Homozygous genotypes don't have a phasing state at all, so plink2 renders homozygous genotypes in VCFs with '|' iff the previous heterozygous genotype was phased.

(There are plans to add support for VCF "phase sets" in the future, but for now, you should stick to bcftools when it's important to preserve that information.)

[Christopher Chang]

Clarification: "No" was in response to the biostars question of whether this behavior can be changed; I cut-and-pasted too quickly.  Phase information should be retained properly; if you have an example where that isn't true, please post it and I'll fix the issue.

[JMR JMR]

Thank you for the reply.

Yes, homozygous genotypes don't have a phasing state so having a "|" is not gonna affect my analysis.

I have another question now. If you notice from the table of frequency of genotypes before the --ref-allele force command, all genotypes where unphased. So, there where no previous phased heterozygous genotypes on that sample. Why did plink2 still gave me phased genotypes?

[Christopher Chang]

Because, if the dataset is partially or fully phased, the homozygous calls *before* the first het on each chromosome are rendered as phased.

[JMR JMR]

Got it! Thanks!