8.15 GWAS of one SNP with PLINK

First, we’ll replicate the GWAS that we did in R with just the first SNP of the VCF, rs9699599.

# replicate first SNP association with PLINK
system(command = "./plink --file genotypes --linear --allow-no-sex --snp rs9699599 --pheno GS451_IC50.txt --pheno-name GS451_IC50")

Breakdown of the PLINK command

./plink: Use the PLINK software
--file genotypes: Genotype data files (genotypes.map, genotypes.ped) begin with the string “genotypes”
--linear: Run a linear additive association test for each SNP
--allow-no-sex: Include samples we don’t have sex metadata for
--snp rs9699599: Only run the analysis for a single SNP (rs9699599)
--pheno GS451_IC50.txt: Phenotype data is located in a file called GS451_IC50.txt
--pheno-name GS451_IC50: The column heading of the phenotype to use in the phenotype file is GS451_IC50

After running PLINK, we get an output file called plink.assoc.linear.

Now look at the output of the plink.assoc.linear output file that PLINK produced.

snp1 <- read.table("plink.assoc.linear",
                   header = TRUE)

head(snp1)

##   CHR       SNP     BP A1 TEST NMISS  BETA  STAT       P
## 1   1 rs9699599 558185  G  ADD    85 1.385 2.036 0.04493

How do these results compare to performing the GWAS by hand?

Notice that the beta (i.e., the “slope” or coefficient) and p-value perfectly matches the results we obtained previously with R.