EPACTS for DIAGRAM

Motivation and Rationale

EPACTS is a software pipeline developed to perform various statistical tests for analysis of whole-genome / whole-exome sequencing data.  The main motivation for using EPACTS is to use a consistent analysis framework for association analysis in the DIAGRAM consortium.  In addition, for analysis of low frequency variants (minor allele frequency [MAF] < 5%), standard logistic regression Wald or likelihood ratio tests found in existing association software are conservative or anti-conservative respectively.  We implemented two statistical tests we will use for for analysis of low frequency variants: (1) logistic regresion-based score test and (2) Firth bias-corrected logistic regression (Firth, 1993).  For analysis of common variants, any asyptotic logistic regression test has well-controlled type I error rates and asymptotically equivalent power.  For simplicity and consistency, we propose the use of both score and Firth tests for testing all allele frequencies.

Outline of analysis protocol

This is an overview of the analysis protocol for analyzing imputed DIAGRAM datasets using the EPACTS pipeline.  We assume that your dataset has been imputed using minimac or Impute2.  Starting with minimac or impute2 output:

1. Download and install EPACTS
2. Prepare VCF file with genotypes / dosages
3. Prepare PED file with phenotypes and covariates
4. Run EPACTS association pipeline
5. Report association results in appropriate format

1.  Download and install EPACTS

For External Users

For external users, follow the instruction at EPACTS page, summarized below.

• Download EPACTS binary at http://www.sph.umich.edu/csg/kang/epacts/download/epacts.v2.2.0.20121026.tar.gz (94MB)
• Uncompress EPACTS package to the directory you would like to install

 tar xzvf epacts.v2.2.0.20121026.tar.gz

• Download the reference FASTA files from 1000 Genomes FTP automatically by running the following commands
  

cd epacts.v2.2.0.20121026/
./ref_download.sh
(For advanced users, to save time for downloading the FASTA files (~900MB), you may copy a local copy of GRCh37 FASTA file and the index file to ${EPACTS_DIR}/ext/ref/)

• Perform a test run by running the following command

epacts.v2.2.0.20121026/example/test_run_epacts.sh

For Local Users in CSG

For users in CSG, EPACTS can be found here:

/net/fantasia/home/hmkang/sw/epacts2.1/

Getting started with an example

Once installed, test out the software by running a quick example using the test data provided in the "example" directory. The example VCF and PED files are:

$ epacts.v2.2.0.20121026/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz

$ epacts.v2.2.0.20121026/example/1000G_dummy_pheno.ped

Run the single variant score test on the example data using this command:

$ epacts.v2.2.0.20121026/epacts single 
--vcf epacts2.1/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz 
--ped epacts2.1/example/1000G_dummy_pheno.ped 
--min-maf 0.001 --chr 20 --pheno DISEASE --cov AGE --cov SEX --test b.score --anno 
--out {OUTPUT_DIR}/test --run 2 &

This command will run the single variant test on the input VCF and PED files, with a minimum MAF threshold of 0.001.  The phenotype is "DISEASE" and we are adjusting the analysis with covariates AGE and SEX.  The output file directory prefix is {OUTPUT_DIR}/test.  Finally, EPACTS will run the analysis in parallel on 2 CPUs.

A more detailed description of the example can be found here.

2.  Prepare VCF file with genotypes / dosages

EPACTS requires input genotype / dosage information in VCF format.  From minimac or Impute2, you wil start with your imputed dosage file.

A.  Convert dosage file into VCF format

Use the wrapper program "dose2vcf" to convert your doseage output to VCF format.  Download the tool from here. If you used rs numbers during imputation, you can find mapping tables ready for dose2vcf here (214 Mb)

To run the wrapper program, use the following command

> dose2vcf/dose2vcf --dose [ .dose file ]  --info [ .info file ] --out [ output file prefix ]

For example, for imputed chromosome 1 data from FUSION:

> dose2vcf/dose2vcf --dose FUSION.GWAS.1KG.imp.chr1.dose --info FUSION.GWAS.1KG.imp.chr1.info --out out/FUSION.GWAS.1KG.imp.chr1

The expected output file is:

out/FUSION.GWAS.1KG.imp.chr1.vcf

Note that for longer chromosomes, the program is quite memory intensive.  In this case, please convert dosages in shorter sections of the chromosome.  For example, if the imputation was performed by sections, then convert these sections to vcf first, and then merge the vcf files together using vcftools vcf-concat:

B.  bgzip and tabix VCF files

Input VCF file must be bgzipped and tabixed before running association to allow efficient random access of the file. Below is an example command to conver plain VCF into bgzipped and tabixed VCF

bgzip input.vcf ## this command will produce input.vcf.gz
tabix -pvcf -f input.vcf.gz ## this command will produce input.vcf.gz.tbi

If the VCF file is separated by chromosome, the VCF file must contain the string "chr1" in the chromosome 1 file, and corresponding chromosome name for other chromosomes.
Sample IDs in the VCF file must be consistent to those from PED file

3.  Prepare PED file for phenotypes and covariates

EPACTS accepts the PED format supported by MERLIN or PLINK to represent the phenotypes and covariates.  You may prepare either (1) a PED file without column headers + accompanying DAT file, or (2) a PED file with column headers.  The standard PED format has 6 mandatory columns:

1. Family ID
2. Individual ID
3. Paternal ID
4. Maternal ID
5. Sex (1=male; 2=female; other=unknown)
6. Phenotype (1 = control; 2 = case)

Columns 7 and onwards are additonal covariates and or phenotypes.  For example

1. QT
2. AGE
   

etc.

Note:  Categorical covariates must be coded as dummy variables!  For example, Sex cannot be coded as "M" or "F".

An example PED file with a header is as follows.  Note that the header must start with a "#" symbol.

#FAM_ID IND_ID FAT_ID MOT_ID SEX DISEASE QT AGE
13281 NA12344 NA12347 NA12348 1 1 94.17 66.1
13281 NA12347 0 0 1 1 109.54 44.0
13281 NA12348 0 0 2 2 119.40 46.6
1328 NA06984 0 0 1 2 87.72 39.3
1328 NA06989 0 0 2 1 100.60 41.7
1328 NA12329 NA06984 NA06989 2 1 100.85 46.4
13291 NA06986 0 0 1 2 91.94 61.9
13291 NA06995 NA07435 NA07037 1 2 104.36 57.4
13291 NA06997 NA06986 NA07045 2 2 107.53 53.1

Alternatively, you can prepare a PED file without a header, and include a corresponding DAT file describing the column headers.

13281 NA12344 NA12347 NA12348 1 1 94.17 66.1
13281 NA12347 0 0 1 1 109.54 44.0
13281 NA12348 0 0 2 2 119.40 46.6
1328 NA06984 0 0 1 2 87.72 39.3
1328 NA06989 0 0 2 1 100.60 41.7
1328 NA12329 NA06984 NA06989 2 1 100.85 46.4
13291 NA06986 0 0 1 2 91.94 61.9
13291 NA06995 NA07435 NA07037 1 2 104.36 57.4
13291 NA06997 NA06986 NA07045 2 2 107.53 53.1

The corresponding DAT file is:

A DISEASE
T QT
C AGE

Key:  A = binary trait; T = quantitative trait; C = covariate

4.  Run EPACTS association pipeline

For detailed description of options, use:

epacts.v2.2.0.20121026/epacts single -man

Primary analyses

There are 4 separate association analyses to be completed.

Association Analysis	Statistical Test	Subset of SNPs	Output File Type	Output Filename Format
1.  Typical DIAGRAM analysis using existing association pipeline (with or without BMI)	Wald or likelihood ratio	All SNPs with MAC >= 1	Custom file based on DIAGRAM format	A. DIAGRAMv4_iSNPs_XXX_1000G_KKK_TTT_YYY_ZZZ.txt B. DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_TTT_YYY_ZZZ.txt
2.  Analysis of low frequency variants using Firth bias-corrected logistic regression (with or without BMI)	Firth bias-corrected	SNPs with 200 >= MAC >= 1	EPACTS output file	A. DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz B. DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

Analysis for QC

For quality control, please run an additional analysis using EPACTS on all SNPs for chromosome 20 only using the SCORE test without BMI adjustment.  These results will be used to compare with results from the primary analyses, to ensure the new EPACTS software has been run correctly.

Association Analysis	Statistical Test	Subset of SNPs	Output File Type	Output Filename Format
3.  Analysis of chromosome 20 using logistic regression score test (without BMI)	Score test	Chromosome 20 SNPs with MAC >= 1	EPACTS output file	A. DIAGRAMv4_iSNPs_XXX_1000G_KKK_SCR_YYY_ZZZ.epacts.gz

Filename key:

iSNPs indicates the file contains imputed SNPs only
XXX indicates a uniquely identifiable STUDY NAME: (e.g.WTCCC, DGI, DGDG, FUSION, ERGO, DUNDEE, NHS,FHS, TYROL, EUROSPAN etc.)
“_adjBMI” – should be present for the model with BMI adjustment
KKK indicates date of the 1000Genomes map, Month and year (defined as MMMYY) of the 1000Genomes map that was used e.g. MAR12
TTT indicates the test used to evaluate significance

• SCR = Score
• WLD = Wald
• LHR = Likelihood ratio
• FBC = Firth Bias Corrected Logistic Regression

YYY indicates the DATE of file generation (MMDDYY format, e.g. 021710 – apologies in advance to our European colleagues)
ZZZ indicates the name + other initials of the uploader (e.g., BFV, LJS, ABC, etc.)

Comparative runtimes for Firth and score tests.

Firth test is relatively computationally intensive compared to the score test.  However, Firth has better test properties (i.e. type I error rates and power) compare to score test when analyzing low frequency variants in studies with unequal numbers of cases and controls.

For example, using a single CPU, analyzing 2333 FUSION individuals, chromosome 20, 22755 SNPs (adjusting for age, sex and 11 dummy variables for birth province):

• Firth test:  56 minutes (3360 seconds)
• Score test:  32 seconds

Hence, we only ask for Firth test results for SNPs with MAC <= 200.

1. Typical DIAGRAM analysis using existing association pipeline

This is the typical DIAGRAM analysis using your current association pipeline and software.   File:1000Genomes march2012 imputation analysis plan 08312012.pdf

Alternative:  Analyze VCF and PED files using the Wald test with the EPACTS software:

This will be a good way to check if your VCF and PED files are correctly formatted for EPACTS.

epacts.v2.2.0.20121026 /epacts single -vcf [INPUT VCF FILENAME] -ped [INPUT PED FILENAME] -out [OUTPUT FILENAME PREFIX] \
-test b.wald -pheno DISEASE -cov AGE -sepchr -anno -min-mac 1 -field EC -run 10

Important: To analyze dosages (not genotypes), you must specify the dosage field with the "--field EC" option. Without this option, you will be analyzing the hard genotypes (i.e. --field option defaults to "GT" or "genotypes")!

2. Analysis of low frequency variants using Firth bias-corrected logistic regression

The Firth bias-corrected test has well-controlled type I error rate and good power for analysis of balanced and unbalanced studies.  However, it is more computationally intensive.  We only run Firth on the subset of variants with 1<= MAC <= 200.

To run the Firth test using the EPACTS software:

epacts.v2.2.0.20121026 /epacts single -vcf [INPUT VCF FILENAME] -ped [INPUT PED FILENAME] -out [OUTPUT FILENAME PREFIX] \
-test b.firth -pheno DISEASE -cov AGE -sepchr -anno -min-mac 1 -max-mac 200  -field EC -run 10

Important:  To analyze dosages (not genotypes), you must specify the dosage field with the "--field EC" option.  Without this option, you will be analyzing the hard genotypes (i.e. --field option defaults to "GT" or "genotypes")!

3. Analysis of chromosome 20 using logistic regression score test

The score test has well-controlled type I error rate and good power for meta-analysis of balanced (equal numbers of cases and controls) studies.  It is also very computationally efficient.  Please run the score test using the EPACTS software.

The EPACTS command for the score test analysis of chromosome 20 is:

epacts.v2.2.0.20121026 /epacts single -vcf [INPUT VCF FILENAME] -ped [INPUT PED FILENAME] -out [OUTPUT FILENAME PREFIX] \
-test b.score -pheno DISEASE -cov AGE -chr 20 -anno -min-mac 1 -field EC -run 10

This command will run single variant analysis using the score test logistic regression on the DISEASE phenotype adjusting for AGE. Add the relevant additional covariates with additional "-cov" options. This assumes that the VCF files are separated by chromosomes (option -sepchr). All variants with at least one minor allele count will be analyzed (option -min-mac 1). It will annotate results by functional category (option -anno) and run the analysis on 10 parallel CPUs (option -run 10).

5.  Report EPACTS results

For analyses 2 and 3, please upload the three .epacts.gz files to the FTP server:

1. Firth test (no BMI):  DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz
2. Firth test (with BMI):  DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_FBC_YYY_ZZZ.epacts.gz
3. Score test (chr 20, no BMI):  DIAGRAMv4_iSNPs_XXX_1000G_KKK_SCR_YYY_ZZZ.epacts.gz

Here's an example score test .epacts file

$ head test.epacts
#CHROM  BEGIN   END     MARKER_ID       NS      AC      CALLRATE        MAF     PVALUE  SCORE   N.CASE  N.CTRL  AF.CASE AF.CTRL
20      68303   68303   20:68303_A/G_Upstream:DEFB125   266     1       1       0.0018797       NA      NA      NA      NA      NA      NA
20      68319   68319   20:68319_C/A_Upstream:DEFB125   266     0       1       0       NA      NA      NA      NA      NA      NA
20      68396   68396   20:68396_C/T_Nonsynonymous:DEFB125      266     1       1       0.0018797       NA      NA      NA      NA      NA      NA
20      76635   76635   20:76635_A/T_Intron:DEFB125     266     0       1       0       NA      NA      NA      NA      NA      NA
20      76689   76689   20:76689_T/C_Synonymous:DEFB125 266     0       1       0       NA      NA      NA      NA      NA      NA
20      76690   76690   20:76690_T/C_Nonsynonymous:DEFB125      266     1       1       0.0018797       NA      NA      NA      NA      NA      NA
20      76700   76700   20:76700_G/A_Nonsynonymous:DEFB125      266     0       1       0       NA      NA      NA      NA      NA      NA
20      76726   76726   20:76726_C/G_Nonsynonymous:DEFB125      266     0       1       0       NA      NA      NA      NA      NA      NA
20      76771   76771   20:76771_C/T_Nonsynonymous:DEFB125      266     3       1       0.0056391       0.68484 0.40587 145     121     0.013793        0.0082645

Each column represents

• CHROM : Chromosome Name
• BEGIN, END : Base position of the variant on each side
• MARKER_ID : [CHROM]:[POS]_[REF]/[ALT]_[ANNOTATION] formatted marker ID. [ANNOTATION] information will be available only with --anno option
• NS : Number of samples with non-missing genotypes
• AC : Non-reference allele count
• CALLRATE : Genotype call rate
• MAF : Minor allele frequency
• PVALUE : P-value

The rest of columns varies by statistical tests. For example, in b.score test, SCORE represents score test statistics, N.CASE and N.CTRL represents the case/control counts, and AF.CASE and AF.CTRL represents the case/control allele frequencies.

Troubleshooting Common Issues

Errors when running EPACTS

Rscript execution error: No such file or directory

epacts2.1/epacts single \

--vcf epacts2.1/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz \
--ped epacts2.1/example/1000G_dummy_pheno.ped \
--min-maf 0.001 --chr 20 --pheno DISEASE --cov AGE --cov SEX --test b.score --anno \
--out out/test --run 1 &
[4] 13569

$ Detected phenotypes with 2 unique values - 1 and 2 - considering them as binary phenotypes... re-encoding them into 1 and 2
Successfully written phenotypes and 2 covariates across 266 individuals
Processing chromosome 20...
Finished generating EPACTS Makefile
Running 1 parallel jobs of EPACTS
epacts2.1/bin/..//bin/make -f out/test.Makefile -j 1

Rscript execution error: No such file or directory

The issue seems to be the absence with Rscript, which is usually installed with R.

If you can find Rscript (e.g. /usr/bin/Rscript, /usr/local/bin/Rscript), or if you can re-install the full Rscript, you can simply avoid the problem by setting your environment variable.

Otherwise, Hyun will modify EPACTS to not requiring this (so you can run R CMD BATCH instead of Rscript).

ERROR: No overlapping IDs between VCF and PED file. Cannot proceed.

Check that your individual ID's in your PED file are the same as those in your VCF file.

For example, if your VCF individual ID's include the family ID's (i.e. ABCD->ABCD001), the individual ID's in the PED file must match it exactly.

Estimated allele frequencies and analysis results do not exactly match results from my existing association software

Check that you have included the same set of covariates (with categorical variables encoded as dummy variables).

Check that you have the same number of cases and controls analyzed.

FInally, check that you used dosages by adding the appropriate "-field" option.  For example, suppose your VCF is:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT A001 B001 C001
11 180567 11:180567 C G 0 PASS . GT:EC 1/1:2.0000 1/1:2.0000 1/1:2.0000
11 186458 11:186458 G A 0 PASS . GT:EC 1/1:1.9850 1/1:1.9750 1/1:1.9840
11 186462 11:186462 C A 0 PASS . GT:EC 1/1:2.0000 1/1:2.0000 1/1:2.0000
11 192958 11:192958 G T 0 PASS . GT:EC 1/1:2.0000 1/1:2.0000 1/1:2.0000
11 192995 11:192995 C T 0 PASS . GT:EC 1/1:2.0000 1/1:2.0000 1/1:2.0000
11 193065 11:193065 G A 0 PASS . GT:EC 1/1:1.9980 1/1:1.9990 1/1:1.9960
11 193096 11:193096 C T 0 PASS . GT:EC 0/1:0.7840 0/1:0.6280 1/1:1.6550
11 193146 11:193146 G A 0 PASS . GT:EC 1/1:1.8550 1/1:1.8460 1/1:1.7940

The genotype information has FORMAT "GT:EC".  For the first SNP (chr11:180567) and individual A001, the genotype is 1/1 and dosage is 2.0000.  To access the dosages, you must specify the option "-field EC"