This page describes some basic ways in which the PLINK/SEQ library can be used to analyse GWAS datasets. In comparison to the PLINK package, advantages include:

  • the ability to handle larger datasets (e.g. for millions of SNPs directly genotyped)
  • a more flexible analysis of post-imputation, probabilistically-called genotypes
  • that it is much easier to extend functionality via R

Disadvantages are that it:

  • is potentially quite a bit slower than PLINK
  • does not currently contain a similarly large range of GWAS functions

In other words, if PLINK or a similar package works for your GWAS data, use it. Otherwise, you might want to consider using the approaches described below

On this page, we 1) introduce the data, 2) show a basic GWAS using PLINK, 3) repeat that analysis in PSEQ and 4), in R using the PLINK/SEQ library, and finally, 5) illustrate how soft-called genotype data can be analysed.

The data

We will use the data from the PLINK resources page. Specifically, download this ZIP archive and extract the files of interest (bolded) as follows:

   creating: example/
  inflating: example/       
  inflating: example/plink.exe       
  inflating: example/       
  inflating: example/pop.cov         
  inflating: example/Haploview.jar   
  inflating: example/wgas1.ped       
  inflating: example/extra.ped       
  inflating: example/gPLINK.jar      
  inflating: example/command-list.txt  

We can convert this PED/MAP file to a PLINK binary PED (BED) file-format (so that we can subsequently load it into a PLINK/SEQ project), by use of the following PLINK command:

plink --file example/wgas1 --make-bed --out example/wgas2

Some key lines from the LOG file are shown below:

 228694 (of 228694) markers to be included from [ example/ ]
 90 individuals read from [ example/wgas1.ped ] 
 Assuming a disease phenotype (1=unaff, 2=aff, 0=miss)
 49 cases, 41 controls and 0 missing
 45 males, 45 females, and 0 of unspecified sex
 Writing pedigree information to [ example/wgas2.fam ] 
 Writing map (extended format) information to [ example/wgas2.bim ] 
 Writing genotype bitfile to [ example/wgas2.bed ] 

Now that we've created a PLINK binary fileset, we will perform some of the steps of a basic GWAS analysis first using PLINK and then PLINK/SEQ.

Basic GWAS analysis using PLINK

Here, we will perform only two steps: first, to generate MDS components:

plink --bfile example/wgas2 --mds-plot 2 --out example/mds1

Second, to perform a logistic regression of disease on SNP genotype, covarying for the MDS components:

plink --bfile example/wgas2 --logistic --covar example/mds1.mds --covar-name C1-C2 --out example/assoc1

By way of illustrating some of the top results, we note that only two SNPs are significant at a p-value threshold of 1e-4:

awk ' NR == 1 || ( $5 == "ADD" && $9 < 1e-4 ) ' example/assoc1.assoc.logistic
 CHR         SNP         BP   A1   TEST  NMISS       OR     STAT         P 
   8  rs11204005   12895576    A    ADD     90   0.0601   -4.217  2.48e-05
   8   rs2460338   12914531    G    ADD     89  0.06723   -4.179  2.93e-05

The first the SNPs have the following results:

grep -w ADD example/assoc1.assoc.logistic | head -3
 CHR         SNP         BP   A1   TEST  NMISS       OR     STAT         P 
   1   rs3094315     792429    G    ADD     89    2.383    1.253    0.2102
   1   rs6672353     817376    A    ADD     89    3e+08  0.00081    0.9994
   1   rs4040617     819185    G    ADD     90    2.234    1.135    0.2564

Basic GWAS analysis using PSEQ

We can now create a PLINK/SEQ project and load the binary-PED file into it. First, we must create a simple project. In this particular case, we do not link to any existing reference datasets and so the syntax is simply:

pseq proj new-project

Next, we load the GWAS data into the project, using the load-plink command:

pseq proj load-plink --file example/wgas2 --id gwas
 expecting 228694 variants on 90 individuals
 inserted 90 individuals, 228694 variants

By default, the phenotype in the FAM file (that corresponds to the sixth column in the original PED file) is called phe1 and automatically stored in the PLINK/SEQ project's INDDB.

We must also load the covariates (the two MDS components) into the project, conforming to the format described here. An appropriately-formatted file can be downloaded from this link: mds1.phe. Alternatively, it can be generated with the following awk command:

awk ' BEGIN { printf "##mds1,Float,X,\"MDS component 1\"\n"; \ printf "##mds2,Float,X,\"MDS component 2\"\n"; \ printf "#ID\tmds1\tmds2\n" } \ NR>1 { printf $1"\t"$4"\t"$5"\n" } ' example/mds1.mds \ > example/mds1.phe

To load this covariate file into the INDDB:

pseq proj load-pheno --file example/mds1.phe
 Processed 90 rows

To run a logistic regression of disease on allele count and covariates for each SNP, use the glm command. This will adopt a logistic or linear regression, depending on the nature of the phenotype specified.

pseq proj glm --phenotype phe1 --covar mds1 mds2

The ouput from this command lists the variant position and ID, the reference and alternate alleles, the number of non-missing individuals included in analysis, the alternate allele frequency F (for binary traits with no covariates, F_A and F_U represent the case and controls frequencies), then the odds ratio (defined with respect to the alternate allele, e.g. for the first variant below, A is protective), its standard error (on log scale), test statistic and p-value.

VAR		       REF  ALT	  N	 F     OR     SE     STAT       P
chr1:792429:rs3094315    G    A  89  0.876  0.419  0.693   -1.253  0.2102
chr1:817376:rs6672353    A    G  89  0.994  2e-09  24378  -0.0008  0.9994
chr1:819185:rs4040617    G    A  90  0.883  0.447  0.708  -1.1350  0.2564

Note that in this example the designation of reference and alternate are arbitrary, reflecting the annotation in the GWAS dataset of A1 and A2, rather than with respect to any known reference genome sequence. As a consequence, the odds ratios will not necessarily be for the same allele in PSEQ as in PLINK, which only has the concept of the minor allele, rather than reference and non-reference.

Examining the output, it has reproduced the basic logistic regression results seen above: e.g. the same two SNPs ranked highest:

VAR                       REF  ALT   N      F     OR    SE  STAT       P
chr8:12895576:rs11204005    A    G  90  0.522 16.638 0.667 4.217 2.5e-05
chr8:12914531:rs2460338     G    C  89  0.522 14.874 0.646 4.179 2.9e-05

In general, we see as expected the identical set of results being generated by PSEQ as PLINK for the analysis of GWAS data by logistic regression with two covariates.

Basic GWAS analysis using the R interface to PLINK/SEQ

As discussed here, PLINK/SEQ projects can be accessed via R. In R,

library(Rplinkseq) pseq.project("proj")

This assumes that proj is in the working directory. We will use the var.iterate() function to obtain for each marker the odds ratio, standard error and p-value from a logistic regression of disease on allele count, controlling for the two MDS components -- that is, the same analysis as above.

Phenotype and covariate information

The ind.fetch() function can be used to extract phenotype information in the INDDB, for example:

p <- ind.fetch( c("phe1","mds1","mds2") ) str(p)
List of 4
 $ ID  : chr [1:90] "CH18526" "CH18524" "CH18529" "CH18558" ...
 $ phe1: int [1:90] 1 1 1 1 1 1 1 2 2 1 ...
 $ mds1: num [1:90] -0.021 -0.0234 -0.0215 -0.0172 -0.0215 ...
 $ mds2: num [1:90] 0.0054 0.0138 0.0175 0.0206 -0.0108 ...

The full form of the ind.fetch() command is:

 ind.fetch( pheno = { (vector of) phenotype name(s) in INDDB } , 
            mask  = { mask specification string } )

We can confirm that the total number of cases (phe1=2) and controls (=1) is the same as in the original dataset.

table( p$phe1 )
 1  2 
41 49 

We can also confirm that the MDS components look as expected (in this dataset, we expect two clusters, corresponding to Chinese and Japanese HapMap samples):

plot( p$mds1 , p$mds2 )

Defining the analysis function to apply

We need to define a function that takes a variant list object (the format of which is described here). To demonstrate the development of this, we'll use the var.fetch() function to get data on the first variant in the project.

k <- var.fetch( limit = 1 )

This returns a list object:

str( k )
List of 1
 $ :List of 7
  ..$ CHR: int 1
  ..$ BP1: int 792429
  ..$ BP2: int 792429
  ..$ ID : chr "rs3094315"
  ..$ NS : int 1
  ..$ CON:List of 8
  .. ..$ IDX   : int 0
  .. ..$ FSET  : int 0
  .. ..$ REF   : chr "G"
  .. ..$ ALT   : chr "A"
  .. ..$ QUAL  : num -1
  .. ..$ FILTER: chr ""
  .. ..$ META  : list()
  .. ..$ GENO  :List of 1
  .. .. ..$ GT: int [1:90] 2 2 1 2 2 2 1 1 2 2 ...

In this case, the CON element (the consensus sample-variant) contains an element called GENO, which contains genotype information (and any meta-information on genotypes too). In this example, the GENO sub-element contains counts of the number of alternate alleles for each individual (recall that designation as reference and alternate is arbitrary in this particular GWAS-centric example). Thus, the element CON$GENO$GT effectively contains the genotype calls for all 90 individuals:

table( k[[1]]$CON$GENO$GT , useNA="always" )
   1    2 <NA> 
  22   67    1 

In other words, given ALT is A and REF is G in this example, this implies 22 A/G heterozygotes, 67 A/A homozygotes and 1 individual with a missing genotype. So G is the minor allele in this example (although it happens to be the "reference" here).

To perform a logistic regresion for this variant, we can use the glm() function:

summary( glm( p$phe1 == 2 ~ k[[1]]$CON$GENO$GT + p$mds1 + p$mds2 , family="binomial" ) )
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          1.8436     1.2767   1.444    0.149    
k[[1]]$CON$GENO$GT  -0.8682     0.6929  -1.253    0.210    
p$mds1              70.7133    13.6709   5.173 2.31e-07 ***
p$mds2             -14.1353    17.0560  -0.829    0.407    

Here we see that for this first SNP, rs3094315, the association analysis yields similar results, as expected, to the PLINK and PSEQ analyses: an odds ratio of 2.383 (i.e. 1/exp(-0.8682)) for the minor allele, a standard error of 0.693 and a p-value of 0.210.

To apply such an analysis to all variants in a dataset (optionally, after applying various masks), we do not need to load all variants into memory. Rather, here we will define a function to perform the analysis, in this instance that is performed a variant at a time. We will then use the var.iterate() function, which is similar in spirit to the family of apply() functions in R.

Note that var.iterate() does not currently return anything. Rather, it is the job of the user's function to collate any results of interest in one or more global objects. Depending on the nature of the data and the analysis, one may want to adopt different approaches to collecting results (i.e. one could directly write to a file with the sink() command, or use approaches other than rbind() which is unlikely to be very efficient).

f <- function(v) { t <- summary( glm( p$phe1 == 2 ~ v$CON$GENO$GT + p$mds1 + p$mds2 , family="binomial" ) )$coef[c(2,6,14)] ; r <<- rbind( r , as.numeric(t) ) }

We generate an object that will be used to collect the results:

r <- numeric()

In other words, we expect a variant object (called v here), perform the regression described above, and then use a short-cut to return the three quantities of interest (coefficient, SE and p-value, assuming there are 3 predictor variables). Applying this function to the first SNP:

f( k[[1]] ) r
[1] -0.8682405  0.6928636  0.2101624

We now test whether or not the function works, applying it to only the first few SNPs (first clearing the results matrix):

r <- numeric() var.iterate( f , "limit=10" )

We can then observe these results for the first 10 SNPs, in the 10-by-3 matrix that has been populated:

             [,1]         [,2]      [,3]
 [1,]  -0.8682405    0.6928636 0.2101624
 [2,] -14.0002760 1455.3976471 0.9923248
 [3,]  -0.8038097    0.7081790 0.2563596
 [4,]  71.2993081   16.5166310        NA
 [5,]  71.2993081   16.5166310        NA
 [6,]   1.0899151    0.8817309 0.2164182
 [7,]   0.2318913    0.3687724 0.5294672
 [8,]  70.7493405   16.5320981        NA
 [9,]  -0.1884215    1.1037734 0.8644544
[10,]  -0.5227010    0.3989326 0.1901116

To apply this function to all variants,

r <- numeric() var.iterate( f )

We now have results for all SNPs in the r matrix. We can plot the p-values as follows:

plot( -log10( r[,3] ) )

We can also compare these p-values to those generated by PLINK:

d <- read.table("example/assoc1.assoc.logistic" , header=T) d <- d[ d$TEST == "ADD" , ] plot( r[,3] , d$P )

As expected, we can confirm that similar results have been produced by both analysis pipelines. Although it is more cumbersome to use R for this particular analysis, it is easier to use this framework to extend to more complex analyses that aren't necessarily implemented in a genetics package.

Comparison across approaches: speed and filesize

To give a sense of how these different tools perform on this particular task, for this particular struture of dataset:

Tool Size of dataset Time: MAF Time: logistic regression
PLINK, PED file 80M 14.6s 1m15s
PLINK, BED file 12.5M 3.7s 50.8s
PSEQ 33M 14.1s 1m33s
R and PLINK/SEQ 33M 31s 1hr12m

Note: the MAF column means that we simply calculate the MAF for every variant; the logistic function is the analysis described above, including two covariates, for every variant in this example dataset.

Some limited conclusions that can be drawn from this particular comparison: 1) as noted above, PSEQ is generally a little slower than PLINK for basic GWAS analysis, i.e. the price paid for greater scalability and extensibility; 2) the logistic regression analysis in R is particuarly slow, although this reflects the particular glm() function. That is, in this particular example, the price paid for getting data from a PLINK/SEQ project into R, per se is about the same as the cost of initially reading the data into PSEQ -- i.e. the simple MAF function in R only takes twice as long as PSEQ. As mentioned above, for certain purposes, e.g. prototyping new methods, this will not be an issue.

Analysis of imputed GWAS data

In this example, we've directly loaded basic SNP data from a PLINK binary-PED fileset. The glm command can be used with any source of data, however, including VCFs. To base the PSEQ glm tests on soft-called genotypes instead of the hard-calls in the GT field, you can add one of the following options:

--options postprobs=PP

to indicate the presence of a tag in the FORMAT fields called PP that contains genotype posterior probabilities; alternatively,

--options dosages=EC

which can be set to point to dosages (expected counts) in the VCF to be used instead, e.g.

  GT:EC  0/1:1.023  ./.:0.234  0/0:0.001 

Note: These two options not available in current release.