Performing post-genome-wide association study analysis: overview, challenges and recommendations

Single-variant approach

Genome wide significant threshold. This approach provides the statistics at a SNP level. This approach aims to score the association between SNPs and the target trait. Usually, the score defined by a $\beta$ value or Odd Ratio (OR), and its standard error (SE) gives an idea of the genotype’s effect on the phenotype. Effect strength is computed using OR. Also, scientists use p-value as a measure of how likely this effect is to occur by chance. Statistically, many researchers use p-value to evaluate the hypothesis that there is no statistical evidence for SNP-trait associations. However, this hypothesis will be rejected when p-value is less than a predetermined threshold that is adjusted for a multi-test. Many methods are used for multiple test p-value correction. These methods include Bonferroni correction and false discovery rate (FDR).²⁰ However, some of these methods do not consider linkage disequilibrium (LD) and they may be too stringent. In general, a threshold of $5e‐08$ is acceptable for human association studies.^21,22 This value has been computed using independent signals in genomes. However, this value is not absolute and it depends on many factors, including LD (and diversity), array type, whole genome sequencing or whole exome sequencing, position number in array or sequencing, imputation panels, and positions finally analyzed.²³ For trans-ethnic populations, it is recommended to estimate the threshold based on population diversity and LD. For instance, analyzing 1000 Genomes data, the suggested significance thresholds were $3.24e‐08$ for Africa, $9.26e‐08$ for Europe, $1.83e‐07$ for Mixed America, $1.61e‐07$ for East Asia and $9.46e‐08$ for South Asia.²⁴ Although, a recent study using an African population proposed a threshold of $5e‐09$.²⁵

Suggestive threshold. Some authors consider GWA threshold very stringent, so studies often include a “suggestive threshold”. This is superior to a genome-wide significant threshold, and values found in the literature are $1e‐05$, $1e‐06$ or lower. Studies estimated on “independant LD block” from 1000 Genomes used $1e‐05$ for the Affymetrix $500K$ and Illumina $317K$ GWAS SNP panels, and $1e‐06$ for HapMap CEPH Utah and Yoruba populations.²⁶ The code below contains an R command to select significant SNPs using a cutoff value of $5e‐08$.

## library need : data.table
#open files with fread function
data.gwas<-fread("result.gwas")
# head to obtain the header of gwas file
head(data.gwas)
# selected lines using threshold value of 5E-8
data.gwas[data.gwas$p.value<5*10**-8,]
# obtained information about min pvalue
data.gwas[which.min (data.gwas$p.value),]

Inflation factors. Population and cryptic relatedness can cause spurious associations in GWAS with p-values higher than random leading to false positive signals. Genomic control (GC) approach is extensively used to effectively control false positive signals. Genomic inflation factors ($\lambda$) can be computed as the median of the resulting chi-squared (${\chi }^2$) test statistics divided by the expected median of the chi-squared distribution.²⁷ Refer to equation 1 below

Z ($\mathit{\text{Beta}}/\mathit{Se}$), ${\chi }^2$ and p-value can be used to compute inflation factor, using $Z^2$ for $Z$, quantile of 1 - p-value at 1 degrees of freedom. The code below demonstrates how to compute inflation factors using the build-in R function $\mathit{\text{qchisq}}$ that can be used to calculate value of quantile for a ${\chi }^2$ distribution.

# compute inflation factors
## use quantile function of chisq to p.value
data.gwas$p.value.qchisq <- qchisq(data.gwas$p.value, 1, lower.tail=FALSE)
## computed lambda
median(data, na.rm=TRUE)/qchisq(0.5, df)

Global visualisation of results. A common way to visualize GWAS results is the Manhattan plot. The Manhattan plot demonstrates the physical location of SNPs distributed by chromosome in the x-axis with the degree to which a SNP is associated with the target trait, i.e. p-value scaled by log10 in the y-axis. A Manhattan plot can be done in R software using the qqman package, which includes functions for creating Manhattan plots and q-q plots from GWAS results.²⁸ The code below demonstrates how to create Manhattan plots, and q-q plots the qqman R package.

## do a qqplot and Manhattan plot using library qqman
library('qqman')
## used function qqplot from qqman
qqplot(data.gwas$p.value)
## used function manhattan from qqman
manhattan(data.gwas, chr="chr", bp="bp", p="p.value",snp="rsid")

Local visualisation of results. Local visualisation can be done using the method described by Pruim et al.²⁹ In this method, information of LD, genes, and previous results of GWAS can be added. Zoom version 2 offers a virtual analysis of local GWAS results. On the other hand, LocusTrack from the UCSC genome-browser adds more annotation compared to LocusZoom.³⁰ Furthermore, BigTop³¹ is capable of providing three-dimensional visualisation of data using allele frequency as a third dimension. The code below demonstrates how to use LocusZoom to visualize SNPs considering LD information.

## use R to reformat your file to be used by locuszoom
R -e "library(data.table);
data.gwas<-fread('result.gwas');
data.gwas<-data.gwas[,c('chr','bp','bp','rsid','af','p.value')];
names(data.gwas)<-c('#CHROM', 'BEGIN', 'END', 'MARKER_ID', 'MAF', 'PVALUE');
write.table (data.gwas, file='gwas.epacts',row.names=F, col.names=T, sep='\t', quote=F)"
## use locus zoom to do a plot,
locuszoom/bin/locuszoom --epacts gwas.epacts \
--delim tab --refsnp rssnp \
--flank 10000 --pop EUR \
--build hg19 --source 1000G_Nov2014\
-p rs7412 --no-date \

Gene-scoring approach

This approach considers the association between a trait and all SNPs within a predefined window around genes rather than each marker individually. In many cases, this approach is more powerful than traditional individual-SNP-based GWAS.

A gene score (GS) is a value given to a gene representing some measures related to a genetic trait. Thus, all the statistical summary values such as p-values and fold changes could be considered gene scoring values. For pGWAS analysis, GS is defined as the sum of all statistically significant alleles (i.e. the risk alleles) of the selected SNPs present in each individual under investigation.³² Various algorithms have been developed to calculate GCs based on GWAS summary statistics.^33-35 It is a common approach to encode and adjust SNPs values prior to calculating GC during analysis.³² The process of SNP encoding aims to ensure that all SNPs are positively correlated with the outcome.³² One more essential aspect to consider while calculating GS is the effect size. Variation in the effect size reflects reduction of predictive power of GS.³² The GS method mostly used in pGWAS analysis is the gene-based p-value.

Gene level p-value. In many pGWAS pipelines, the first step in downstream GWAS analysis is to assign the SNPs to functional genomic features. The latter includes: coding genes, non-coding RNAs, 5’UTR, 3’UTR, proximal promoters, regulatory element, and enhancer elements.

Two common methods for assigning SNPs to their corresponding genes are GLOSSI³⁶ and VEGAS.^34,37 GLOSSI is available as an R package, while VEGAS is available as an online tool as well as a stand-alone tool to be run on a local machine. Besides VEGAS and GLOSSI, many pGWAS tools provide procedures to calculate genes’ p-values. For instance, ancGWAS provides four different methods to calculate the genes p-values: Simes, Smallest, Fisher, and Gwbon. Similar to ancGWAS, MAGMA provides three methods to calculate genes’ p-values. However, one of the MAGMA methods is similar to the VEGAS method. The other two methods are by considering either the smallest SNP p-value, or the highest SNP p-values.

Pathway-sub-network-based approaches

This approach considers the fact that complex biological phenomena addressed by GWAS, including the molecular basis of the rare disease, often arise due to gene interaction rather than single gene effect. Therefore, scientists analyse GWAS based on the biological network theory to understand the disease-causing genes and mechanisms involved in traits and complex diseases, such as in rare diseases. This approach is based on the results obtained from the gene-based association test and provides a higher level of complexity by considering biological networks and/or genes ontology. Analysing GWAS based on biological networks allows us to capture biological interactions between various molecules such as proteins, functional DNA motifs, coding and non-coding RNA, as well as disease mechanisms and to consider epigenetic changes, including methylation states or other modifications (phosphorylation, acetylation, etc). Also, this approach aims to map genes that are associated with significant SNPs into known pathways/Gene Ontology terms. The result of this approach provides information about the over-represented pathways in a given set of genes/SNPs.

Fine-mapping and lead SNPs

GWAS will often identify a number of SNPs in LD with each other as being associated with the phenotype. The lead SNPs are those with the most significant p-value – they may be causal or not. The other SNPs in this region may only have an association because they are in LD with the causal SNP or they may be independently associated. This GWAS limitation may be observed when integrating information of variability in allele frequency, SNPs in LD block and imputation.^38,39 For instance, when simulating GWAS results using OR value of $1.5$ and allele frequency of $0.5$, only $21$% of the simulations demonstrated that the identified causal variants were not the most strongly associated variants.⁴⁰ On the other hand, changing the previous parameters OR value of $1.1$ and allele frequency of $5$%, resulted in $2.4$% concordance between causal variants and lead SNPs.⁴⁰ Therefore, changing the previous parameters OR value of $1.1$ and allele frequency of $5$%, resulted in $2.4$% concordance between causal variants and lead SNPs.⁴⁰

Therefore, fine-mapping methods are used to identify causal variants that are associated with the target trait and the number of putative causal variants from GWAS data. The fine-mapping approach integrates summary statistics from GWAS data, LD and functional annotations. There are two fine-mapping methods: (i) heuristic method that penalizes regression models, and (ii) Bayesian fine-mapping methods.^38,39

Lead SNPs. Lead SNPs are defined as independent SNPs that have reached a minimum p-value threshold, i.e. they are independent of each other at the LD threshold. It is common to measure LD as $r^2$ where the square of the correlation coefficient between any two indicator variables is $r^2$.⁴¹ PLINK has implemented a function called’ld-clump’ that clumps independent SNPs.¹⁵ Moreover, FUMA tool⁴² identifies lead SNPs by double clumping method. The first clumping is used for clumping SNPs with p-value $<$ $0.05$ at genome-wide significant p-value, i.e. p-value $<$ $5e‐08$ and independent at $r^2<$ $0.6$. This first clumping function reports significant independent SNPs. The second clumping is of significant independent SNPs at $r^2<$ 0.1 and reports lead SNPs. The code below demonstrates how to use the Plink tool to perform clumping to get lead SNPs.

## define lead snps using plink
# --clump-p1 Significance threshold for index SNPs
# --clump-p2 Secondary significance threshold for clumped SNPs
# --clump-r2 LD threshold for clumping
# --clump-kb Physical distance threshold for clumping
# --clump-snp-field your snp field must be same than in your assoc file
# --clump-field: p value header
plink --bfile plinkbase --clump result.gwas --clump-snp-field rsid\
--clump-field p.value --clump-p1 0.001 --clump-p2 0.1 \
--clump-r2 0.1 --clump-kb 250 --out assoc

Heuristic fine-mapping approaches. These are used to identify potential causal SNPs.³⁸ They are proposed to filter SNPs around lead SNPs considering the value of their pairwise correlation $r^2$. They consider a hierarchical clustering technique to cluster all SNPs in a given region based on their pairwise $r^2$. Another way is to use LD block by direct extraction of block and the selection of the same position on the block or by the visualisation representation of LD using LocusZoom²⁹ or Haploview.⁴³ Nonetheless, the described method is not a statistical approach to define putative causal variants, such as Bayesian methods or regression models. The code below demonstrates how to combine Plink tool and Haploview to perform heuristic fine-mapping.

### use haploview to plot your data
#### haploview has also an web interface
plink --keep list.ind -bfile $baseplink --recode tab --out fileres \
--chr chr --from-kb begin --to-kb end --maf 0.01 \
# execute Haploview in the previous range
java -jar Haploview.jar -n -minMAF $cut_maf -missingCutoff 0.01 \
-pedfile fileres.ped -map fileres.map -png -blockoutput

Bayesian methods: framework. The Bayesian method is commonly used to identify causal variants in a predefined SNPs window containing $n$ number of SNPs. Knowing data ($D$) the Bayesian method aimed to maximize statistical model ($M$) using the following conditional probability

Users first define an initial number of causal variant $c$ (between $1,m$ SNPs). This number $c$ is defined using genome-wide significant SNPs. To model Bayesian statistics, software will define a model $M$ that contains $c$ SNPs. This model often will restrain choice to the significant SNPs or suggest significant SNPs ($m$) in the windows. Considering the rule of combination $C_{m,c}=\frac{m!}{c!\left(n-c\right)!}$, the probability model $P\left(D|M\right)$ is relatively easier to compute. The following Bayesian rules can be used

Most of the tools are aimed to find a maximum probability to have the best combination of the causal variants. To model Bayesian statistics, many tools integrate various data such as LD and GWAS summary statistics.

Bayesian methods: posterior inclusion probability. This approach is used to compute post inclusion probability (PIP) at each SNP $i$. The PIP is computed by the sum of the posteriors over all models that include SNP $i$ as a causal variant.³⁸

Using the rank of PIP is an excellent way to select putative causal variants.³⁸ The PIP approach should be used with caution to identify causal variants in the high LD regions. Therefore, it is recommended to estimate the posterior expected number of causal SNPs by summing the estimated PIPs for all SNPs in the region.

Bayesian methods: credible sets. This approach is used to define a set of variants that could have a good candidate. One way to estimate a credible set is to use PIP by ranking the values and doing the cumulative sum of PIP from the largest to the smallest. Then select all variants where the sum is less than the predefined $\alpha$ cutoff value. In general, researchers use 99% or 95% as a recommended value for $\alpha$. Studies have shown that credible set and coverage probabilities are over-conservative in most fine-mapping situations as data sets are not randomly selected from among all causal variants. Therefore, an adjusted coverage is proposed to reduce such over-conservation in this approach.⁴⁴ The code below demonstrates how to combine Plink tool and FINEMAP to perform Bayesian-based fine-mapping analysis.

# extract a range of interest
plink -bfile $plk --keep-allele-order --chr chr
--from-kb begin --to-kb end --make-bed -out plink_range
#compute LD
plink --r2 square0 yes-really -bfile plink_range -out ld_range
# format ld
sed 's/\t/ /g' tmp.ld > $outld
## Finemaping
#create a config file
echo "z;ld;snp;config;cred;log;n_samples" \
> fileconfig.finemap
echo "$filez;$ld;${out}.snp;${out}.config;${out}.cred;${out}.log
;${params.n_pop}" >> $fileconfig
FINEMAP --cond --in-files fileconfig  --log \
--cond-pvalue 0.0000001 --n-causal-snps $ncausalsnp
### caviarbf
caviarbf -z ${filez} -r $ld -t 0 -a ${params.caviarbf_avalue} \
-c $ncausalsnp -o ${output} -n ${params.n_pop}
nb=`cat ${filez}|wc -l `
modelsearch -i $output -p 0 -o $output -m \$nb

Bayesian methods: trans ethnic fine-mapping. This is used to perform trans-ethnic fine-mapping studies using simulation results in similar fine-mapping resolution among the European and Asian ancestries. However, the inclusion of samples with African ancestry in meta-analysis leads to a significant improvement in fine-mapping resolution due to the lowest LD in the African ancestry population. The probability that the lead GWAS variants were also the causal variants increased using trans-ethnic GWAS data.^40,45 It is a process that relies on disparate LD patterns in populations of diverse genetic ancestries to localize the causal variants. This approach has been successfully implemented to fine-map and leads to several common GWAS findings.⁴⁶

Integrating annotation into fine-mapping. Functional annotations have been shown to improve the discovery power and fine-mapping accuracy. Therefore, some tools such as CausalDB,⁴⁷ PAINTOR⁴⁸ and BIMBAM⁴⁹ have integrated expression quantitative trait locus (eQTL) information in fine-mapping approach.

Bayesian methods: software. Several Bayesian-based fine-mapping tools have been developed using summary statistics, LD and eQTL (Table 1). Also, some pipelines have integrated different Bayesian-based software in order to compare their results (fine-mapping of h3agwas, fine-mapping in FinnGen, and FM pipeline).

Other fine-mapping approaches. Other approaches, such as regression models, are used with all SNPs in the lead SNPs region to analyze SNPs jointly. The comparison of various approaches including elastic net, ridge, Lasso, MCP, and the normal-exponential shrinkage prior, have shown that penalized methods outperform single marker analysis.⁵³ Furthermore, a forward/stepwise regression can be used to test the independence of multi SNPs using the following algorithm:

The R library SusieR has implemented a method of regression fine-mapping analyses.

Other resources for fine-mapping analyses. Several review articles have been published for fine-mapping analyses. Nevertheless, we recommend the following scientific articles as a good source for beginners: “From genome-wide associations to candidate causal variants by statistical fine-mapping”,³⁸ “A practical view of fine-mapping and gene prioritization in the post-genome-wide association era”,³⁹ and “Fine-mapping genetic associations”.⁵⁴

Conditional association and imputation using summary statistics

The aim of performing conditional association and imputation using summary statistics is to evaluate the association between SNPs and biological trait by combining various GWAS summaries from different studies. This method requires a reference population to estimate LD information. The imputation method performs meta-analysis to infer the missing genotypes among the studies before evaluating the association between the SNPs and the biological trait. This method estimates the effects of many variants that are not directly genotyped.

Polygenic predictions of disease risk

This method is used to predict disease risk using GWAS summaries.⁵⁵ Polygenic risk score (PRS) could be used to predict an individual’s likelihood to develop a specific trait or to estimate the level of predictive power that the trait is associated with a particular set of variants.⁵⁶ Although PRS methods are classified into Bayesian-based methods and non-Bayesian methods, there are more classifications of underlying PRS methods.⁵⁷ PRS is calculated by aggregating effects from a large set of causal SNPs. Several tools were developed to calculate PRS. For performing PRS studies, we highly recommend this recent review paper.⁵⁸ PRS is usually computed after the challenges associated with GWAS are carefully addressed. Here, we demonstrate key quality control (QC) measures and include sample bash scripts. Table 2 summarizes the seven QC measures and contains the guidelines on specific thresholds. Thresholds can differ depending on the study’s unique features.

GWAS effect allele. As datasets for PRS come from different GWAS experiments, it is critical to ensure consistency. Knowing which allele is considered the effect allele, it is vital to get an accurate PRS score. However, the effect allele is not labeled clearly in many datasets.⁵⁹ Different allele coding schemes exist, including Illumina’s TOP/BOTTOM coding concept, ALT/REF, effect/other, HapMap’s forward allele coding, Illumina’s A/B allele coding, Affymetrix’s A/B allele coding, REF (reference)/ALT (alternative), PLINK’s 1/2 allele coding, A1 (allele1)/A2 (allele2), A0 (allele 0)/A1 (allele1), effect allele/non-effect allele, effect allele/other allele, and many others.^58,59 To avoid allele inconsistency, researchers should carefully read the documentation of GWAS datasets.

SNPs level errors. The QC assessment at the SNP level is crucial to avoid misleading PRS. SNPs level errors include (i) mismatching SNPs, i.e., inconsistent SNPs due to position difference in genomic position or nucleotide type, (ii) existence of duplicate SNP, (iii) ambiguous SNPs, i.e., researchers have no idea about SNP strand (ambiguous SNPs usually are C/G or A/T SNPs), and (iv) missing alleles.

Chip heritability. The chip heritability ($h_{\mathit{SNP}}^2$)is also known as SNP-based heritability, which is defined as the portion of the phenotypic variation that the genotyped genetic marker can explain.^60,61 Higher values of heritability indicate that the phenotype is explained best by the genotype, i.e. set of SNPs. Choi et al.⁵⁸ recommend $h_{\mathit{SNP}}^2>0.05$ to perform PRS analysis. To estimate $h_{\mathit{SNP}}^2$, researchers should use LD score regression that could be used to distinguish polygenicity (SNPs effects) and confounding biases, including cryptic relatedness and population stratification.⁶²

The code below demonstrates how to use Plink to perform quality control checks and calculating PRS.

# $bfile: Target data set in Plink binary format.
# $QC: Base data set i.e., base GWAS summary statistic. This file contains P-value information
# $p1: P-value threshold for a SNP to be included as an index SNP. Choose value of 1 if you want to
##include all SNPs are include for clumping.
# $r2: Cutoff for r2 value i.e., SNPs having value higher than given r2 will be removed
# $kb: cutoff value window size in kilobase, i.e., SNPs within $kb of the index SNP are considered for clumping.
# $SNP: Name the column SNP that containing the SNP IDs
# $P: Name of the column that containing the P-value information
# $Output: Output file

################### Quality Control of Target Samples
plink\
  --bfile $bfile \
  --maf 0.05 \
  --mind 0.1 \
  --geno 0.1 \
  --hwe 1e-6 \
  --make-just-bim \
  --make-just-fam \
  --out $Output.qc

################### Clumping
plink \
  --bfile $bfile \
  --clump-p1 $p1 \
  --clump-r2 $r2 \
  --clump-kb $kb \
  --clump $QC \
  --clump-snp-field $SNP \
  --clump-field $P \
  --out $Output

############### calculating PRS
# $listpvalue: Threshold values
# $valid.snp: valid SNPs
# Here we are assuming that 3 column for SNP ID; 4 for effective allele information;
### the 12 for effect size estimate.
#### Moreover, this file has a header.

plink \
  --bfile $bfile \
  --score $QC 3 4 12 header \
  --q-score-range $listpvalue \
  --extract $valid.snp \
  --out $Output

Meta-analysis

The meta-analysis approach can be used to evaluate the association between SNPs and biological traits by combining various GWAS summaries from different studies. The following paragraphs will provide the key concepts for performing a meta-analysis. To have concrete information about the meta-analysis approach, we recommend this review article.⁶³

Heterogeneity of source. Heterogeneity in data could be derived using GWAS summary statistics. The standard variables to estimate heterogeneity in data include odds ratios, standardized effect sizes, other metrics along with their uncertainty (e.g. variance or 95% confidence interval) and the accompanying p-values. However, there might be many other variables for each dataset that are important to deal with in order to estimate heterogeneity in data.

Standard meta-analysis. This is used to perform the meta-analysis approach which is to sum the Z-scores across all studies and weigh them appropriately using the sample sizes. See equation (3) below.

where $Z_k$ is $Z$-score from $K^{t}h$ study, and $w_k$ weight of studies relative to population size.

Independence of the samples. Conventional meta-analysis has an assumption that assumes that effect sizes are independent. Simulation studies demonstrate that failure to account for overlapping samples could greatly inflate type I error.⁶⁴ If accounting for the overlap is unavoidable, the overlap/covariance can be estimated using Z’s covariance between summary statistics. Some tools can account for overlapping samples, such as METAL software⁶⁵ and ASSET.⁶⁶

Correcting for population structure with genomics control. The presence of population structure in the GWAS study can impact an over-dispersion of the corresponding association test statistics. One approach to limit this problem is to correct statistics of each summary using genomic control. This correction factor is given as the inflation factor ($\lambda$) which is the test statistics’ median divided by its expectation under the null hypothesis.⁶⁷

p-values versus Z scores. Meta-analysis methods based on p-values were widely used in different scientific fields until the $1980s$. Then became unpopular and almost abandoned in biomedical sciences. Nowadays, the meta-analysis approach is performed using Z-score. There are two methods to estimate Z-score for GWAS data. The first method is demonstrated in equation (4).

In the second method (equation (5)), Z-score is estimated using p-value and the effect of allele.

where $\text{sign}\left({\mathrm{\Delta }}_i\right)$ is a sign of relation.

Random effects versus fixed effects. In the presence of variability of allelic effect as in trans-ethnic studies, it is common to perform a random-effect meta-analysis to correct variability of $\beta$ effect between different studies. For instance, GWAMA tool⁶⁸ computes a random-effects variance component using Cochran’s statistic (Q-value) to balance weight used in meta-analysis. On the other hand, Metasoft⁶⁹ proposed two other different methods to take into account the heterogeneity, which are Random Effects model⁷⁰ and binary effects model with m-value, i.e. weight of each summary study in summary statistics.⁷¹

Heterogeneity test. Heterogeneity at a locus can be reflected in the variability in population or environment. It can be relevant to gene-environment interaction and the reason behind the variability in GWAS approach of each data-set, i.e. covariable, model and approximation of GWAS between summary statistics. The heterogeneity is computed between two sets of summary statistics rather than one locus.

Cochran’s statistic provides a test of heterogeneity of allelic effects at SNPs $j$ using equation (6) below.

where $N$ denotes study number.

Alternatively, we use Q statistic to quantify the extent of heterogeneity in allelic effects across studies.⁷² See equation (7) below

Other meta-analysis approaches. Traditionally, meta-analyses of GWAS have focused on combining results of multiple studies for similar traits. The Bayesian framework has been tested to estimate $\beta$ on different phenotypes.⁷³ MetABF tool⁷⁴ has implemented a method to perform meta-analysis across genome-wide association studies of diverse phenotypes. It is important to note that a recent review on cancer suggests that it is possible to obtain “noteworthy” Bayesian results at higher p-values that are not considered statistically significant in GWAS.⁷⁵

Study alignment and error trapping. Meta-analysis aggregates various summary statistics. Therefore, any error to designate the effect allele and other allele or strand issue can cause an error in estimate $\beta$. Such errors might lead to misleading meta-analysis results as it increases $Q$ of Cohran and heterogeneity between studies.

Effect size. By conducting a meta-analysis, researchers often neglect the sample size variation among different studies “true effect sizes are the same across studies”. However, in some cases, researchers introduce the correct effect size by considering the posterior probability for each study. Some software for estimating effect size are given in Table 3.

Meta-analysis output. Meta-analysis provides a new set of summary statistics. For each position that has not been discarded, new statistics will be calculated. These statistics include new values for $\beta$, $\sigma$ and p-value. Users should be aware of the LD and the reference population used.

More resources for meta-analysis. We recommend h3bionet/h3agwas for meta-analysis pipeline. In addition, those interested in meta-analysis are directed to a review published by Zeggini et al.⁷⁷ and another review article by Evangelou and Ioannidis.⁶³ The code below demonstrates how to use the Metasoft tool to perform meta analysis.

# example with metasoft
## reformatting input files.
### File 1
R -e "library(data.table);data.gwas<-fread('$File 2');
data.gwas<-data.gwas[,c('CHR','SNP','BP','ALLELE1',
'ALLELE0','BETA','SE','P_BOLT_LMM')];
names(data.gwas)<-c('CHR', 'SNP','BP','A1',
'A2','BETA', 'SE','P');
data.gwas<-data.gwas[data.gwas[['P']]<1];
write.table(data.gwas, file='gwas_file2.qassoc',
row.names=F, col.names=T, sep='\t', quote=F)"
### File 2
R -e "library(data.table);data.gwas<-fread('$File 1');
data.gwas<-data.gwas[,c('CHR','SNP','BP','ALLELE1',
'ALLELE0','BETA','SE','P_BOLT_LMM')];
names(data.gwas)<-c('CHR', 'SNP','BP','A1',
'A2','BETA', 'SE','P');
data.gwas<-data.gwas[data.gwas[['P']]<1];
write.table(data.gwas, file='gwas_file1.qassoc',
row.names=F, col.names=T, sep='\t', quote=F)"
## metasoft
### Metasoft with binary effect
java -jar Metasof/Metasoft.jar \
-input file_merge_all.meta\
-output meta.meta \
-pvalue_table Metasof/HanEskinPvalueTable.txt\
-binary_effects\

Colocalization analysis

Colocalization. Colocalization is an approach used to integrate annotations with GWAS results. The annotations resources include gene expression (eQTLs), protein expression (pQTLs), exon splicing (sQTLs), DNA methylation (mQTLs), and chromatin acetylation and chromatin accessibility (caQTLs).⁷⁸

Statistics for colocalization studies. Several parametric and non-parametric statistics can be done for colocalization studies.⁷⁹

Resources for colocalization studies. We recommend the following scientific resources for those who are beginners in this field:

Using summary statistics from multiple phenotypes and traits. The methods for GWAS are mostly focused on single variant analysis with a single phenotype or trait. Increasing evidence shows that pleiotropy, one gene’s effect on multiple phenotypes, plays a pivotal role in many complex traits. Therefore, associating different GWAS results for multiple phenotypes can provide an extensive power by aggregating multiple weak signals.⁸³ Different approaches have been developed to integrate dependent p-values to assess the association between a gene and multiple correlated phenotypes.⁸³ Several tools exist to perform multi traits analysis, including Multi-Trait Analysis of GWAS (MTAG)⁸⁴ and CPASSOC package.⁸⁵

Mendelian randomisation

Mendelian randomisation (MR) is a statistical approach that can be defined as “the use of genetic variants as instrumental variables to investigate the effects of modifiable risk factors for disease”.⁸⁶ For instance, one trait (phenotype or disease) might be affected by confounding or reverse causation rather than a conventional observational variable. Therefore, MR aims to provide a statistical frame to verify the causality between locus and phenotype and exclude pleiotropy. Such methods will provide a reliable explanation of the results.

Assumption of MR. There are three main assumptions for MR.^87-90 These are (i) the genetic variant that is associated with the exposure (significant association), (ii) the genetic variant that is independent of the outcome given to the exposure and all confounders (measured and unmeasured) of the exposure-outcome association, (iii) the genetic variant that is independent of factors (measured and unmeasured) that confound the exposure-outcome relationship.

Statistical methods for MR. MR general strategy is to compare $\mathit{\text{beta}}$ values or p-value of the same position of two or more different GWAS using related/confounding phenotype. Several methodologies have been developed for MR analysis, an example is the ratio of coefficients estimator, which can be modeled using equation (8):

where ${\beta }_e$ represents the change in exposure per variant allele, and ${\beta }_o$ represents the change in outcome per variant allele. Another model is the two-stage least squares. It employs a two-stage regression approach with two regression models where the first stage regression’s output is used as the input of the second stage regression. More methods for MR exist, including control function estimator, limited information maximum likelihood method, verse variance weighted method, and MR-Egger method.⁹¹

MR software. Several tools have been developed to compute MR, e.g GSMR (Generalised Summary-data-based Mendelian Randomisation) from GCTA, TwoSampleMR (version 0.4.20), MR-PRESSO, and MR-LDP which is integrated LD information. We recommend the following tutorial: https://bioconductor.org/packages/release/bioc/vignettes/GMRP/inst/doc/GMRP.pdf.⁹²

The code below demonstrates how to use gcta64 tool to perform Mendelian randomisation analysis.

#reformat file 1
R -e "library(data.table);
data.gwas<-fread('$File');
data.gwas<-data.gwas[,c('SNP','ALLELE1','ALLELE0','A1FREQ','BETA','SE','P_BOLT_LMM')];
names(data.gwas)<-c('SNP','A1', 'A2','freq','b', 'se','p');
data.gwas[['N']]<-10000;
data.gwas<-data.gwas [data.gwas[['p']]<1];
SNPtmp<-table(data.gwas[['SNP']]);
UniqSNP<-names(SNPtmp)[SNPtmp==1];
write.table(data.gwas[data.gwas[['SNP']] %in% UniqSNP],file='FilePheno1',row.names=F,
col.names=T, sep='\t', quote=F)"

#reformat file 2
R -e "library(data.table);
data.gwas<-fread('$File');
data.gwas<-data.gwas[,c('SNP','ALLELE1','ALLELE0','A1FREQ','BETA','SE','P_BOLT_LMM')];
names(data.gwas)<-c('SNP','A1', 'A2','freq','b', 'se','p');
data.gwas[['N']]<-10000;
data.gwas<-data.gwas[data.gwas[['p']]<1];
SNPtmp<-table(data.gwas[['SNP']]);
UniqSNP<-names(SNPtmp)[SNPtmp==1];
write.table(data.gwas [data.gwas[['SNP']] \%in% UniqSNP],
file='FilePheno2',row.names=F, col.names=T,
sep='\t', quote=F)
# create
echo "P1 FilePheno1" > exposure
echo "Exp FilePheno2" > outcome

## gcta need a plink file or bgen file
## GSMR analyses,
###forward-GSMR analysis (coded as 0),
###reverse-GSMR analysis (coded as 1)
###and bi-GSMR analysis (both forward- and reverse-GSMR analyses, coded as 2).
gcta64 --bfile plkfile --gsmr-file exposure outcome \
--gsmr-direction 0 --out test_gsmr_result

MR and gene expression. Recently, some methods have integrated MR into GWAS and eQTL to test if the effect of gene expression is zero on the trait.^93-95 These methods provide a promising way to combine GWAS summary statistics and expression data.