## Genome-wide association analysis

Now that our data is loaded, filtered, and additional SNP genotypes imputed we are ready to perform genome-wide association analysis. This involves regressing each SNP separately on a given trait, adjusted for sample level clinical, environmental, and demographic factors. Due to the large number of SNPs and the generally uncharacterized relationships to the outcome, a simple single additive model will be employed.

The GWAA function requires two arguments. The genodata argument should specify the entire genotype data object in SnpMatrix format. The phenodata argument should be a data frame with a column of sample IDs, corresponding to the row names of genodata, and a columns for the continuous outcome variable. These columns must be named “id” and “phenotype”, respectively. In order to fit the model, genotype data is converted to numeric format using the as function from snpStats. The genotypes of each SNP are then coded as continuous, thereby taking on the value of 0, 1, and 2. For this example, we wish for the value of the genotype to reflect the number of minor alleles. However, following conversion our values will reflect the opposite. To fix this a flip.matrix procedure is included in our GWAA function, which can be turned on or off using the flip argument.

Due to the large number of models that require fitting, the GWA analysis can be deployed in parallel across multiple processors or machines to reduce the running time. Here we demonstrate two basic methods for performing parallel processing using the doParallel package. This will be carried out differently depending on whether or not the analysis is run on a UNIX based system, though the arguments are the same. The user can specify the number of processes using the worker argument (set to 2 by default). Additional arguments include select.snps and nSplits. The former allows the user to subset the analysis via a vector of SNP IDs. The latter specifies a number of SNP-wise splits that are made to the genotype data. The function runs the GWA analysis on these smaller subsets of the genotype data one at a time. After each subset has finished running the function will print a progress update onto the R console. By default this is set to 10.

### GWAA function

# Genome-wide Association Analysis
# Parallel implementation of linear model fitting on each SNP

GWAA <- function(genodata=genotypes,  phenodata=phenotypes, family = gaussian, filename=NULL,
append=FALSE, workers=getOption("mc.cores",2L), flip=TRUE,
select.snps=NULL, hosts=NULL, nSplits=10)
{
if (!require(doParallel)) { stop("Missing doParallel package") }

#Check that a filename was specified
if(is.null(filename)) stop("Must specify a filename for output.")

#Check that the genotype data is of class 'SnpMatrix'
if( class(genodata)!="SnpMatrix") stop("Genotype data must of class 'SnpMatrix'.")

#Check that there is a variable named 'phenotype' in phenodata table
if( !"phenotype" %in% colnames(phenodata))  stop("Phenotype data must have column named 'phenotype'")

#Check that there is a variable named 'id' in phenodata table
if( !"id" %in% colnames(phenodata)) stop("Phenotype data must have column named 'id'.")

#If a vector of SNPs is given, subset genotype data for these SNPs
if(!is.null(select.snps)) genodata<-genodata[,which(colnames(genodata)%in%select.snps)]

#Check that there are still SNPs in 'SnpMatrix' object
if(ncol(genodata)==0) stop("There are no SNPs in the 'SnpMatrix' object.")

#Print the number of SNPs to be checked
cat(paste(ncol(genodata), " SNPs included in analysis.\n"))

#If append=FALSE than we will overwrite file with column names
if(!isTRUE(append)) {
columns<-c("SNP", "Estimate", "Std.Error", "t-value", "p-value")
write.table(t(columns), filename, row.names=FALSE, col.names=FALSE, quote=FALSE)
}

# Check sample counts
if (nrow(phenodata) != nrow(genodata)) {
warning("Number of samples mismatch.  Using subset found in phenodata.")
}

# Order genodata rows to be the same as phenodata
genodata <- genodata[phenodata$id,] cat(nrow(genodata), "samples included in analysis.\n") # Change which allele is counted (major or minor) flip.matrix<-function(x) { zero2 <- which(x==0) two0 <- which(x==2) x[zero2] <- 2 x[two0] <- 0 return(x) } nSNPs <- ncol(genodata) genosplit <- ceiling(nSNPs/nSplits) # number of SNPs in each subset snp.start <- seq(1, nSNPs, genosplit) # index of first SNP in group snp.stop <- pmin(snp.start+genosplit-1, nSNPs) # index of last SNP in group if (is.null(hosts)) { # On Unix this will use fork and mclapply. On Windows it # will create multiple processes on localhost. cl <- makeCluster(workers) } else { # The listed hosts must be accessible by the current user using # password-less ssh with R installed on all hosts, all # packages installed, and "rscript" is in the default PATH. # See docs for makeCluster() for more information. cl <- makeCluster(hosts, "PSOCK") } show(cl) # report number of workers and type of parallel implementation registerDoParallel(cl) foreach (part=1:nSplits) %do% { # Returns a standar matrix of the alleles encoded as 0, 1 or 2 genoNum <- as(genodata[,snp.start[part]:snp.stop[part]], "numeric") # Flip the numeric values of genotypes to count minor allele if (isTRUE(flip)) genoNum <- flip.matrix(genoNum) # For each SNP, concatenate the genotype column to the # phenodata and fit a generalized linear model rsVec <- colnames(genoNum) res <- foreach(snp.name=rsVec, .combine='rbind') %dopar% { a <- summary(glm(phenotype~ . - id, family=family, data=cbind(phenodata, snp=genoNum[,snp.name]))) a$coefficients['snp',]
}

# write results so far to a file
write.table(cbind(rsVec,res), filename, append=TRUE, quote=FALSE, col.names=FALSE, row.names=FALSE)

cat(sprintf("GWAS SNPs %s-%s (%s%% finished)\n", snp.start[part], snp.stop[part], 100*part/nSplits))
}

stopCluster(cl)

return(print("Done."))
}


### Phenotype data preparation

First we create a data frame of phenotype features that is the concatenation of clinical features and the first ten principal components. The HDL feature is normalized using a rank-based inverse normal transform. We then remove variables that we are not including in the analysis, i.e. HDL(non-normalized), LDL, TG, and CAD. Finally, we remove samples with missing normalized HDL data.

library(GenABEL)
source("GWAA.R")

# Merge clincal data and principal components to create phenotype table
phenoSub <- merge(clinical,pcs)      # data.frame => [ FamID CAD sex age hdl pc1 pc2 ... pc10 ]

# We will do a rank-based inverse normal transformation of hdl
phenoSub$phenotype <- rntransform(phenoSub$hdl, family="gaussian")

# Show that the assumptions of normality met after transformation
par(mfrow=c(1,2))
hist(phenoSub$hdl, main="Histogram of HDL", xlab="HDL") hist(phenoSub$phenotype, main="Histogram of Tranformed HDL", xlab="Transformed HDL")


# Remove unnecessary columns from table
phenoSub$hdl <- NULL phenoSub$ldl <- NULL
phenoSub$tg <- NULL phenoSub$CAD <- NULL

# Rename columns to match names necessary for GWAS() function
phenoSub <- rename(phenoSub, replace=c(FamID="id"))

# Include only subjects with hdl data
phenoSub<-phenoSub[!is.na(phenoSub$phenotype),] # 1309 subjects included with phenotype data print(head(phenoSub))  ## id sex age pc1 pc2 pc3 pc4 ## 2 10004 2 50 -0.012045108 -0.007231015 -0.003001290 -0.0107972693 ## 3 10005 1 55 -0.016702930 -0.005347697 0.014449836 -0.0006151058 ## 4 10007 1 52 -0.009537235 0.004556977 0.002683566 0.0166255657 ## 5 10008 1 58 -0.015392106 -0.002446933 0.020508791 -0.0057241772 ## 6 10009 1 59 -0.015123858 -0.002353917 0.021360452 0.0069156529 ## 7 10010 1 54 -0.012816157 0.005126124 0.014654847 -0.0147082270 ## pc5 pc6 pc7 pc8 pc9 ## 2 -0.0077705400 -0.0046457510 0.0018061075 -0.003087891 -0.001833242 ## 3 0.0345170160 0.0387085513 0.0205790788 -0.012265508 0.003592690 ## 4 -0.0002363142 0.0055146271 0.0159588869 0.027975455 0.029777180 ## 5 -0.0039696226 0.0053542437 -0.0007269312 0.027014714 0.010672162 ## 6 0.0400677558 0.0232224781 0.0152485234 0.013296852 0.022746352 ## 7 -0.0008190769 -0.0003831342 -0.0131606658 -0.013647709 -0.008912913 ## pc10 phenotype ## 2 -0.004538746 -2.2877117 ## 3 -0.002287043 -0.4749316 ## 4 -0.007461255 0.8855512 ## 5 -0.003352997 -0.1644639 ## 6 0.013143889 0.3938940 ## 7 -0.056187339 1.7109552  ### Parallel model fitting Using this phenotype data, we perform model fitting on each of the typed SNPs in thegenotype object and write the results to a .txt file. # Run GWAS analysis # Note: This function writes a file, but does not produce an R object start <- Sys.time() GWAA(genodata=genotype, phenodata=phenoSub, filename=gwaa.fname)  ## Loading required package: doParallel ## Loading required package: foreach ## Loading required package: iterators ## Loading required package: parallel  ## 656890 SNPs included in analysis. ## 1309 samples included in analysis. ## socket cluster with 2 nodes on host 'localhost' ## GWAS SNPs 1-65689 (10% finished) ## GWAS SNPs 65690-131378 (20% finished) ## GWAS SNPs 131379-197067 (30% finished) ## GWAS SNPs 197068-262756 (40% finished) ## GWAS SNPs 262757-328445 (50% finished) ## GWAS SNPs 328446-394134 (60% finished) ## GWAS SNPs 394135-459823 (70% finished) ## GWAS SNPs 459824-525512 (80% finished) ## GWAS SNPs 525513-591201 (90% finished) ## GWAS SNPs 591202-656890 (100% finished) ## [1] "Done."  end <- Sys.time() print(end-start)  ## Time difference of 2.259378 hours  # Add phenosub to saved results save(genotype, genoBim, clinical, pcs, imputed, target, rules, phenoSub, support, file=working.data.fname(7))  ### Model fitting of non-typed SNPs We also perform association testing on additional SNPs from genotype imputation. Here we use thesnp.rhs.tests function from snpStats to perform the analysis based on the imputation “rules” we calculated previously. We need to specify the variables from the phenoSub data frame that we are including in the model with row names corresponding to the sample IDs. The resulting SNPs are combined with the chromosome position information to create a table of SNPs, location and p-value. Finally, we take $$-log_{10}$$ of the p-value for plotting. # Carry out association testing for imputed SNPs using snp.rhs.tests() rownames(phenoSub) <- phenoSub$id

imp <- snp.rhs.tests(phenotype ~ sex + age + pc1 + pc2 + pc3 + pc4 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10,
family = "Gaussian", data = phenoSub, snp.data = target, rules = rules)

# Obtain p values for imputed SNPs by calling methods on the returned GlmTests object.
results <- data.frame(SNP = imp@snp.names, p.value = p.value(imp), stringsAsFactors = FALSE)
results <- results[!is.na(results$p.value),] #Write a file containing the results write.csv(results, impute.out.fname, row.names=FALSE) # Merge imputation testing results with support to obtain coordinates imputeOut<-merge(results, support[, c("SNP", "position")]) imputeOut$chr <- 16

imputeOut$type <- "imputed" # Find the -log_10 of the p-values imputeOut$Neg_logP <- -log10(imputeOut$p.value) # Order by p-value imputeOut <- arrange(imputeOut, p.value) print(head(imputeOut))  ## SNP p.value position chr type Neg_logP ## 1 rs1532624 9.805683e-08 57005479 16 imputed 7.008522 ## 2 rs7205804 9.805683e-08 57004889 16 imputed 7.008522 ## 3 rs12446515 1.430239e-07 56987015 16 imputed 6.844591 ## 4 rs17231506 1.430239e-07 56994528 16 imputed 6.844591 ## 5 rs173539 1.430239e-07 56988044 16 imputed 6.844591 ## 6 rs183130 1.430239e-07 56991363 16 imputed 6.844591  ### Mapping associated SNPs to genes Using a separate data file containing the chromosome and coordinate locations of each protein coding gene, we can locate coincident genes and SNPs. We use the following function to extract the subset of SNPs that are near a gene of interest. # Returns the subset of SNPs that are within extend.boundary of gene # using the coords table of gene locations. map2gene <- function(gene, coords, SNPs, extend.boundary = 5000) { coordsSub <- coords[coords$gene == gene,] #Subset coordinate file for spcified gene

coordsSub$start <- coordsSub$start - extend.boundary # Extend gene boundaries
coordsSub$stop <- coordsSub$stop + extend.boundary

SNPsub <- SNPs[SNPs$position >= coordsSub$start & SNPs$position <= coordsSub$stop &
SNPs$chr == coordsSub$chr,] #Subset for SNPs in gene

return(data.frame(SNPsub, gene = gene, stringsAsFactors = FALSE))
}


The SNP with the lowest p-value in both the typed and imputed SNP analysis lies within the boundaries of the cholesteryl ester transfer protein gene, CETP. We can call the map2gene function for “CETP” to filter the imputed genotypes and extract only those SNPs that are near CETP. This will be used for post-analytic interrogation to follow.

# Read in file containing protein coding genes coords
genes <- read.csv(protein.coding.coords.fname, stringsAsFactors = FALSE)

# Subset for CETP SNPs
impCETP <- map2gene("CETP", coords = genes, SNPs = imputeOut)

# Filter only the imputed CETP SNP genotypes
impCETPgeno <- imputed[, impCETP\$SNP ]

save(genotype, genoBim, clinical, pcs, imputed, target, rules,
phenoSub, support, genes, impCETP, impCETPgeno, imputeOut, file = working.data.fname(8))