Skip to contents

Read PLINK genotypes into the genotypes= shape

Source: R/io_qc.R
read_plink.Rd

Reads the binary trio prefix.bed / prefix.bim / prefix.fam (SNP-major .bed), or a prefix.raw from plink --recode A. Dosages COUNT THE A1 ALLELE of the .bim (in the .bed encoding: 00 = A1A1 -> 2, 10 = het -> 1, 11 = A2A2 -> 0, 01 = missing -> NA). Stated because half the silent bugs in genomic pipelines are an allele-counting convention nobody wrote down.

Usage

read_plink(
  prefix,
  storage = c("double", "integer", "raw"),
  id = c("iid", "fid_iid")
)

Arguments

prefix

path without extension

storage

type of the dosage matrix: "double" (8 bytes per genotype, NA for missing), "integer" (4 bytes) or "raw" (1 byte, 5 for missing, the BLUPF90 code). Every genotypes= argument takes the three without converting, so "raw" is the one for a large genotyped set

id

"iid" (the default) takes the IID as the animal id and refuses an IID repeated in the file; "fid_iid" takes paste(FID, IID, sep = "_"), for files where the same IID appears in different families

Value

list with ids (character, from the .fam or .raw IID, or FID_IID), m (dosage matrix with marker names), ready for genotypes=, and, from a .bed, map (the .bim as a data.frame chr, snp, cm, pos, a1, a2, one row per column of m; a .raw has no map)

Details

The .fam, .bim and the identifier columns of a .raw are read as TEXT, exactly as written: an IID 0001 stays 0001 (not 1), T and F stay letters (not TRUE/FALSE), NA is the animal named NA, and # or a quote inside an id is part of it. PLINK tells samples apart by FID and IID together, so the same IID in two families is legal there and would be the same animal twice here: it is refused, naming the IID and its families, unless id = "fid_iid" names each sample FID_IID. A line with the wrong number of fields is refused with its line number.

The .bim comes back as map, so the markers of the sex chromosomes can be told apart: PLINK codes a male hemizygous X as homozygous, which makes a sire and his son opposite homozygotes on about 2pq of the X markers without any error. preGSf90 leaves the sex chromosomes out of its parentage test (sex_chr) and recommends dropping them from the evaluation (excludeCHR); the chromosome code of the X depends on the species set PLINK was run with (23 in the human default, X with --output-chr).

Examples

if (FALSE) { # \dontrun{
g <- read_plink("chip")
# the autosomes only, for a species whose X is coded 19
g$m <- g$m[, !g$map$chr %in% c("19", "X")]
} # }