Load GAPIT
library(GAPIT)Before the gwaspr package can be used, we first need to
run GWAS with GAPIT.
Genotype data (myG)
Our population, the LDP, was genotyped using an exome caputure
array and filtered with the following criteria to yield a data set
of 336,367 SNPs, stored in our object myG.
- Only Bi-allelic = TRUE
- Minimum Read Depth = 3
- Minor Allele Frequency = >5%
- Maximum Missing Frequency = 25%
- Heterozygous = <25%
Note:
header = T
# Load our genotype file in hapmap format (note: header = F)
myG <- read.csv("gwaspr_myG_hmp.csv", header = F)
# Map + marker Info
myG[1:10,1:11]## V1 V2 V3 V4 V5 V6 V7 V8
## 1 V1 V2 V3 V4 V5 V6 V7 V8
## 2 rs alleles chrom pos strand assembly center protLSID
## 3 Lcu.1GRN.Chr1p853882 A/G 1 853882 <NA> <NA> <NA> <NA>
## 4 Lcu.1GRN.Chr1p854117 G/A 1 854117 <NA> <NA> <NA> <NA>
## 5 Lcu.1GRN.Chr1p854159 A/C 1 854159 <NA> <NA> <NA> <NA>
## 6 Lcu.1GRN.Chr1p854174 C/G 1 854174 <NA> <NA> <NA> <NA>
## 7 Lcu.1GRN.Chr1p870836 A/G 1 870836 <NA> <NA> <NA> <NA>
## 8 Lcu.1GRN.Chr1p870898 C/T 1 870898 <NA> <NA> <NA> <NA>
## 9 Lcu.1GRN.Chr1p870903 G/A 1 870903 <NA> <NA> <NA> <NA>
## 10 Lcu.1GRN.Chr1p871108 A/G 1 871108 <NA> <NA> <NA> <NA>
## V9 V10 V11
## 1 V9 V10 V11
## 2 assayLSID panel QCcode
## 3 <NA> <NA> <NA>
## 4 <NA> <NA> <NA>
## 5 <NA> <NA> <NA>
## 6 <NA> <NA> <NA>
## 7 <NA> <NA> <NA>
## 8 <NA> <NA> <NA>
## 9 <NA> <NA> <NA>
## 10 <NA> <NA> <NA>
# genotrype calls
myG[1:10,12:17]## V12 V13 V14 V15 V16
## 1 V12 V13 V14 V15 V16
## 2 X3156.11_AGL CDC_Asterix_AGL CDC_Cherie_AGL CDC_Glamis_AGL CDC_Gold_AGL
## 3 G A G A G
## 4 G G G A G
## 5 C C C C C
## 6 G G G G G
## 7 G G N G G
## 8 T C N C T
## 9 G G N A G
## 10 A A A A A
## V17
## 1 V17
## 2 CDC_Greenstar_AGL
## 3 A
## 4 G
## 5 A
## 6 C
## 7 A
## 8 C
## 9 G
## 10 APhenotype data (myY)
# Load our phenotype file
myY <- read.csv("gwaspr_myY.csv")
myY[1:15,]## Name Cotyledon_Color Cotyledon_RedvsYellow DTF_Sask_2017
## 1 CDC_Asterix_AGL Yellow 0 54.7
## 2 CDC_Rosie_AGL Red 1 59.0
## 3 X3156.11_AGL Red 1 60.7
## 4 CDC_Greenstar_AGL Yellow 0 56.7
## 5 CDC_Cherie_AGL Red 1 54.3
## 6 CDC_Glamis_AGL Yellow 0 59.0
## 7 CDC_Gold_AGL Yellow 0 54.0
## 8 CDC_Imax_AGL Red 1 53.0
## 9 CDC_Impower_AGL Yellow 0 57.0
## 10 CDC_KR.1_AGL Red 1 57.0
## 11 CDC_LeMay_AGL Yellow 0 55.7
## 12 CDC_Maxim_AGL Red 1 54.3
## 13 CDC_QG.1_AGL Green NA 57.7
## 14 CDC_Red_Rider_AGL Red 1 57.0
## 15 CDC_Redcoat_AGL Red 1 57.7
## DTF_Nepal_2017
## 1 128.0
## 2 123.3
## 3 125.3
## 4 121.0
## 5 125.3
## 6 123.0
## 7 125.3
## 8 128.7
## 9 123.3
## 10 121.7
## 11 126.7
## 12 125.7
## 13 126.7
## 14 123.5
## 15 126.7
myY <- myY[,-2] Covariate data (myCV)
myCV[1:10,]## Name b
## 1 CDC_Asterix_AGL -1.201707e-03
## 2 CDC_Rosie_AGL -1.752062e-03
## 3 X3156.11_AGL -1.114926e-03
## 4 CDC_Greenstar_AGL -5.607040e-05
## 5 CDC_Cherie_AGL -7.040853e-04
## 6 CDC_Glamis_AGL -1.234765e-03
## 7 CDC_Gold_AGL -9.853332e-05
## 8 CDC_Imax_AGL -1.210779e-03
## 9 CDC_Impower_AGL -1.396606e-03
## 10 CDC_KR.1_AGL -7.986330e-04Run GWAS
Run GWAS on the 3 traits myY
# Run GWAS
myGAPIT <- GAPIT(
# Phenotype data
Y = myY,
# Genotype data
G = myG,
# Set PCA number
PCA.total = 4,
# Select Models
model = c("GLM","MLM","MLMM","FarmCPU","BLINK") )GWAS with covariate
Run GWAS on the DTF_Sask_2017 trait with b
as a covariate.
# Prep data
myY2 <- myY[,c("Name","DTF_Sask_2017")]
# Rename the trait to be rerun with a covariate
colnames(myY2)[2] <- "DTF_Sask_2017_b"
# Run GWAS with a covariate
myGAPIT <- GAPIT(
# Phenotype data
Y = myY2,
# Genotype data
G = myG,
# Covariate data
CV = myCV,
# Set PCA number
PCA.total = 0,
# Select Models
model = c("GLM","MLM","MLMM","FarmCPU","BLINK") )