Data formats

Genotype data formats

ADMIXTOOLS 2 can read genotype data in four formats.

  1. Binary PLINK format (PACKEDPED), described here
  2. Binary PLINK 2 format (PFILE), described here
  3. Binary PACKEDANCESTRYMAP format
  4. Text based EIGENSTRAT format

PACKEDANCESTRYMAP and EIGENSTRAT are described here

In all four formats a dataset consists of three files: one file for the genotype matrix, one file for the SNP metadata, and one file for the sample metadata.

The following shows what a data set for three SNPs and five samples grouped into two populations (French.DG, Onge.DG) could look like:

EIGENSTRAT/PACKEDANCESTRYMAP

EIGENSTRAT and PACKEDANCESTRYMAP use the same format for the metadata. The only difference is that the genotype data is stored in a text file in EIGENSTRAT, but in a binary file in PACKEDANCESTRYMAP.

The .ind file has 3 out of the 6 columns in the .fam file, in different order:

S_French-1.DG M French.DG
S_French-2.DG F French.DG
B_French-3.DG M French.DG
BR_Onge-2.DG F Onge.DG
BR_Onge-1.DG F Onge.DG

.snp files are the same as .bim files, except that the first two columns are swapped:

rs3094315 1 0.02013 752566  G   A
rs12124819 1    0.020242    776546  A   G
rs28765502 1    0.022137    832918  T   C

EIGENSTRAT .geno files look like this:

29101
10211
22901

One SNP per row, one sample per column, and missing data is denoted by 9.

Reading genotype files

While there are other R packages for reading and writing binary PLINK files (plink2R, genio), no such packages exist to date for PACKEDANCESTRYMAP format files, so some of the functions shown here may be useful outside of the ADMIXTOOLS framework.

Note that the functions below impose the restriction than individual IDs may not be duplicated across populations or families.

geno = read_plink("plink/prefix")
geno = read_eigenstrat("eigenstrat/prefix")
geno = read_packedancestrymap("packedancestrymap/prefix")

# read SNPs 1000 to 2000
geno = read_packedancestrymap("packedancestrymap/prefix", first = 1000, last = 2000)
# read only ind1 and ind2
geno = read_packedancestrymap("packedancestrymap/prefix", inds = c('ind1', 'ind2'))

For PFILE, there is no read_pfile — PFILE inputs are consumed directly by the f-statistic pipeline (extract_f2, f2_from_geno, f4blockdat_from_geno, qpfstats) which auto-detect the format from the file extensions at the prefix.

Extracting populations

Populations can be extracted from PLINK files (and written to another set of PLINK files) like this:

extract_samples("input/prefix", "output/prefix", pops = c("pop1", "pop2"))

extract_samples can also be used to extract samples, similar to plink --keep.

Combining data sets

ADMIXTOOLS 2 can’t merge data sets. Data in PLINK format can be merged using plink --bmerge.

Computing allele frequencies

afs = plink_to_afs("plink/prefix")
afs = eigenstrat_to_afs("eigenstrat/prefix")
afs = packedancestrymap_to_afs("packedancestrymap/prefix")
afs = pfile_to_afs("plink2/prefix")

# pfile_to_afs supports the same range/subset args as the others
afs = pfile_to_afs("plink2/prefix", inds = c('ind1', 'ind2'))
afs = pfile_to_afs("plink2/prefix", pops = c('pop1', 'pop2'))
afs = pfile_to_afs("plink2/prefix", first = 1000, last = 2000)

# When .pvar lacks centimorgan positions (the default), supply a map
# via cm_file (TSV with variant_id + cm columns)
afs = pfile_to_afs("plink2/prefix", cm_file = "ref_genmap.tsv")

f2 cache metadata

When extract_f2() writes to disk it emits two sidecar files alongside the per-pair .rds arrays:

  • .f2_cache_id — a short text file containing a sha256:<64hex> hash over the inputs and filter args. Cheap to read; useful as a cache key for orchestrators that want to skip re-extraction when nothing has changed.
  • cache_metadata.json — a schema-versioned JSON file with n_snps, n_blocks, pops, build timestamp, filter args, and the same cache_id. Written atomically (tempfile + rename) before .f2_cache_id, so a crash between the two writes leaves a recoverable directory.

Read the metadata programmatically:

meta = read_f2_cache_metadata("my_f2_dir/")
meta$n_snps     # SNPs that actually contributed to f2 blocks
meta$pops       # populations in the cache
meta$cache_id   # matches .f2_cache_id

To check whether an existing cache directory still matches your current inputs without re-extracting, recompute the cache id from the genotype prefix and compare:

fresh  = compute_f2_cache_id("geno_prefix", pops = pops, ...)
cached = read_f2_cache_metadata("my_f2_dir/")$cache_id
identical(fresh, cached)  # TRUE if the cache is still valid for these inputs

If only the genotype prefix and the f2 directory are available (e.g. an orchestrator probe), compute_f2_cache_id("my_f2_dir/") reads the hash directly from the sidecar — no genotype-file access needed.

Admixture graph formats

There are several ways in which admixture graphs can be represented. Some of them are useful because they can be stored as human readable text files, and others are useful because they make it easier for R to process graphs. ADMIXTOOLS 2 has several function for importing, exporting and converting between formats.

Before fitting an admixture graph, the graph only has nodes and edges. After fitting a graph, each edge has a weight associated with it, and maybe some other information. Here I mostly show examples of graphs that are not fitted yet.

Edge list format

The simplest way to represent a graph in a text file is by writing each edge in the graph in a separate line:

R      O
R      n1
n1     n2
n1     n3
n2     A
n2     admix
n3     B
n3     n4
n4     C
n4     admix
admix  D

We can save this in a text file graph.txt and read it into R like this:

graph_el = read_table('~/Downloads/graph.tsv', col_names = F)
plotly_graph(graph_el)

qpgraph() uses this format to represent fitted graphs. In addition to the first two columns (from and to), there will be additional columns to indicate the type of edge, its estimated weight, and possibly more information.

qpgraph(example_f2_blocks, example_igraph)$edges %>% head
## # A tibble: 6 × 6
##   from  to                   type   weight   low  high
##   <chr> <chr>                <chr>   <dbl> <dbl> <dbl>
## 1 R     Chimp.REF            edge  0.0724     NA    NA
## 2 R     N4N                  edge  0.0724     NA    NA
## 3 N0N   Altai_Neanderthal.DG edge  0.00321    NA    NA
## 4 N0N   N3N0                 edge  0.0111     NA    NA
## 5 N1N   N0N                  edge  0.0841     NA    NA
## 6 N1N   N3N1                 edge  0.0362     NA    NA

igraph format

Many ADMIXTOOLS 2 functions use the igraph package, which has its own graph format. The most important ADMIXTOOLS 2 functions work with graphs in either the edge list or the igraph format, but some only work with graphs in the igraph format. It’s easy to convert between the two formats:

class(graph_el)
## [1] "tbl_df"     "tbl"        "data.frame"
graph = edges_to_igraph(graph_el)
class(graph)
## [1] "igraph"
graph_el2 = as.data.frame(igraph::as_edgelist(graph))
class(graph_el2)
## [1] "data.frame"

Original ADMIXTOOLS format

The text representation of the graph above in the original ADMIXTOOLS software would look like this:

vertex              R
vertex              O
vertex              n1
vertex              n2
vertex              n3
vertex              A
vertex              admix
vertex              B
vertex              n4
vertex              C
vertex              D
label               O                   O
label               A                   A
label               B                   B
label               C                   C
label               D                   D
edge                R_O                 R                   O
edge                R_n1                R                   n1
edge                n1_n2               n1                  n2
edge                n1_n3               n1                  n3
edge                n2_A                n2                  A
edge                n3_B                n3                  B
edge                n3_n4               n3                  n4
edge                n4_C                n4                  C
edge                admix_D             admix               D
admix               admix               n2                  n4

If this is stored in a file called graph.txt, it can be imported into R like this:

graph_el = parse_qpgraph_graphfile('graph.txt')

DOT format

DOT format is often used for plotting admixture graphs. You can export a fitted admixture graph like this:

fit = qpgraph(example_f2_blocks, example_igraph)$edges
write_dot(fit, 'graph.dot')

You can read a DOT file back into an igraph with parse_dot():

graph = parse_dot("graph.dot")

parse_dot() is a minimal reader that recovers the topology from from -> to edges. Edge attributes (weights, types) are not preserved on the round trip; for full round-tripping of fitted graphs, use the edge-list format instead.