# The grammar of population genetics

## Tidy data in population genetics

The fundamental tenet of tidy data that is that each observation should have its own row and each variable its own column, such that each value has its own cell. Applying this logic to population genetic data means that each individual should have its own row, with individual metadata (such as its population, sex, phenotype, etc) as the variables. Genotypes for each locus can also be thought of as variables, however, due to the large number of loci and the restricted values that each genotype can take, it would be very inefficient to store them as individual standard columns.

In `tidypopgen`, we represent data as a `gen_tbl`, a subclass of `tibble` which has two compulsory columns: `id` of the individual (as a `character`, which must be unique for each individual), and `genotypes` (stored in a compressed format as a File-Backed Matrix. The `genotypes` column stores a vector of the row indices of the matrix for each individual and, when printed, shows the first two genotypes for the individual).

The real data reside on disk, and an attribute `bigsnp` of the `genotype` column contains all the information to access it. There is also an additional attribute , `loci` which provides all the information about the loci, including the column indices that represent each locus in the FBM. The vector of row indices and the table of loci can be subsetted and reordered without changing the data on disk; thus, any operation on the `gen_tibble` is fast as it shapes the indices of the genotype matrix rather than the matrix itself.
The `loci` tibble includes columns `big_index` for the index in the FBM, `name` for the locus name (a `character`, which must be unique), `chromosome` for the chromosome (an `integer`, if known, otherwise set to `NA`), `position` for the position on the chromosome (an `integer`, if known, otherwise set to `NA`), `genetic_dist` for the genetic distance on the chromosome (`numeric`, if known, else set to 0), `allele_ref` for the the reference allele (a `character`), and `allele_alt` for the alternate allele (a `character`, which can be `0` for monomorphic loci, following the same convention as plink). Additional individual metadata can be stored as columns in the main `gen_tbl`, whilst additional loci information (such as the position in centimorgans) can be added as columns in the `loci` attribute table.

In principle, it is possible to use use multiple ways to compress the genotypes. `tidypopgen` currently uses a `bigSNP` object from the package `bigsnpr`. It is very fast and well documented, but it is mostly geared towards diploid data. `tidypopgen` expands that object to deal with different levels of ploidy, including multiple ploidy within a single dataset; however, most functions are currently incompatible with ploidy levels other than 2 (but they will return a clear error message and avoid computing anything incorrectly).

## The grammar of population genetics

### The `gen_tibble`

Given information about the individuals, their genotypes, and the loci:

    library(tidypopgen)
    #> Loading required package: dplyr
    #>
    #> Attaching package: 'dplyr'
    #> The following objects are masked from 'package:stats':
    #>
    #>     filter, lag
    #> The following objects are masked from 'package:base':
    #>
    #>     intersect, setdiff, setequal, union
    #> Loading required package: tibble
    example_indiv_meta <- data.frame(
      id = c("a", "b", "c", "d", "e"),
      population = c("pop1", "pop1", "pop2", "pop2", "pop2")
    )
    example_genotypes <- rbind(
      c(1, 1, 0, 1, 1, 0),
      c(2, 0, 0, 0, NA, 0),
      c(1, 2, 0, 0, 1, 1),
      c(0, 2, 0, 1, 2, 1),
      c(1, 1, NA, 2, 1, 0)
    )
    example_loci <- data.frame(
      name = c("rs1", "rs2", "rs3", "rs4", "x1", "x2"),
      chromosome = c(1, 1, 1, 1, 2, 2),
      position = c(3, 5, 65, 343, 23, 456),
      genetic_dist = c(0, 0, 0, 0, 0, 0),
      allele_ref = c("A", "T", "C", "G", "C", "T"),
      allele_alt = c("T", "C", NA, "C", "G", "A")
    )

We can create a simple `gen_tibble` object (of class `gen_tbl`) with:

    example_gt <- gen_tibble(example_genotypes,
      indiv_meta = example_indiv_meta,
      loci = example_loci,
      backingfile = tempfile()
    )
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/file15ba95c4a0906.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/file15ba95c4a0906.rds
    #> with backing file: /tmp/Rtmp6pqvt1/file15ba95c4a0906.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/file15ba95c4a0906.gt')

We are provided information on where the three files underlying the genotype information are stored. As we don’t want to keep the files, we used the tmp directory; normally you will want to use your working directory so that the files will not be cleared by R at the end of the session. It is important that these files are not deleted or moved, as `gen_tibble` stores their paths for future use.

Now let’s have a look at our `gen_tibble`:

    example_gt
    #> # A gen_tibble: 6 loci
    #> # A tibble:     5 × 3
    #>   id    population  genotypes
    #>   <chr> <chr>      <vctr_SNP>
    #> 1 a     pop1        [1,1,...]
    #> 2 b     pop1        [2,0,...]
    #> 3 c     pop2        [1,2,...]
    #> 4 d     pop2        [0,2,...]
    #> 5 e     pop2        [1,1,...]

As discussed above, when this tibble is called, the `genotypes` column prints the first two genotypes of the individual. The `genotypes` column itself contains the indices of the individuals in the FBM as values, and the FBM as an attribute.

To retrieve the rest of the genotypes (which are compressed in the FBM), we use:

    example_gt %>% show_genotypes()
    #>      [,1] [,2] [,3] [,4] [,5] [,6]
    #> [1,]    1    1    0    1    1    0
    #> [2,]    2    0    0    0   NA    0
    #> [3,]    1    2    0    0    1    1
    #> [4,]    0    2    0    1    2    1
    #> [5,]    1    1   NA    2    1    0

If we want to extract the information about the loci for which we have genotypes (which are stored as an attribute of that column), we say:

    example_gt %>% show_loci()
    #> # A tibble: 6 × 8
    #>   big_index name  chromosome position genetic_dist allele_ref allele_alt chr_int
    #>       <int> <chr> <chr>         <dbl>        <dbl> <chr>      <chr>        <int>
    #> 1         1 rs1   1                 3            0 A          T                1
    #> 2         2 rs2   1                 5            0 T          C                1
    #> 3         3 rs3   1                65            0 C          <NA>             1
    #> 4         4 rs4   1               343            0 G          C                1
    #> 5         5 x1    2                23            0 C          G                2
    #> 6         6 x2    2               456            0 T          A                2

Note that, if we are passing a `gen_tibble` to a function that works on genotypes, it is generally not necessary to pass the column `genotypes` in the call:

    example_gt %>% indiv_het_obs()
    #> [1] 0.6666667 0.0000000 0.5000000 0.3333333 0.6000000

However, if such a function is used within a `dplyr` verb such as `mutate`, we need to pass the `genotype` column to the function:

    example_gt %>% mutate(het_obs = indiv_het_obs(.data$genotypes))
    #> # A gen_tibble: 6 loci
    #> # A tibble:     5 × 4
    #>   id    population  genotypes het_obs
    #>   <chr> <chr>      <vctr_SNP>   <dbl>
    #> 1 a     pop1        [1,1,...]   0.667
    #> 2 b     pop1        [2,0,...]   0
    #> 3 c     pop2        [1,2,...]   0.5
    #> 4 d     pop2        [0,2,...]   0.333
    #> 5 e     pop2        [1,1,...]   0.6

Or, more simply:

    example_gt %>% mutate(het_obs = indiv_het_obs(genotypes))
    #> # A gen_tibble: 6 loci
    #> # A tibble:     5 × 4
    #>   id    population  genotypes het_obs
    #>   <chr> <chr>      <vctr_SNP>   <dbl>
    #> 1 a     pop1        [1,1,...]   0.667
    #> 2 b     pop1        [2,0,...]   0
    #> 3 c     pop2        [1,2,...]   0.5
    #> 4 d     pop2        [0,2,...]   0.333
    #> 5 e     pop2        [1,1,...]   0.6

## Standard `dplyr` verbs to manipulate the tibble

The individual metadata can then be processed with the usual `tidyverse` grammar. So, we can filter individuals by population with

    example_pop2 <- example_gt %>% filter(population == "pop2")
    example_pop2
    #> # A gen_tibble: 6 loci
    #> # A tibble:     3 × 3
    #>   id    population  genotypes
    #>   <chr> <chr>      <vctr_SNP>
    #> 1 c     pop2        [1,2,...]
    #> 2 d     pop2        [0,2,...]
    #> 3 e     pop2        [1,1,...]

There are a number of functions that compute population genetics quantities for each individual, such as individual observed heterozygosity. We can compute them simply with:

    example_gt %>% indiv_het_obs()
    #> [1] 0.6666667 0.0000000 0.5000000 0.3333333 0.6000000

Or after filtering:

    example_gt %>%
      filter(population == "pop2") %>%
      indiv_het_obs()
    #> [1] 0.5000000 0.3333333 0.6000000

We can use `mutate` to add observed heterozygosity as a column to our `gen_tibble` (again, note that functions that work on genotypes don’t need to be passed any arguments if the tibble is passed directly to them, but the column `genotypes` has to be provided when they are used within `dplyr` verbs such as `mutate`):

    example_gt %>% mutate(het_obs = indiv_het_obs(genotypes))
    #> # A gen_tibble: 6 loci
    #> # A tibble:     5 × 4
    #>   id    population  genotypes het_obs
    #>   <chr> <chr>      <vctr_SNP>   <dbl>
    #> 1 a     pop1        [1,1,...]   0.667
    #> 2 b     pop1        [2,0,...]   0
    #> 3 c     pop2        [1,2,...]   0.5
    #> 4 d     pop2        [0,2,...]   0.333
    #> 5 e     pop2        [1,1,...]   0.6

There are a number of functions that estimate quantities at the individual level, and they are prefixed by `indiv_`.

# Using verbs on loci

Since the genotypes of the loci are stored as a compressed list in one column, it is not possible to use standard `dplyr` verbs on them. However, `tidypopgen` provides a number of specialised verbs, postfixed by `_loci`, to manipulate loci.

A key operation on loci is their selection (and removal). The compressed nature of genotypes imposes some constraints on the possible grammar. For selection, there are two verbs: `select_loci` and `select_loci_if`. `select_loci` understands the concise-minilanguage spoken by standard `dplyr::select` that allows to easily refer to variables by their names. However, `select_loci` criteria can not be based on the actual genotypes (e.g. on heterozygosity or missingness). For that, we have to use `select_loci_if`, which can operate on the genotypes but is blind to the names of loci.

Let us start by looking at the loci names in our simple dataset:

    loci_names(example_gt)
    #> [1] "rs1" "rs2" "rs3" "rs4" "x1"  "x2"

We can see that there are two categories of loci, one starting with “rs” and the other with “x”. If we wanted to select only loci that have an “rs” code, we would use:

    example_sub <- example_gt %>% select_loci(starts_with("rs"))
    example_sub
    #> # A gen_tibble: 4 loci
    #> # A tibble:     5 × 3
    #>   id    population  genotypes
    #>   <chr> <chr>      <vctr_SNP>
    #> 1 a     pop1        [1,1,...]
    #> 2 b     pop1        [2,0,...]
    #> 3 c     pop2        [1,2,...]
    #> 4 d     pop2        [0,2,...]
    #> 5 e     pop2        [1,1,...]

This gives us a `gen_tibble` with only 4 loci, as expected. We can confirm that we have the correct loci with:

    loci_names(example_sub)
    #> [1] "rs1" "rs2" "rs3" "rs4"

Let us check that this has indeed impacted the individual heterozygosity

    example_sub %>% indiv_het_obs()
    #> [1] 0.5000000 0.0000000 0.1666667 0.1666667 0.4000000

We can also subset and reorder by passing indices:

    example_gt %>%
      select_loci(c(2, 6, 1)) %>%
      show_loci()
    #> # A tibble: 3 × 8
    #>   big_index name  chromosome position genetic_dist allele_ref allele_alt chr_int
    #>       <int> <chr> <chr>         <dbl>        <dbl> <chr>      <chr>        <int>
    #> 1         2 rs2   1                 5            0 T          C                1
    #> 2         6 x2    2               456            0 T          A                2
    #> 3         1 rs1   1                 3            0 A          T                1

This operation could be helpful when merging datasets that do not fully overlap on their loci (more on that later).

    example_gt %>%
      select_loci(c(2, 6, 1)) %>%
      show_genotypes()
    #>      [,1] [,2] [,3]
    #> [1,]    1    0    1
    #> [2,]    0    0    2
    #> [3,]    2    1    1
    #> [4,]    2    1    0
    #> [5,]    1    0    1

The limit of `select_loci` is that it can not directly summarise the genotypes. We can do that separately and then feed the result as a set of indices. For example, we might want to impose a minimum minor allele frequency. `loci_maf()` allows us to inspect the minimum allele frequencies in a `gen_tibble`:

    example_gt %>% loci_maf()
    #> [1] 0.500 0.400 0.000 0.400 0.375 0.200

We can now create a vector of indices of loci with a minimum allele frequency (MAF) larger than 0.2, and use it to select:

    sel_indices <- which((example_gt %>% loci_maf()) > 0.2)
    example_gt %>%
      select_loci(all_of(sel_indices)) %>%
      show_loci()
    #> # A tibble: 4 × 8
    #>   big_index name  chromosome position genetic_dist allele_ref allele_alt chr_int
    #>       <int> <chr> <chr>         <dbl>        <dbl> <chr>      <chr>        <int>
    #> 1         1 rs1   1                 3            0 A          T                1
    #> 2         2 rs2   1                 5            0 T          C                1
    #> 3         4 rs4   1               343            0 G          C                1
    #> 4         5 x1    2                23            0 C          G                2

Note that passing a variable directly to `select` is deprecated, and so we have to use `all_of` to wrap it.

`select_loci_if` allows us to avoid creating a temporary variable to store indices:

    example_gt_sub <- example_gt %>% select_loci_if(loci_maf(genotypes) > 0.2)
    example_gt_sub %>% show_genotypes()
    #>      [,1] [,2] [,3] [,4]
    #> [1,]    1    1    1    1
    #> [2,]    2    0    0   NA
    #> [3,]    1    2    0    1
    #> [4,]    0    2    1    2
    #> [5,]    1    1    2    1

Note that, as we need to tidy evaluate `loci_maf` within the `select_loci_if` verb, we need to provide it with the column that we want to use (even though it has to be `genotypes`). Also note that, with `select_loci_if`, we cannot reorder the loci.

`select_loci_if` is very flexible; for example, we could filter loci with a MAF greater than 0.2 that are also on chromosome 2.

We can use a similar approach to select only alleles on a given chromosome:

    example_gt %>%
      select_loci_if(
        loci_chromosomes(genotypes) == 2 &
          loci_maf(genotypes) > 0.2
      ) %>%
      show_loci()
    #> # A tibble: 1 × 8
    #>   big_index name  chromosome position genetic_dist allele_ref allele_alt chr_int
    #>       <int> <chr> <chr>         <dbl>        <dbl> <chr>      <chr>        <int>
    #> 1         5 x1    2                23            0 C          G                2

Incidentally, `loci_maf()` is one of several functions that compute quantities by locus; they can be identified as they start with `loci_`.

# Grouping individuals in populations

In population genetics, we are generally interested in computing quantities that describe groups of individuals (i.e. populations). Grouping can be used in a number of ways.

As a starting point, we can group by population and get pop sizes:

    example_gt %>%
      group_by(population) %>%
      tally()
    #> # A tibble: 2 × 2
    #>   population     n
    #>   <chr>      <int>
    #> 1 pop1           2
    #> 2 pop2           3

For functions that return one result per individual (such as `indiv_het_obs` that we used before), we can use `summarise`, which returns a new `tibble` with one line per population. For example, we can count the number of individuals per population, as well as their mean heterozygosity with :

    example_gt %>%
      group_by(population) %>%
      summarise(n = n(), mean_het = mean(indiv_het_obs(genotypes)))
    #> # A tibble: 2 × 3
    #>   population     n mean_het
    #>   <chr>      <int>    <dbl>
    #> 1 pop1           2    0.333
    #> 2 pop2           3    0.478

However, note that this is somewhat inefficient, as computing the pop averages requires multiple access to the data. A more efficient approach is to:

    example_gt %>%
      mutate(het_obs = indiv_het_obs(genotypes)) %>%
      group_by(population) %>%
      summarise(n = n(), mean_het = mean(het_obs))
    #> # A tibble: 2 × 3
    #>   population     n mean_het
    #>   <chr>      <int>    <dbl>
    #> 1 pop1           2    0.333
    #> 2 pop2           3    0.478

In this way, we compute all individual heterozygosities in one go (optimising our file access time), and then generate the population summaries.

When working with loci, you can use `reframe` to apply functions across groups. To use `reframe` correctly, be sure to select the `genotype` column within the loci verbs:

    example_gt %>%
      group_by(population) %>%
      reframe(loci_hwe = loci_hwe(genotypes))
    #> # A tibble: 12 × 2
    #>    population loci_hwe
    #>    <chr>         <dbl>
    #>  1 pop1            0.5
    #>  2 pop1            0.5
    #>  3 pop1            0.5
    #>  4 pop1            0.5
    #>  5 pop1            0.5
    #>  6 pop1            0.5
    #>  7 pop2            0.6
    #>  8 pop2            0.5
    #>  9 pop2            0.5
    #> 10 pop2            0.7
    #> 11 pop2            0.6
    #> 12 pop2            0.6

Some functions, such as `loci_maf()`, also have a method for grouped tibbles that allows an even easier syntax:

    example_gt %>%
      group_by(population) %>%
      loci_maf()

These methods are coded in C++ to be very efficient, and they should be preferred to the `reframe()` option. By default, their outputs will return a tidy tibble, with one row per locus per population (i.e. the same as we would get with `reframe()`). However, using the `type` argument, it is also possible to return a list:

    example_gt %>%
      group_by(population) %>%
      loci_maf(type = "list")

Or a matrix:

    example_gt %>%
      group_by(population) %>%
      loci_maf(type = "matrix")

Certain metrics and analyses are naturally defined by a grouped tibble as they refer to populations, such as population-specific Fst. For example:

    example_gt %>%
      group_by(population) %>%
      pop_fst()
    #>        pop1        pop2
    #>  0.34146341 -0.04878049

These type of functions are prefixed with `pop`.

# Functions applies to all pairwise and nwise combinations of individuals or populations

The final group of verbs are prefixed by `pairwise_` and `nwise_`, and they are designed to compute pairwise statistics between all pairs or combinations of n individuals or populations (when targeting populations, we start with `pairwise_pop_` and `nwise_pop_`). For example, to get the Identity by State for all pairs of individuals, we use:

    example_gt %>%
      pairwise_ibs()
    #> # A tibble: 10 × 3
    #>    item1 item2 value
    #>    <chr> <chr> <dbl>
    #>  1 a     b     0.7
    #>  2 a     c     0.75
    #>  3 a     d     0.667
    #>  4 a     e     0.9
    #>  5 b     c     0.6
    #>  6 b     d     0.4
    #>  7 b     e     0.5
    #>  8 c     d     0.75
    #>  9 c     e     0.6
    #> 10 d     e     0.5

# Saving and reading data

We can save a `gen_tibble` with `gt_save()`. This command will save a file with extension `.gt`. Together with the `.rds` and `.bk` files, the `.gt` file include all the information stored in the `gen_tibble`. Note that, whilst the `.rds` and `.bk` file have to share name, the `.gt` file can be named differently (but, by default, if no specific name is given, `gt_save` will use the same pattern as for the `.rds` and `.bk` file).

So, let us save our file:

    gt_file_name <- gt_save(example_gt)
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/file15ba95c4a0906.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/file15ba95c4a0906.rds
    #> with backing file: /tmp/Rtmp6pqvt1/file15ba95c4a0906.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/file15ba95c4a0906.gt')
    gt_file_name
    #> [1] "/tmp/Rtmp6pqvt1/file15ba95c4a0906.gt"
    #> [2] "/tmp/Rtmp6pqvt1/file15ba95c4a0906.rds"
    #> [3] "/tmp/Rtmp6pqvt1/file15ba95c4a0906.bk"

And if we ever need to retrieve the location of the `.bk` and `.rds` files for a gen_tibble, we can use:

    gt_get_file_names(example_gt)
    #> [1] "/tmp/Rtmp6pqvt1/file15ba95c4a0906.rds"
    #> [2] "/tmp/Rtmp6pqvt1/file15ba95c4a0906.bk"

In a later session, we could reload the data with:

    new_example_gt <- gt_load(gt_file_name[1])
    new_example_gt %>% show_genotypes()
    #>      [,1] [,2] [,3] [,4] [,5] [,6]
    #> [1,]    1    1    0    1    1    0
    #> [2,]    2    0    0    0   NA    0
    #> [3,]    1    2    0    0    1    1
    #> [4,]    0    2    0    1    2    1
    #> [5,]    1    1   NA    2    1    0

We can see that our genotypes were recovered correctly.

As we saw at the beginning of this vignette, it is possible to create a `gen_tibble` with data in data.frames and tibbles. We can use that function to wrangle small datasets in custom formats, but more commonly SNP data are stored as PLINK bed files or VCF files. `gen_tibble` can directly read both types of files (including gzipped vcf files), we just need to provide the path to file as the first argument of `gen_tibble`; for example, if we want to read a PLINK bed file, we can simply use:

    bed_path_pop_a <- system.file("extdata/pop_a.bed", package = "tidypopgen")
    pop_a_gt <- gen_tibble(bed_path_pop_a, backingfile = tempfile("pop_a_"))
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/pop_a_15ba93f240620.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/pop_a_15ba93f240620.rds
    #> with backing file: /tmp/Rtmp6pqvt1/pop_a_15ba93f240620.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/pop_a_15ba93f240620.gt')

For this vignette, we don’t want to keep files, so we are using again a temporary path for the backing files, but in normal instances, we can simply omit the `backingfile` parameter, and the `.rds` and `.bk` file will be saved with the same name and path as the original `.bed` file.

We can also export data into various formats with the family of functions `gt_as_*()`. Some functions, such as `gt_as_hierfstat()`, `gt_as_genind()` or `gt_as_genlight()` return an object of the appropriate type; other functions, such as `gt_as_plink()` or `gt_as_geno_lea()` write a file in the appropriate format, and return the name of that file on completion. For example, to export to a PLINK .bed file, we simply use:

    gt_as_plink(example_gt, file = tempfile("new_bed_"))
    #> [1] "/tmp/Rtmp6pqvt1/new_bed_15ba922f5359d.bed"

This will also write a .bim and .fam file and save them together with the .bed file. Note that, from the main tibble, only `id`, `population` and `sex` will be preserved in the .fam file. It is also possible to write .ped and .raw files by specifying `type="ped"` or `type="raw"` in `gt_as_plink()` (see the help page for `gt_as_plink()` for details).

# Merging data

Merging data from different sources is a common problem, especially in human population genetics where there is a wealth of SNP chips available. In `tidypopgen`, merging is enacted with an `rbind` operation between `gen_tibbles`. If the datasets have the same loci, then the merge is trivial. If not, then it is necessary to subset to the same loci, and ensure that the data are coded with the same reference and alternate alleles (or swap if needed). Additionally, if data come from SNP chips, there is the added complication that the strand is not always consistent, so it might also be necessary to flip strand (in that case, ambiguous SNPs have to be filtered). The `rbind` method for `gen_tibbles` has a number of parameters that allow us to control the behaviour of the merge.

Let us start by bringing in two sample datasets (note that we use tempfiles to store the data; in real applications, we will usually avoid defining a backingfile and let the function create backing files where the bed file is stored):

    bed_path_pop_a <- system.file("extdata/pop_a.bed", package = "tidypopgen")
    bigsnp_path_a <-
      bigsnpr::snp_readBed(bed_path_pop_a, backingfile = tempfile("pop_a_"))
    pop_a_gt <- gen_tibble(bigsnp_path_a)
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/pop_a_15ba96667166c.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/pop_a_15ba96667166c.rds
    #> with backing file: /tmp/Rtmp6pqvt1/pop_a_15ba96667166c.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/pop_a_15ba96667166c.gt')
    bed_path_pop_b <- system.file("extdata/pop_b.bed", package = "tidypopgen")
    bigsnp_path_b <-
      bigsnpr::snp_readBed(bed_path_pop_b, backingfile = tempfile("pop_b_"))
    pop_b_gt <- gen_tibble(bigsnp_path_b)
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/pop_b_15ba96ee7d229.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/pop_b_15ba96ee7d229.rds
    #> with backing file: /tmp/Rtmp6pqvt1/pop_b_15ba96ee7d229.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/pop_b_15ba96ee7d229.gt')

And inspect them:

    pop_a_gt
    #> # A gen_tibble: 16 loci
    #> # A tibble:     5 × 3
    #>   id    population  genotypes
    #>   <chr> <chr>      <vctr_SNP>
    #> 1 GRC24 pop_a       [0,0,...]
    #> 2 GRC25 pop_a       [0,0,...]
    #> 3 GRC26 pop_a       [0,1,...]
    #> 4 GRC27 pop_a       [0,0,...]
    #> 5 GRC28 pop_a       [0,0,...]

And the other one:

    pop_b_gt
    #> # A gen_tibble: 17 loci
    #> # A tibble:     3 × 3
    #>   id     population  genotypes
    #>   <chr>  <chr>      <vctr_SNP>
    #> 1 SL088  pop_b       [0,0,...]
    #> 2 SL1329 pop_b       [1,1,...]
    #> 3 SL1108 pop_b       [0,1,...]

Here we are using very small datasets, but in real life, `rbind` operations are very demanding. Before performing such an operation, we can run `rbind_dry_run`:

    report <- rbind_dry_run(pop_a_gt, pop_b_gt, flip_strand = TRUE)
    #> harmonising loci between two datasets
    #> flip_strand =  TRUE  ; remove_ambiguous =  TRUE
    #> -----------------------------
    #> dataset: reference
    #> number of SNPs: 16 reduced to 12
    #> ( 2 are ambiguous, of which 2  were removed)
    #> -----------------------------
    #> dataset: target
    #> number of SNPs: 17 reduced to 12
    #> ( 5 were flipped to match the reference set)
    #> ( 2 are ambiguous, of which 2 were removed)

Note that, by default, `rbind` will NOT flip strand or remove ambiguous SNPs (as they are only relevant when merging different SNP chips), you need to set `flip_strand` to TRUE.

The report object contains details about why each locus was either kept or removed, but usually the report is sufficient to make decisions on whether we want to go ahead. If we are happy with the likely outcome, we can proceed with the `rbind`. Note that the data will be saved to disk. We can either provide a path and prefix, to which ‘.RDS’ and ‘.bk’ will be appended for the `bigSNP` file and its backing file; or let the function save the files in the same path as the original backing file of the first object).

NOTE: In this vignette, we save to the temporary directory, but in real life you want to save in a directory where you will be able to retrieve the file at a later date!!!

    # #create merge
    merged_gt <- rbind(pop_a_gt, pop_b_gt,
      flip_strand = TRUE,
      backingfile = file.path(tempdir(), "gt_merged")
    )
    #> harmonising loci between two datasets
    #> flip_strand =  TRUE  ; remove_ambiguous =  TRUE
    #> -----------------------------
    #> dataset: reference
    #> number of SNPs: 16 reduced to 12
    #> ( 2 are ambiguous, of which 2  were removed)
    #> -----------------------------
    #> dataset: target
    #> number of SNPs: 17 reduced to 12
    #> ( 5 were flipped to match the reference set)
    #> ( 2 are ambiguous, of which 2 were removed)
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/gt_merged.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/gt_merged.rds
    #> with backing file: /tmp/Rtmp6pqvt1/gt_merged.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/gt_merged.gt')

Let’s check the resulting `gen_tibble`:

    merged_gt
    #> # A gen_tibble: 12 loci
    #> # A tibble:     8 × 3
    #>   id     population  genotypes
    #>   <chr>  <chr>      <vctr_SNP>
    #> 1 GRC24  pop_a       [0,0,...]
    #> 2 GRC25  pop_a       [0,0,...]
    #> 3 GRC26  pop_a       [1,1,...]
    #> 4 GRC27  pop_a       [0,0,...]
    #> 5 GRC28  pop_a       [0,1,...]
    #> 6 SL088  pop_b       [0,0,...]
    #> 7 SL1329 pop_b       [1,1,...]
    #> 8 SL1108 pop_b       [0,1,...]

We can look at the subsetted loci (note that we used the first population as reference to determine the strand and order of alleles):

    merged_gt %>% show_loci()
    #> # A tibble: 12 × 8
    #>    big_index name       chromosome  position genetic_dist allele_ref allele_alt
    #>        <int> <chr>      <chr>          <int>        <int> <chr>      <chr>
    #>  1         1 rs3094315  1             752566            0 A          G
    #>  2         2 rs3131972  1             752721            0 G          A
    #>  3         3 rs12124819 1             776546            0 A          G
    #>  4         4 rs11240777 1             798959            0 G          A
    #>  5         5 rs1110052  1             873558            0 G          T
    #>  6         6 rs6657048  1             957640            0 C          T
    #>  7         7 rs2488991  1             994391            0 T          G
    #>  8         8 rs2862633  2           61974443            0 G          A
    #>  9         9 rs28569024 2          139008811            0 T          C
    #> 10        10 rs10106770 2          235832763            0 G          A
    #> 11        11 rs11942835 3          155913651            0 T          C
    #> 12        12 rs5945676  23          51433071            0 T          G
    #> # ℹ 1 more variable: chr_int <int>

Note that the `big_index` values have changed compared to the original files, as we generated a new FBM with the merged data.

By default, `rbind` using the loci names to match the two datasets. The option ‘use_position = TRUE’ forces matching by chromosome and position; note that a common source of problems when merging is the coding of chromosomes (e.g. a dataset has ‘chromosome1’ and the other ‘chr1’). If using ‘use_position = TRUE’, you should also make sure that the two datasets were aligned to the same reference genome assembly.

# Imputation

Many genetic analysis (e.g. PCA) do not allow missing data. In many software implementations, missing genotypes are imputed on the fly, often in a simplistic manner. In `tidypopgen`, we encourage taking this step explicitly, before running any analysis.

Let us start with a dataset that has some missing genotypes:

    bed_file <- system.file("extdata", "example-missing.bed", package = "bigsnpr")
    missing_gt <- gen_tibble(bed_file, backingfile = tempfile("missing_"))
    #>
    #> gen_tibble saved to /tmp/Rtmp6pqvt1/missing_15ba934032838.gt
    #> using bigSNP file: /tmp/Rtmp6pqvt1/missing_15ba934032838.rds
    #> with backing file: /tmp/Rtmp6pqvt1/missing_15ba934032838.bk
    #> make sure that you do NOT delete those files!
    #> to reload the gen_tibble in another session, use:
    #> gt_load('/tmp/Rtmp6pqvt1/missing_15ba934032838.gt')
    missing_gt
    #> # A gen_tibble: 500 loci
    #> # A tibble:     200 × 3
    #>    id     population  genotypes
    #>    <chr>  <chr>      <vctr_SNP>
    #>  1 ind_1  fam_1       [0,.,...]
    #>  2 ind_2  fam_2       [0,.,...]
    #>  3 ind_3  fam_3       [0,0,...]
    #>  4 ind_4  fam_4       [.,0,...]
    #>  5 ind_5  fam_5       [0,0,...]
    #>  6 ind_6  fam_6       [0,0,...]
    #>  7 ind_7  fam_7       [0,0,...]
    #>  8 ind_8  fam_8       [0,0,...]
    #>  9 ind_9  fam_9       [0,0,...]
    #> 10 ind_10 fam_10      [0,0,...]
    #> # ℹ 190 more rows

We can visualise the amount of missingness with:

    missing_gt %>%
      loci_missingness() %>%
      hist()

If we attempt to run a PCA on this dataset, we get:

    missing_pca <- missing_gt %>% gt_pca_autoSVD()
    #> Error: You can't have missing values in 'G'.

Note that using gt_pca_autoSVD with a small dataset will likely cause an error, see manual page for details.

It is possible to obviate to this problem by filtering loci with missing data, but that might lose a lot of loci. The alternative is to impute the missing the data. `tidypopgen` provides a wrapper for fast imputation available in `bigsnpr`, a simple imputation (`gt_impute_simple`) based on the frequency of the alleles at each locus (by random sampling, or using the mean or mode). These methods are fine to impute a few missing genotypes, but they should not be used for any sophisticated imputation (e.g. of low coverage genomes).

We use the simple approach to fix our dataset:

    missing_gt <- gt_impute_simple(missing_gt, method = "mode")

We can now check that our dataset has indeed been imputed:

    gt_has_imputed(missing_gt)
    #> [1] TRUE

However, note that a `gen_tibble` stores both the raw data and the imputed data. Even after imputation, the imputed data are not used by default:

    gt_uses_imputed(missing_gt)
    #> [1] FALSE

And indeed, if we summarise missingness, we still get:

    missing_gt %>%
      loci_missingness() %>%
      hist()

We can manually force a `gen_tibble` to use the imputed data:

    gt_set_imputed(missing_gt, set = TRUE)
    missing_gt %>%
      loci_missingness() %>%
      hist()

However, this is generally not needed, we can keep our `gen_tibble` set to use the raw data:

    gt_set_imputed(missing_gt, set = FALSE)
    missing_gt %>%
      loci_missingness() %>%
      hist()

And let functions that need imputation use it automatically:

    missing_pca <- missing_gt %>% gt_pca_partialSVD()
    missing_pca
    #>  === PCA of gen_tibble object ===
    #> Method: [1] "partialSVD"
    #>
    #> Call ($call):gt_pca_partialSVD(x = .)
    #>
    #> Eigenvalues ($d):
    #>  146.859 106.219 90.352 80.983 69.332 68.427 ...
    #>
    #> Principal component scores ($u):
    #>  matrix with 200 rows (individuals) and 10 columns (axes)
    #>
    #> Loadings (Principal axes) ($v):
    #>  matrix with 500 rows (SNPs) and 10 columns (axes)

Note that, when the function is finished, the `gen_tibble` is back to using the raw genotypes:

    gt_uses_imputed(missing_gt)
    #> [1] FALSE
    missing_gt %>%
      loci_missingness() %>%
      hist()

More details about PCA and other analysis is found in the vignette on population genetic analysis.