A strong peak can hide other peaks. In this case, the potential causative SNP in the peak can be used as covariate and the GWAS can be run again to assess what is the weight of the other SNPs when the covariate SNP weight is removed from the analysis. To do so one needs to generate a file with 2 columns, the first containing 1s and the second a code for the SNP to use as covariate. For example, if the SNP can be coded as 1 and the reference allele as 0, therefore, a set of 4 accessions were only the 2 first accessions have the SNP would yield a covariate file looking like this:

$ cat covariate_file.txt
1	1
1	1
1	0
1	0

This file can then be used as third argument in run_gwas_gemma.sh such as:

bash run_gwas_gemma.sh phenotype.tsv vcf_file.vcf.gz covariate_file.txt

The covariate analysis allows to assess the weight of the SNP set as covariate. To do so, compare the output in the log files with and without the covariate. One can estimate the weight of the SNP set as covariate by subtracting the PVE (phenotypic variance explained in %) values of the analysis without and with the covariate.

For instance:

Analysis without covariate

## Summary Statistics:
## number of total individuals = 83
## number of analyzed individuals = 83
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs = 954216
## number of analyzed SNPs = 10321
## REMLE log-likelihood in the null model = -175.057
## MLE log-likelihood in the null model = -175.918
## pve estimate in the null model = 0.942693
## se(pve) in the null model = 0.0201995
## vg estimate in the null model = 38.6983
## ve estimate in the null model = 0.906486
## beta estimate in the null model =   8.3849
## se(beta) =   0.104506

Analysis with covariate

## Summary Statistics:
## number of total individuals = 83
## number of analyzed individuals = 83
## number of covariates = 2
## number of phenotypes = 1
## number of total SNPs = 954216
## number of analyzed SNPs = 10315
## REMLE log-likelihood in the null model = -121.042
## MLE log-likelihood in the null model = -123.252
## pve estimate in the null model = 0.193322
## se(pve) in the null model = 0.165238
## vg estimate in the null model = 0.622643
## ve estimate in the null model = 1.00113
## beta estimate in the null model =   11.1001  -8.34676
## se(beta) =   0.157083  0.345247

The PVE explained by the SNP set as covariate is 0.942693 - 0.193322 = 0.749371 (75%). The SNP has therefore a large effect on the phenotype.

One can also estimate directly the effect of the SNP set as covariate by looking at the beta estimate given in the log file of the analysis with covariate. Here beta estimate in the null model = 11.1001 -8.346761 indicates that the alternative allele of my SNP set as covariate reduce my phenotype by -8.34676 (in phenotype unit).

NB: When PLINK is generating the bed, bim, and fam files, it considers the reference allele of each variant as the allele with the highest frequency. This can change the interpretation of the beta in the log file from an analysis using a covariate. The beta gives an intercept and a slop value. For instance, in the example above, the value -8.34 would have actually been 8.34 if my SNP alternative allele had a frequency > 50% in the individuals analyzed. In case of doubt, plot your phenotype by allele and you will clearly see the direction of your phenotype.

One can add several covariates in the covariate_file.txt, as long as the first column is made of 1s. The beta estimate will then include one intercept and a coefficient for each covariate in the same order as indicated in the covariate_file.txt and with the same caveat concerning the definition of the reference allele (see above).

SNP_significant <- GWAS_run(gwas.results, threshold_pvalue = "bonferroni", highlighted_SNP="4:10420088")

SNP_significant <- GWAS_run(path.file, threshold_pvalue = "5E-8", highlighted_SNP="4:10420088")