Invariant sites and indels

What logically follows from this description is that genome polarisation is meaningful only for variant sites (markers). Invariant sites will have low support and, though they will not change the diem outcome, they will slow down convergence and obscure the nature of change between taxa. During diem iterations their obscuring effect is removed by rescaling the hybrid indices onto \([0,1]\). Including invariant sites in genome polarisation is acceptable, and likely unavoidable. A proportion of truly invariant sites will be mis-diagnosed as variant due to sequencing errors.

At this time, genome polarisation uses variant markers where alleles differ in nucleotide substitutions. Insertions or deletions (indels) are not allowed to be classified among the most frequent alleles in the vcf2diem function.

Sites not informative for genome polarisation

Along side invariant sites, some variant sites are also uninformative for genome polarisation. Sites that include only missing data and heterozygous genotypes are unpolarisable. In general applications, one might wish to exclude sites that are only included in a vcf due to the presence of a single heterozygous individual.

Uninformative sites will have support equal to 0, and their polarity will thus be equal to the null polarity. Including uninformative sites has similar consequences to including invariant sites. When a user wishes to retain such sites, we advice to then filter polarised markers with high support (or simply high diagnostic index, which is strongly correlated with support) for subsequent analyses and interpretation.

Haploid data

Haploid data, including markers on sex chromosomes of the heterogametic sex, are usually shown as sequence alignments. To illustrate the principle of genome representation used in diem, let’s have an 8bp-long alignment with five individuals.

5 8
ind1   AACCTTGG
ind2   TACGATGG
ind3   TACC-TGT
ind4   CACGTTGG
ind5   AACCTTGT

We can see that the example alignment contains four variant markers at sites 1, 4, 5, and 8. Marker 5 contains a deletion in ind3, while other markers vary in different nucleotide substitutions. In marker 1, two individuals have an allele A, two individuals have an allele T and one has a C. This marker can be resolved in terms of two most frequent alleles A and T, where the individual ind4 has a third allele and is therefore considered as having an unknown genomic state. To convert the marker into diem format, we must make a decision on which allele will be encoded as 0. The choice is arbitrary. Good practice is to keep a record of genotype encoding. Let’s decide that all alleles in ind1 will be encoded as 0. Marker 1 will then be 022_0.

Contrary to DNA sequence alignments, diem format transposes the data. Rows in diem input files are markers and columns represent individual samples. The diem output for the example alignment will thus have four lines representing the four variant markers, and five columns representing the individuals.

Marker1   022_0
Marker4   02020
Marker5   02_00
Marker8   00202

The diem format does not store information on marker identification or location. This metadata must be saved separately in a ‘bed-like’ file (with at minimum the scaffold and position for each site). diem uses a prefix S at line starts to ensure that all operating systems always read genotypes as strings without attempting to convert them into numbers. Encoding of markers 4 and 8 shown above is at risk of being interpreted as numeric. The final file in diem format will be:

S022_0
S02020
S02_00
S00202

Diploid data

# Import the first converted file for all individuals 
# and without changing marker polarity
myotis <- importPolarized("myotis-01.txt", 
                          changePolarity = rep(FALSE, 11), 
                          ChosenInds = 1:14)
# Check whether a marker includes heterozygotes 
# or that there is more than one genotype
apply(myotis, MARGIN = 2, 
      FUN = \(x) any(grepl("1", x)) | sum(levels(factor(x)) %in% c("0", "1", "2")) > 1) 
#    m1    m2    m3    m4    m5    m6    m7    m8    m9   m10   m11 
# FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
         

Polyploid data

The same principles apply to higher ploidy datasets. Preparing data in diem format for genome polarisation includes these steps:

1. Remove indels.
2. Identify two most frequent alleles for each site.
3. Encode homozygous and heterozygous (mixed) genotypes that include only the two most frequent alleles.

Note that vcf2diem only converts diploid datasets at this time, and should not be used to convert polyploid vcf files.