Loading GWAS Data

Loading GWAS data takes the statistics and effects from a GWAS and stores them together. In the process, non-biallelic data is dropped. The two files are described below.

• association file – GWAS is performed using trait data (phenotypes measured on all individuals in an association panel) and genotypic data, usually high density SNP data sets. Following GWAS analysis using the General Linear Model (GLM) or Mixed Linear Model (MLM), TASSEL [2] generates output files presenting associations between the genetic markers used in the study and the trait under study (with correction for population structure and relatedness within the population used). For each SNP-trait association, the F-statistics and p-values are displayed, along with degrees of freedom, error mean square for the model, R^2 of the model (the portion of the total variation explained by the full model), and R^2 of the marker (the portion of total variation explained by the marker but not by the other terms in the model). The p-value and R^2 values are used from every marker-trait association as inputs into PAST. PAST only accepts bi-allelic markers; those with more than 2 alleles are dropped during the analysis.

• effects file - For every marker/trait association in the association file, the number of observations for taxa carrying that allele (Obs), the chromosomal location of the marker, and the estimate of the effect of that allele is calculated for every marker allele and presented in the effects file. Because of the way that TASSEL codes alleles, the last allele estimate for a marker is always zero and the other allele estimates are relative to that.

gwas_data <- load_GWAS_data(demo_association_file, 
                            demo_effects_file)

Loading LD Data

LD data is loaded from the linkage disequilibrium file. In the process, incomplete cases are dropped and the data is split into a data.frame for each chromomsome. Your LD data’s Chromosome information should match your GFF annotation’s chromosome column (the first column). The LD file is described below.

• LD file – Linkage disequilibrium is calculated between pairs of markers that were used in a GWAS study. TASSEL does this and will calculate D’, the standardized disequilibrium coefficient, to determine recombination between pairs of alleles (polymorphic loci). This will give a genetic distance between these pairs; TASSEL also calculates r2 and p values. The genetic distance between markers is not a physical distance, which is generally measured in base pairs or Kbp, and although it is correlated, will vary across chromosomal regions within a species. Thus, associations can be found between a SNP and a causal gene that may be thousands or tens of thousands of base pairs away.

LD <- load_LD(demo_LD_file)

Assigning SNPs to Genes

PAST uses the linkage disequilibrium output from TASSEL between each marker SNP (denoted as the reference SNP) and its closest neighboring SNPs (50 upstream and 50 downstream). Within this window, linkages between SNPs are calculated. The threshold for linkage can be determined from a plot of linkage disequilibrium values (-log(pDiseq) against r2). Based on this plot, Tang et al. [1] defined linkage when the two SNPs being compared had R^2 > 0.8 [3]. PAST uses linkage data to determine which SNP represents the linkage group (the tagSNP), and then uses the tagSNP to determine the linked gene(s) within a window of ± 1Kb. The rationale behind the decision to use 1 Kb is that most genes are regulated within 1 Kb upstream and downstream of the start and stop codons, respectively. At this point, the association and effects data of the tagSNP are transferred to the linked gene. If more than one gene is equally linked to a tagSNP, the attributes of the tagSNP are transferred to both (or all) linked genes.

• genes_file - an annotations file in GFF format. The chromosome in the first column should match the chromosomes in your LD data.

genes <-assign_SNPs_to_genes(gwas_data, 
                             LD, 
                             demo_genes_file,
                             c("gene"),
                             1000, 
                             0.8, 
                             2)

Finding Pathway Significance

Pathway scores are obtained from gene-set enrichment calculations [1, 3]. First, all genes found by tagSNPs are ranked by their effect values from negative to positive in the case of traits where a reduction is beneficial, as in disease resistance, or vice-versa where increase is beneficial, as in the case of yield. Enrichment is based on gene membership in pathways, as assigned by a pathways database supplied by the user. Only pathways with a certain number of genes (supplied by the user) or more genes are considered to reduce bias from small sample size. Next, a running sum is calculated in a manner similar to that used for a weighted Kolmogorov-Smirnov statistic. The running sum statistic increases or decreases if genes are or are not, in the pathway, respectively. The score increases by the fraction of genes in the pathway weighted by the absolute value of the gene effect value or decreases by the fraction of genes not in the pathway. It is a running sum statistic because each gene is considered sequentially in order of their rank among all genes. The final enrichment score (ES) for the pathway is the maximum positive deviation from zero and can be visualized by plotting the values of the running sum statistic against the rank order of genes (see section on rug plots). The significance of a pathway is determined by running 1000 permutations of all genes and their gene effect values to generate a null distribution for the ES. The null distribution mean (μ) and standard deviation (σ) serve to normalize the ES for the pathway. The values of p are then corrected for the false discovery rate as calculated by the QVALUE package [4] in R.

• pathways_file - the pathways file is a tab-delimited file that, for every gene in a pathway contains a line formatted as described below. The gene names in your pathways file must match the gene names in your GFF.

PWY-ID\tPathway Description\tGene

rugplots_data <- find_pathway_significance(genes, 
                                           demo_pathways_file, 
                                           5,
                                           "increasing", 
                                           1000, 
                                           2)

Plotting Selected Pathways

A rug plot is generated for each pathway of interest to visualize the gene-set enrichment calculation. All genes found by the tagSNPs are ranked based on their effect values and their ranks are projected along the x axis of the plot. Hatch marks along the top of the graph denote the rank position of the genes that have membership in the pathway. Values of the running sum statistic enrichment score for each pathway gene are then plotted against their rank. The highest point in the curve is the ES for the pathway and is denoted by the vertical dashed line.

plot_pathways(rugplots_data, 
              "pvalue", 
              0.02, 
              "increasing", 
              tempdir())