Practical 5: Conditional and Joint Analysis (COJO)

Last updated: 2023-07-25

Checks: 7 0

Knit directory: SISG2023_Association_Mapping/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20230530)

The command set.seed(20230530) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 70aa7ab

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 70aa7ab. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


working directory clean

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/SISGM15_prac5.Rmd) and HTML (docs/SISGM15_prac5.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	70aa7ab	Joelle Mbatchou	2023-07-25	fix paths
html	2b96c60	Joelle Mbatchou	2023-07-25	Build site.
Rmd	63f527f	Joelle Mbatchou	2023-07-25	fix paths
html	c09af89	Joelle Mbatchou	2023-06-28	Build site.
Rmd	5163a3d	Joelle Mbatchou	2023-06-28	update year
Rmd	1f3a4dc	Joelle Mbatchou	2023-06-25	add files for website

This practical runs on the SISG2023 server.

This practical aims at

Familiarizing you with GCTA-COJO (or COJO in short)
Exploring different factors influencing COJO outcomes

The practical is run in R but uses the R function system() to run PLINK and GCTA from the terminal. If you have PLINK and GCTA installed on your own computer then you could also run it locally if your prefer. In that case you’d need to update a few links provided below.

The COJO algorithm is not designed for fine-mapping per se. However, many of the challenges illustrated and discussed in this practical are relevant for any method using external linkage disequilibrium (LD) reference.

Part I: the data

We provide an R code that simulates a 1 Mb long chromosome with M=2000 SNPs organized within 20 LD blocks. Each block contains 100 SNPs, among which 5 causal variants. SNPs within a block are numbered between 1 and 100, such that the squared correlation \(r^2_{i_{k}j_{k}}\) of allele counts at SNP \(i_k\) and \(j_k\) within LD block \(k\) is

\[ r^2_{i_{k}j_{k}} = \rho_k^{2|i_k-j_k|} \]

LD blocks are characterized by the parameters \(\rho^2_k\), which varies from 0.1 (when \(k=1\); low LD locus) to 0.9 (when \(k=20\); high LD locus).

The code below generates the LD correlation structure between SNPs in each block.

Run the following commands

set.seed(28072023)
nblocks <- 20
rhos    <- sqrt(seq(0.1,0.9,len=nblocks))
m       <- 100 # number of SNPs per LD block
mcBlock <- 5   # number of causals LD per block
M       <- m * nblocks 
R       <- matrix(0,nrow=M,ncol=M)
icausal <- c()
for(k in 1:nblocks){
  l <- ((k-1)*m + 1):(k*m);
  R[l,l]  <- outer(1:m,1:m,FUN=function(i,j) rhos[k]^abs(i-j))
  icausal <- c(icausal,sample(l,mcBlock))
}

The figure code and figure below shows the LD correlation matrix for SNPs the 20-th LD block (\(\rho^2_{20}=0.9\))

k=20
l=((k-1)*m + 1):(k*m)
heatmap(R[l,l],Rowv=NA,Colv=NA)
print("Extract of LD structure for 20-th LD block")
print( R[l,l][1:5,1:5] )

Run the following commands

chr     <- 10 # Chromosome number
pos     <- 1234557 + sort(sample(0:1e6,M)) # Random Position for SNPs
a1a2    <- do.call("rbind",lapply(1:M,function(j) sample(c("A","C","G","T"),2))) # alleles
snps    <- paste0("SNP",1:M) # SNP ID
ldblock <- rep(1:nblocks,each=m) # LD block ID
names(ldblock) <- snps

The R code below generates and shows the LD score of each SNP on the chromosome (x-axis: genomic position in Mb; y-axis: LD score)

Cols <- sample(colors(),nblocks)
ldscores <- diag(crossprod(R))
plot(pos/1e6,ldscores,pch=19,col=Cols[ldblock],
     axes=FALSE,xlab="Genomic Position (Mb)",
     ylab="LD scores")
axis(1);axis(2)
legend("topleft",legend=paste0("Block #",1:nblocks),
       box.lty=0,pch=19,cex=0.5,col=Cols)
cat(paste0("mean LD score = ",round(mean(ldscores),3),
           " - SD LD score = ",round(sd(ldscores),3)))

Run the following commands. This is a function to generate genotypes corresponding to the specified LD structure. For simplicity, we simulate all SNPs with an allele frequency equal to 0.5.

library(MASS)
simGeno <- function(R,n){
  z1 <- do.call("cbind",lapply(1:nblocks,function(i){
    l <- ((i-1)*m + 1):(i*m)
    mvrnorm(n=n,mu=rep(0,m),Sigma = R[l,l])
  }))
  z2 <- do.call("cbind",lapply(1:nblocks,function(i){
    l <- ((i-1)*m + 1):(i*m)
    mvrnorm(n=n,mu=rep(0,m),Sigma = R[l,l])
  }))
  x <- (z1>0) + (z2>0)
  return(x)
}

Run the following commands to genarate genotypes and phenotypes of GWAS participants. GWAS sample size is Ngwas=100000

Ngwas <- 5e4

## Simulate genotypes
Xgwas <- simGeno(n = Ngwas,R)

## Simulate phenotype
mc    <- length(icausal) # total number of causal variants
q2    <- 0.01 #variance explained by all SNPs on the chromosome
b     <- rnorm(n=mc,mean=0,sd=sqrt(q2/mc))
g     <- sqrt(2)*c((Xgwas[,icausal]-1)%*%b)
e     <- rnorm(n=Ngwas,mean=0,sd=sqrt(1-q2))
Ygwas <- g + e

## Running GWAS
var_x <- apply(Xgwas,2,var)
beta  <- cov(Xgwas,Ygwas) / var_x # estimated regression coefficients
se    <- sqrt( (var(Ygwas) - beta*beta*var_x)/((Ngwas-2)*var_x) ) # standard errors
pval  <- 2 * pt(q=abs(beta/se),df=Ngwas-2,lower.tail = F) # T-distribution

## GWAS data - COJO format
gwas  <- cbind.data.frame(SNP=snps,A1=a1a2[,1],A2=a1a2[,2],
                          Freq=colMeans(Xgwas)/2,beta=beta,
                          se=se,P=pval,N=Ngwas)
print(head(gwas,3))

# folder where to store the data
# default is ".", i.e. current directory
# this can be changed
datPath <- "." 

write.table(gwas,paste0(datPath,"/GWAS.ma"),
            quote=FALSE,row.names=FALSE,
            col.names=TRUE,sep="\t")
causals <- snps[icausal]
write(causals,paste0(datPath,"/causals.snplist")) ## list of causal SNPs

Run the following commands to simulate a LD reference (i.e., set of genotypes in PLINK format) from the same population.

## Set path for PLINK
plink   <- "/data/SISG2023M15/exe/plink"

## Simulate and write LD ref
simLDref <- function(Nldref){
  Xldref <- simGeno(n = Nldref,R)
  refGeno <- t(sapply(1:M,function(j) {
    c(paste0(a1a2[j,1],"\t",a1a2[j,1]),
      paste0(a1a2[j,1],"\t",a1a2[j,2]),
      paste0(a1a2[j,2],"\t",a1a2[j,2]))
  }))
  ped <- do.call("cbind",lapply(1:M,function(j){
    refGeno[j,1+Xldref[,j]]}
  ))
  ## fam file
  iid    <- paste0("IID",1:Nldref)
  fid    <- iid
  pid    <- rep(0,Nldref)
  mid    <- rep(0,Nldref)
  sex    <- sample(1:2,Nldref,replace=TRUE)
  pheno  <- rep(-9,Nldref)
  fam    <- cbind.data.frame(fid,iid,pid,mid,sex,pheno)
  
  ## ped/geno
  mapData <- cbind.data.frame(chr,snps,0,pos)
  pedData <- cbind.data.frame(fam,ped)
  
  write.table(mapData,paste0(datPath,"/ldRef.map"),
              quote=FALSE,row.names=FALSE,col.names=FALSE,sep="\t")
  write.table(pedData,paste0(datPath,"/ldRef.ped"),
              quote=FALSE,row.names=FALSE,col.names=FALSE,sep="\t")
  system(paste0(plink," --file ldRef --make-bed --out ldRef"))
}
simLDref(Nldref = 5000)

Part II: running COJO

If you have run all the commands above then the following files must be available in your current directory. To check type the following command in the terminal.

ls -lt GWAS.ma
ls -lt ldRef.*

You can now run COJO. Either from the terminal

GCTA=/data/SISG2023M15/exe/gcta-1.94.1
${GCTA} --bfile ldRef --cojo-file GWAS.ma --chr 10 --cojo-slct --cojo-p 2.5e-5 --out test1

or from R (calling terminal using the system() command)

gcta <- "/data/SISG2023M15/exe/gcta-1.94.1"

system(paste0(gcta," --bfile ldRef ",
              "--cojo-file GWAS.ma --chr 10 ",
              "--cojo-slct --cojo-p 2.5e-5 --out test1"))

Question 1. How many SNPs are detected? How many of those are causal SNPs? (Note that you can obtain causal SNPs in your currrent R session as causals = snps[icausal], or in the file named causals.snplist).

Question 2. Regenerate LD reference data with a lower sample size Nldref=2000, 1000 and 500 and rerun 1). What do you observe? Are all LD blocks affected the same?

Question 3. Set the variance explained by SNPs on the chromosome to 3% (q2=0.03) and re-run 1) and 2). What can you conclude regarding the number of SNPs detected and the proportion of non-causal SNPs detected?

Part III: fixing COJO?

There is no simple way to fix the inflation of false positive observed when the LD reference is too small. As a rule of thumb, Yang et al. (GCTA website) recommend using sample sizes of at least 4000. Nevertheless, we observe that using a more stringent threshold for detecting collinearity might help.

Question 4. Set the variance explained by SNPs on the chromosome to 3% (q2=0.03) and the size of the LD reference to 1000. Re-run COJO adding the following flag --cojo-collinear 0.1. Quantify the improvement in the number of false positives.

sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] vctrs_0.6.3       cli_3.6.1         knitr_1.43        rlang_1.1.1      
 [5] xfun_0.39         stringi_1.7.12    promises_1.2.0.1  jsonlite_1.8.7   
 [9] workflowr_1.7.0   glue_1.6.2        rprojroot_2.0.3   git2r_0.32.0     
[13] htmltools_0.5.5   httpuv_1.6.11     sass_0.4.6        fansi_1.0.4      
[17] rmarkdown_2.23    jquerylib_0.1.4   evaluate_0.21     tibble_3.2.1     
[21] fastmap_1.1.1     yaml_2.3.7        lifecycle_1.0.3   whisker_0.4.1    
[25] stringr_1.5.0     compiler_4.3.1    fs_1.6.2          Rcpp_1.0.11      
[29] pkgconfig_2.0.3   rstudioapi_0.15.0 later_1.3.1       digest_0.6.33    
[33] R6_2.5.1          utf8_1.2.3        pillar_1.9.0      magrittr_2.0.3   
[37] bslib_0.5.0       tools_4.3.1       cachem_1.0.8

Practical 5: Conditional and Joint Analysis (COJO)

Summer Institute of Statical Genetics (Module 15)

2023-07-25

Part I: the data

Part II: running COJO

Part III: fixing COJO?