Basic Commands and Options

PLINK is a comprehensive program with an enormous range of functionalities and options. We will introduce some basic commands here to get you started, but inevitably, you will want to visit the PLINK documentation.

Regarding the PLINK version, you may see three versions available currently (Jan 2020). They are: v1.07, v1.9, and v2. PLINK v1.9 and v2 are under active developement whereas v1.07 is not. Throughout this lab, we will use PLINK v1.9. The reason is two-fold: 1) v1.9 has many enhancements over v1.07 and it is still getting even better; 2) As statistical genetics a fast growing field, new features are actively incorporated into newer version.

In the following, we list some commonly used commands in PLINK. To get to know more about how to use PLINK v1.9, check here and here for online documentation.

In addition to the commands which generate files on their own, the following basic options are important. In any analysis, it is important to perform quality control on the input data (so that we reduce the chance to be placed at “garbage in, garbage out” situation). And equally importantly, we should report how the QC is done in the manuscript (e.g. what thresholds/limits are used to filter out outliers, etc). Often, you may want to use many values to explore whether the perceived association or relationship is robust to your a priori limits.

PLINK syntax

PLINK syntax follows basic Unix style commands. However, one notable element of the syntax is that PLINK generally takes the file name without extension. For example, one of the first steps in PLINK is to make a BED file from a PED file. In such an example, a command could be:

plink --file g_data --make-bed --out g_data_out

Here, we are calling the PLINK command and providing the root for the input file “g_data.ped”, commanding PLINK to make a BED file (“–make-bed”), and naming an output “g_data_out.bed”.

PED files

It is helpful to know the correct format for PED files in case you want to troubleshoot or design an automated script to modify an existing PED file. All PED files are white-space (space or tab) delimited files, arranged such that the first six columns are mandatory:

• 1) Family ID
• 2) Individual ID
• 3) Paternal ID
• 4) Maternal ID
• 5) Sex (1 = male, 2 = female; other = unknown)
• 6) Phenotype

All IDs are alphanumeric, and a PED file must have only 1 phenotype in the sixth column, and may be quantitative or qualitative. Every two columns after the first six are genotypes of SNPs listed in .map file in the same order. These SNPs should be biallelic so that can be represented by numbers or letters (1,2,3,4 or A,T,G,C), as long as 0 is not used (this is default for missing data). So, each of the two columns represent the genotype of a biallelic locus. For instance, 7th and 8th column are allele calls for the first variant in the .map file. Therefore, the number of columns in any PED file is equal to 2 times the number of SNPs (genotypic data) plus the leading six columns.

If you’d like to get to know more on PED format, PLINK documentation has detailed description at here. More importantly, there are many more formats that could be input or generated by PLINK. Whenever you encounter a new format, you can get to know about it using PLINK documentation File formats page.

Although BED files (binary PED files) are often used for analyses to reduce computational time, they are much harder to work with since they are in binary format, and thus generally modifications are made to PED files and then they are then converted to BED files using the command in Section 1.2.

Basic Operations (HapMap Example)

curl zzz.bwh.harvard.edu/plink/hapmap1.zip > hapmap1.zip

mkdir plink                   ##Make a new directory for the PLINK tutorial
unzip hapmap1.zip -d plink    ##Unzip the HapMap data into the PLINK directory
cd plink                      ##Go to the PLINK directory

module avail                   ##Step 1
module avail plink             ##Step 2
module load plink/1.90b6.9     ##Step 3

module avail #to list all packages installed on the cluster 
module avail <insert_package_name> #to check if a specific package is installed

plink --file hapmap1 --make-bed --out hapmap1

plink --bfile hapmap1 --missing --out miss_stat

wc miss_stat.imiss

plink --bfile hapmap1 --freq --out freq_stat

plink --file hapmap1 --hardy

plink --file hapmap1 --het

Reading in data files

Although data can be simulated, generated, and analyzed all within R, often we may want to visualize or analyze data from some other source. We can learn R syntax by learning how to import new data. Use the following commands to set the directory and import the downloaded text file lab1_r_example_data.txt as the data object “new_data”. Don’t forget to set your directory to where you saved it:

# To change the working directory, you modify the following line
setwd("data")
# read table
new_data <- read.table("lab1_r_example_data.txt", header=TRUE)

Exploring objects

We can use functions such as “dim” (prints object dimensions) and “class” (prints object class) to investigate attributes of an object. An object’s class is important to know since different functions apply differently to different classes, which may cause errors. Let’s look at a few attributes of new_data, using the following commands:

dim(new_data)

Problem 5
What are the dimensions of new_data?

class(new_data)

Problem 6
What class is new_data?

While dim and class can be used to explore objects, you may want to explore these functions themselves. This is done by this simple syntax:

# One question mark before a phrase opens its specific page of R documentation
?dim

# Two question marks before a phrase searches R documentation
??class

Simulating genotype data

For simple simulation of genotype data, we’ll assume the SNP is biallelic and therefore use a binomial distribution. By passing our desired arguments for the parameters of the rbinom function, we can randomly generate genotypes. We set the parameter \(n\) equal to 1000 for the number of individuals in our simulation. Assuming the individuals are diploid, the number of trials, \(size\), will be 2 per individual. Finally, the probability of “success” for each trial, the minor allele frequency, is \(p\). Run the following simulation of genotype data:

# Simulate random SNP genotypes for 1000 diploid individuals, given a minor allele frequency of 0.2
num_individuals <- 1000
ploidy_level <- 2
maf <- 0.2
geno <- rbinom(n=num_individuals, size=ploidy_level, p=maf)

We now have the object “geno” containing genotypes for 1000 individuals. Given the minor allele is designated as 1 and major allele as 0, each genotype is either homozygous dominant (0+0=0), heterozygous (0+1=1), or homozygous recessive (1+1=2). Let’s rename these as AA, Aa, and aa to make this more intuitive:

# The gsub function acts as a find and replace, type ?gsub in the console for more info
geno <- gsub("0", "AA", geno)
geno <- gsub("1", "Aa", geno)
geno <- gsub("2", "aa", geno)

Check the number of times each genotype occurs by using the the “table” function:

table(geno)

Problem 7
We set our minor allele frequency as 0.2, but the proportion of genotypes containing the minor allele is roughly double that. Why is this the case?