Quick start example 1: Simulated data

J Huisman

15 Oct 2020

1 Quick start example 1: Simulated data

1.1 Get started

For a succinct version of this example where you can copy/paste all code in one go, see the main vignette

In this example genotypes are simulated from a fictional pedigree, which together with the associated life history data is provided with the package. This pedigree consists of 5 discrete generations with interconnected half-sib clusters (Pedigree II in the paper1).

# Install the package. This is only the required the first time, or if you wish to update
install.packages("sequoia")   

# Load the package. This is required at the start of every new R session.
library(sequoia) 

# Load the example pedigree and life history data
data(Ped_HSg5, LH_HSg5)

# Take a look at the data structure
tail(Ped_HSg5)
head(LH_HSg5)

# or, in Rstudio, view the full dataframe:
View(Ped_HSg5)

1.2 Simulate genotypes

Simulate some genotype data to use for this try-out.

Geno <- SimGeno(Ped = Ped_HSg5, nSnp = 200)

Function SimGeno() can simulate data with a specific genotyping error rate, call rate and proportion of non-genotyped parents, see its helpfile (?SimGeno) for details. Here we use the defaults; 40% of parents are presumed non-genotyped (their genotypes are removed from the dataset before SimGeno() returns its result).

1.3 Parentage assignment

Run sequoia with the genotype data, lifehistory data, and otherwise default values. It is often advisable to first only run parentage assignment, and check if the results are sensible and/or if any parameters need adjusting. Full pedigree reconstruction, including sibship clustering etc., is much more time consuming.

To only run a data check, duplicate check (always performed first) and parentage assignment, specify Module = 'par' (in older versions: MaxSibIter=0):

ParOUT <- sequoia(GenoM = Geno,
                  LifeHistData = LH_HSg5,
                  Module = 'par')
## Genotype matrix looks OK! There are  920  individuals and  200  SNPs.
## Ageprior: Flat 0/1, overlapping generations, MaxAgeParent = 6,6

## 
## ~~~ Duplicate check ~~~
## 
## ~~~ Parentage assignment ~~~
## Parentage ...
## Initial total LL :
## [1] -78021.56
## Post-parentage total LL :
## [1] -60845.1
## Estimating birth years ...
## Calculating parental LLR ...
## assigned 540 dams and 534 sires to 920 individuals

You will see several plots appearing:

In addition, you will see several messages, including about the initial and post-parentage total LL (log10-likelihood). This is the probability of observing the genotype data, given the (current) pedigree; initially it is assumed that all individuals are unrelated. This number is negative, and gets closer to zero when the pedigree explains the data better.

The result is a list, with the following elements:

names(ParOUT)
##  [1] "Specs"         "ErrM"          "AgePriors"     "LifeHist"     
##  [5] "DupLifeHistID" "NoLH"          "PedigreePar"   "TotLikPar"    
##  [9] "AgePriorExtra" "LifeHistPar"

which are explained in detail in the helpfile (?sequoia) and the main vignette.

you find the assigned parents in list element PedigreePar:

tail(ParOUT$PedigreePar)
##         id    dam   sire LLRdam LLRsire LLRpair OHdam OHsire MEpair
## 915 b05187 a04045   <NA>   9.06      NA      NA     0     NA     NA
## 916 a05188 a04045   <NA>   7.07      NA      NA     0     NA     NA
## 917 a05189 a04006 b04177   6.12    6.72    8.33     0      0      0
## 918 b05190 a04006 b04177   6.00    6.19    8.79     0      0      0
## 919 b05191 a04006 b04177   5.16    7.75    8.49     0      0      0
## 920 b05192 a04006 b04177   5.86    6.36    9.21     0      0      0

we can compare these to the true parents, in the original pedigree from which we simulated the genotype data:

chk_par <- PedCompare(Ped1 = Ped_HSg5, Ped2 = ParOUT$PedigreePar)

chk_par$Counts["GG",,]   
##           parent
## class      dam sire
##   Total    540  534
##   Match    540  534
##   Mismatch   0    0
##   P1only     0    0
##   P2only     0    0

Here ‘GG’ stands for Genotyped offspring, Genotyped parent.

1.4 Pedigree reconstruction

Then, we can run full pedigree reconstruction. By default, this re-runs the parentage assignment.

SeqOUT <- sequoia(GenoM = Geno, LifeHistData = LH_HSg5, Module = "ped")

If you don’t want to re-run parentage assignment (e.g. because it took quite a bit of time), or if you have specified several non-default parameter values you want to use again, , you can provide the old output as input.

Re-used will be:

The last is generated by MakeAgePrior(), which has detected that all parent-offspring pairs have an age difference of \(1\), and all siblings an age difference of \(0\), i.e. that generations do not overlap. This information will be used during the full pedigree reconstruction.

ParOUT$AgePriors
##   M P FS MS PS
## 0 0 0  1  1  1
## 1 1 1  0  0  0
## 2 0 0  0  0  0

So run sequoia() with the old output as input,

SeqOUT <- sequoia(GenoM = Geno,
                  SeqList = ParOUT,
                  Module = "ped")
## Genotype matrix looks OK! There are  920  individuals and  200  SNPs.
## using LifeHistData in SeqList
## using AgePriors in SeqList
## using PedigreePar in SeqList
## 
## ~~~ Duplicate check ~~~
## 
## ~~~ Full pedigree reconstruction ~~~
## Sibships - Initial Total LL :
## [1] -60845.1
## Round 01 end, Total LogLik; time (sec):
## [1] -47920.713     13.875
## No. dams, sires for real indiv.:
## [1] 894 895
## Round 02 end, Total LogLik; time (sec):
## [1] -47338.94420      7.59375
## No. dams, sires for real indiv.:
## [1] 895 895
## Round 03 end, Total LogLik; time (sec):
## [1] -47322.34604      1.21875
## No. dams, sires for real indiv.:
## [1] 895 895
## Round 04 end, Total LogLik; time (sec):
## [1] -47322.349678      0.890625
## No. dams, sires for real indiv.:
## [1] 895 895
## Estimating birth years ...
## Calculating parental LLR ...
## assigned 960 dams and 960 sires to 920 + 80 individuals (real + dummy)

Full pedigree reconstruction will take at least a few minutes to run on this fairly simple dataset with one thousand individuals. It may take a few hours on larger or more complex datasets, and/or if there is much ambiguity about relationships due to a low number of SNPs, genotyping errors, and missing data.

You will get several messages:

Again you will see some plots:

For additional plots and a few tables to inspect the pedigree, you can use SummarySeq():

summary_seq1 <- SummarySeq(SeqOUT, Panels=c("sibships", "D.parents", "LLR"))

names(summary_seq1)
## [1] "PedSummary"  "ParentCount" "GPCount"     "SibSize"

And again we can also compare the results to the true parents:

chk <- PedCompare(Ped1 = Ped_HSg5, Ped2 = SeqOUT$Pedigree)

chk$Counts
## , , parent = dam
## 
##     class
## cat  Total Match Mismatch P1only P2only
##   GG   540   540        0      0      0
##   GD   355   355        0      0      0
##   GT   895   895        0      0      0
##   DG    36    36        0      0      0
##   DD    29    29        0      0      0
##   DT    65    65        0      0      0
##   TT   960   960        0      0      0
## 
## , , parent = sire
## 
##     class
## cat  Total Match Mismatch P1only P2only
##   GG   534   534        0      0      0
##   GD   361   361        0      0      0
##   GT   895   895        0      0      0
##   DG    42    42        0      0      0
##   DD    23    23        0      0      0
##   DT    65    65        0      0      0
##   TT   960   960        0      0      0

See ?PedCompare for a more interesting example with some mismatches.

If you wish to count e.g. the number of full sibling pairs that are assigned as full siblings, paternal half siblings, maternal halfsiblings, or unrelated, or similar comparisons, use ComparePairs().

1.5 Save results

Lastly, you often wish to save the results to file. You can do this as an Rdata file, which can contain any number of R objects. You can later retrieve these in R using load, so that you can later resume where you left of. The disadvantage is that you cannot open the Rdata files outside of R (as far as I am aware). Therefore, sequoia also includes a function to write all output to a set of plain text files in a specified folder.

save(SeqOUT, Geno, file="Sequoia_output_date.RData")
writeSeq(SeqList = SeqOUT, GenoM = Geno, folder = "Sequoia-OUT")

  1. Huisman, Jisca. “Pedigree reconstruction from SNP data: parentage assignment, sibship clustering and beyond.” Molecular ecology resources 17.5 (2017): 1009-1024.↩︎