Roeland E. Voorrips, 2024-10-16

This vignette shows how to use the main functions provided in R package PolyHaplotyper and explains their output.

The PolyHaplotyper package contains tools to derive phased haplotypes from unphased bi-allelic marker (SNP) data in collections of individuals of any even ploidy level: diploid, tetraploid, hexaploid. The haplotyping does not rely on linkage mapping and often is very fast (within a few seconds per haploblock).

The markers must have been grouped into “haploblocks” that are so tightly linked that recombination may be assumed not to occur within the population. Typically such a haploblock would be the collection of all SNPs within a single contig. Since the number of possible haplotypes and hence the computation time increases dramatically with increasing number of markers it is best to split larger haploblocks to a maximum of 7-8 markers (in tetraploids) or 5-6 markers (in hexaploids).

The haplotyping is more reliable if one or more FS (Full-Sib) families, with parents, are present, but will also proceed if no such FS families are available, or if the parents have missing data.

Haploblock definitions

The markers must be assigned to haploblocks which are stored as a list. In the example the haploblock information is contained in an extra column “contig” in the marker dosage data:

demo_snpdos[1:6, 1:8]

##           marker contig s488 s649 s487 s639 s1_3 s1_4
## 1 WagCmDlf_18464 ctg001    1    1    0    0    1    1
## 2 WagCmDlf_12697 ctg001    1    1    1    1    0    0
## 3 WagCmDlf_61726 ctg001    6    6    6    6    6    6
## 4 WagCmDlf_26276 ctg001    1    1    1    1    0    0
## 5 WagCmDlf_33716 ctg002    6    6    6    6    6    6
## 6 WagCmDlf_19610 ctg002    5    5    6    6    6    6

We convert the haploblock information to list format:

hblist <- haploblock_df2list(demo_snpdos, mrkcol=1, hbcol=2)
# number of markers in each haploblock:
sapply(hblist, length)

## ctg001 ctg002 ctg003 ctg004 ctg005 ctg006 ctg007 
##      4      4      5      5      4      4      4

This small example has data for 7 haploblocks, each with 4 or 5 markers.

Marker data

The marker data must be supplied as a matrix or data.frame, with markers in rows and individuals in columns, with each cell containing the dosage (0 … ploidy) of the target marker allele. Our data are in data.frame demo_snpdos, which, as we saw above, contains an extra column “contig” defining the haploblocks. We remove the contig column from demo_snpdos to get it in the required format:

demo_snpdos <- demo_snpdos[, -2]
demo_snpdos[1:6, 1:8]

##           marker s488 s649 s487 s639 s1_3 s1_4 s1_5
## 1 WagCmDlf_18464    1    1    0    0    1    1    1
## 2 WagCmDlf_12697    1    1    1    1    0    0    0
## 3 WagCmDlf_61726    6    6    6    6    6    6    6
## 4 WagCmDlf_26276    1    1    1    1    0    0    0
## 5 WagCmDlf_33716    6    6    6    6    6    6    6
## 6 WagCmDlf_19610    5    5    6    6    6    6    6

The SNP dosages are all in the range 0…6, as these demo data are obtained from a hexaploid crop (chrysanthemum). NA indicates an unknown dosage.

Pedigree

The pedigree must be a data.frame or matrix with in the first 3 columns the names of the individual, its mother and its father. Parents may be missing (NA) but individuals may not, and there may not be more than one line for the same individual. Additional columns may be present. In our case demo_ped has a fourth column containing the sample nr:

head(demo_ped)

##   genotype mother father sample_nr
## 1     2096     NA     NA      s629
## 2     4431     NA     NA      s630
## 3     6352   8420     NA      s631
## 4     7688     NA     NA      s632
## 5     7783   4431     NA      s201
## 6     8309     NA     NA      s611

In this example several individuals are represented by more than one line, because here the pedigree is also used to specify the (sometimes duplicated) samples for each individual.
In order to remove the duplicates from the pedigree we build a list with the samples representing each replicated individual.

# create a list of of replicates:
tb <- table(demo_ped$genotype)
replist <- list()
for (dup in names(tb[tb>1])) {
  replist[[dup]] <- demo_ped$sample_nr[demo_ped$genotype == dup]
}

We remove the replicated individuals from the pedigree, keeping only the first row for each set of duplicates:

dim(demo_ped)

## [1] 661   4

for (dup in seq_along(replist)) {
  dupsamp <- as.character(replist[[dup]])[-1]
  demo_ped <- demo_ped[!(demo_ped$sample_nr %in% dupsamp),]
}
dim(demo_ped)

## [1] 631   4

We see that 30 lines containing duplicates have been discarded.

In demo_snpdos, the data.frame containing the SNP dosages, the column names are the sample names, and duplicated samples are still present. First we merge the dosage data of the replicates for each individual, again using replist (the list of replicates):

# merge all the duplicated samples:
dim(demo_snpdos)

## [1]  30 662

demo_snpdos <- mergeReplicates(mrkDosage=demo_snpdos, replist=replist,
                               solveConflicts=TRUE)
dim(demo_snpdos)

## [1]  30 631

Function mergeReplicates has converted the demo_snpdos data.frame to a matrix, with the rownames taken from the column “marker”. For each individual present in multiple replicates only the first column name is retained, and the merged data are the consensus scores over all replicates. We see that demo_snpdos now has 31 less columns: 1 was the removed “marker” column and 30 were columns for the duplicated samples, corresponding to the 30 duplicates removed from the pedigree.

Now we must change the sample numbers to the corresponding individual names:

colnames(demo_snpdos) <- 
  demo_ped$genotype[match(colnames(demo_snpdos), demo_ped$sample_nr)]

FS families

Finally we specify the four FS families and their parents (note that individual 39287 is the father of two FS families):

parents <- cbind(c(36451, 41234, 9656, 32141),
                 c(39287, 40360, 9541, 39287))
parents

##       [,1]  [,2]
## [1,] 36451 39287
## [2,] 41234 40360
## [3,]  9656  9541
## [4,] 32141 39287

# find the FS individuals by looking in the pedigree for their mother and father:
FS <- list()
for (p in seq_len(nrow(parents))) {
  FS[[p]] <- 
    demo_ped$genotype[!is.na(demo_ped$mother) & demo_ped$mother==parents[p, 1] &
                      !is.na(demo_ped$father) & demo_ped$father==parents[p, 2]]
}
# sizes of the FS families:
sapply(FS, length)

## [1] 405  53  76  37

Overviews by FS family

The results of the haplotyping are now in a list named results. We first see how well the different FS families have been haplotyped, using function overviewByFS:

ovwFS <- overviewByFS(haploblock=hblist, parents=parents, FS=FS,
                      hapresults=results)
# The full table is too wide to show;
# for FS family 1 the results are:
knitr::kable(ovwFS$ovw[, 1: 8], digits=c(0,0,0,0,3,0,0,0))

# the final columns for the non-FS individuals and all individuals:
knitr::kable(ovwFS$ovw[, c(1:2, 27:31)])

#

overviewByFS returns a list with two items: ovw and messages. Both are matrices with one row per haploblock. Matrix ovw first gives the number of markers (nmrk) and of assigned haplotypes (nhap) for the haploblock over all individuals, and then for each FS family specific information: parmrk (0, 1 or 2: the number of parents with complete marker dosages), fit (1 if a polysomic segregation could be fitted, else 0), P (the chi-squared P-value of the best fitting segregation of haplotypes, even if this segregation was rejected), mrk (the number of FS individuals with complete marker dosages), imp (the number of FS individuals for which marker dosages were imputed) and hap (the number of FS individuals with an assigned haplotype combination). After the last FS family we have two columns for “rest” (all individuals that are not in the FS families or their parents) and three columns for “all” (all individuals), again with mrk, imp and hap indicating the number of individuals with complete marker data, with imputed marker genotypes and with assigned haplotypes. See the help file for further details (?overviewByFS).

The second item in the OverviewByFS result, messages, is also a matrix with one row per haploblock. It has first a message for the entire haploblock and then one message per FS family.

Pedigree check

Next we check each individual for a match between the possible gametes of its parents and its own haplotype composition, using function pedigreeHapdosCheck:

pedchk <- pedigreeHapCheck(ped=demo_ped, mrkDosage=demo_snpdos,
                           haploblock=hblist,
                           hapresults=results)

## checking haploblock 1 of 7: ctg001
## checking haploblock 2 of 7: ctg002
## checking haploblock 3 of 7: ctg003
## checking haploblock 4 of 7: ctg004
## checking haploblock 5 of 7: ctg005
## checking haploblock 6 of 7: ctg006
## checking haploblock 7 of 7: ctg007

#show part of ped_arr for haploblock 1:
knitr::kable(pedchk$ped_arr[1:8,,1])

#show parents_arr for parents with more than 10 offspring and haploblock 1:
knitr::kable(pedchk$parents_arr[pedchk$parents_arr[,3,1]>10,,1])

#

pedigreeHapdosCheck returns a list with two items, the 3-dim arrays ped_arr and parents_arr. In both cases the first dimension are individuals and the third dimension are haploblocks; the second dimension are the columns as shown above.

Array ped_arr shows for each individual in the pedigree whether it has complete marker data (mrk), whether some of its marker dosages were imputed (imp) and whether it has an assigned haplotype combination (hap), and if its haplotype combination is compatible with that of its parents assuming Double Reduction to be impossible (noDR) or possible (withDR); NA indicates that the individual itself or both its parents have no haplotypes assigned.
Array parents_arr gives information for each individual that is a parent, whether it has complete marker data (par_mrk) and an assigned haplotype combination (par_hap), how many of its progeny have complete marker data (mrk) and haplotype data (hap), and how many of its progeny are compatible with it, assuming Double Reduction to be impossible(nonDRmatch) or possible (DRmatch). For details see the help (?pedigreeHapdosCheck).

Summary statistics

The results of overviewByFS and pedigreeHapdosCheck can be summarized using function calcStatistics:

cst <- calcStatistics(pedchk=pedchk, ovwFS=ovwFS)
knitr::kable(cst$pedstats)

knitr::kable(cst$FSstats, digits=c(0,0,0,2,2,2))

#

calcStatistics returns a list with two matrices. Matrix pedstats gives per haploblock the totals over all individuals from pedchk$ped_arr. The total of columns match.NA + noDR.TRUE + noDR.FALSE and of match.NA + withDR.TRUE + withDR.FALSE is the total number of individuals. Matrix FSstats gives the totals or means per FS family over all haploblocks from ovwFS$ovw. For details see the help (?calcStatistics).

Number of markers vs number of haplotypes

Finally we can get a table of the number of haplotypes vs the number of markers per haploblock:

calcMrkHaptable(ovwFS=ovwFS)

##    haplotypes
## mrk 5 6 7 8
##   4 3 1 0 1
##   5 0 1 1 0

In this example there are 5 haploblocks with 4 markers and 2 haploblocks with 5 markers, so the table contains two rows.

Segregation in one FS, one haploblock

With function showOneFS it is possible to study the segregation of one haploblock in one FS family. Both the segregation of markers and of haplotypes is shown.

# show the segregation for FS number 1 (FSnr=1) and 
# haploblock number 1 (hbresults=results[[1]])
showOneFS(FSnr=1, hbresults=results[[1]], mrkDosage=demo_snpdos, 
          FS=FS, parents=parents)

## $mrkdat
##                (36451) (39287) 43 93 386 436 729 779 1072 <NA>
## frq                386     436 39 65 103 104  47  31    1   15
## WagCmDlf_18464       1       1  0  0   1   1   2   2    3   NA
## WagCmDlf_12697       0       1  0  1   0   1   0   1    0   NA
## WagCmDlf_61726       6       6  6  6   6   6   6   6    6   NA
## WagCmDlf_26276       0       1  0  1   0   1   0   1    0   NA
## 
## $hapdat
##           (36451) (39287) 43 93 386 436 729 779 1072 <NA>
## frq           386     436 39 65 103 104  47  31    1   15
## ctg001_03       5       4  6  5   5   4   4   3    3   NA
## ctg001_08       0       1  0  1   0   1   0   1    0   NA
## ctg001_11       1       1  0  0   1   1   2   2    3   NA
## 
## $usedhap
##                3 8 11
## WagCmDlf_18464 0 0  1
## WagCmDlf_12697 0 1  0
## WagCmDlf_61726 1 1  1
## WagCmDlf_26276 0 1  0

The result of showOneFS is a list with 3 elements. The first two are matrices that show the segregation in terms of marker dosages and haplotype dosages, respectively. The row names of these matrices are first “frq” (frequency, the number of FS individuals with each genotype) followed by the marker names or the haplotype names. The columns of both matrices are a bit complex. The first two columns represent the parents. Their column names are the names of the parents in (brackets), and their first row does not contain a frequency but the marker dosage IDs of the parents. All following columns have as names one of the marker dosage IDs occurring in the FS, and the first row has the number of FS individuals with that marker dosage ID. (A marker dosage ID or mrkdid is a number encoding the dosages of all the SNP markers in an individual). In all columns, the next rows contain the dosages of the markers or the haplotypes. The third element of the result is a matrix listing the composition of the haplotypes present, in terms of their marker alleles (with a 0 being the reference allele and a 1 the alternative allele).