Genome Analysis Wiki - User contributions [en]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:50:22Z

Clement Ma: /* Single Variant Association Analysis */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --min-mac 1 --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE BETA SEBETA ZSTAT
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.28783 0.80648 0.75876 1.0629
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.99286 -17.704 1978.1 -0.0089502
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.643 -0.44056 0.95048 -0.46351
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.096774 12 4 0 0.3989 0.76461 0.90638 0.84358
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.076266 1.9741 1.1136 1.7728
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.048387 3 1 0 0.9904 16.668 1385.4 0.012031
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 5.9095e-16 1.4376 4.1107e-16
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:49:12Z

Clement Ma: /* Single Variant Association Analysis */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score-min-mac 1 --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE BETA SEBETA ZSTAT
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.28783 0.80648 0.75876 1.0629
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.99286 -17.704 1978.1 -0.0089502
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.643 -0.44056 0.95048 -0.46351
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.096774 12 4 0 0.3989 0.76461 0.90638 0.84358
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.076266 1.9741 1.1136 1.7728
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.048387 3 1 0 0.9904 16.668 1385.4 0.012031
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 5.9095e-16 1.4376 4.1107e-16
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:23:12Z

Clement Ma: /* Run Groupwise Test */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE BETA SEBETA ZSTAT
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.28783 0.80648 0.75876 1.0629
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.99286 -17.704 1978.1 -0.0089502
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.643 -0.44056 0.95048 -0.46351
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.096774 12 4 0 0.3989 0.76461 0.90638 0.84358
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.076266 1.9741 1.1136 1.7728
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.048387 3 1 0 0.9904 16.668 1385.4 0.012031
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 5.9095e-16 1.4376 4.1107e-16
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:22:21Z

Clement Ma: Undo revision 12372 by Clement Ma (talk)

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:21:31Z

Clement Ma: /* Run Groupwise Test */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE BETA SEBETA ZSTAT
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.28783 0.80648 0.75876 1.0629
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.99286 -17.704 1978.1 -0.0089502
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.643 -0.44056 0.95048 -0.46351
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.096774 12 4 0 0.3989 0.76461 0.90638 0.84358
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.076266 1.9741 1.1136 1.7728
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.048387 3 1 0 0.9904 16.668 1385.4 0.012031
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 5.9095e-16 1.4376 4.1107e-16
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-12T12:16:22Z

Clement Ma: /* Single Variant Association Analysis */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
Interpretation for top associated variant chr22:36995620_A/G:
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516

The score test PVALUE = 5.6717e-09 is strongly significant.
For the notation chr22:36995620_A/G, the reference allele is the first allele "A", and the non-reference (effect) allele is the second allele "G".
The SCORE statistic = 0.29032 is positive (e.g. SCORE>0), meaning that the effect allele INCREASES your risk for disease. Note that the SCORE is NOT the BETA or log odds ratio. To calculate the BETA and SE(BETA), you must run the b.wald test.
There is a large difference between the allele frequencies for cases and controls. AF.CASE = 0.51613 >> AF.CTRL = 0.0064516.
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

== Return to Ancestry on your own genome==
Let's go see if Ancestry finished: [[SeqShop:_Ancestry_On_Your_Own_Genome,_December_2014#Checking_if_Pileup_finished]]

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-11T17:04:27Z

Clement Ma: /* Single Variant Association Analysis */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

head $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

SeqShop: Genetic Association Analysis Practical, December 2014

2014-12-11T16:47:18Z

Clement Ma: /* Setup your run environment */

== Introduction ==
[[Media:Seqshop association 2014 06.pdf|View Lecture Slides]]

[[Media:Seqshop association practice 2014 06.pdf|View Introductory Slides for Practical Session]]

== Goals of This Session ==

* Understand how to annotate variants using EPACTS
* Understand how to run single variant association analysis using EPACTS
* Understand how to run rare variant association test using EPACTS
* Understand how to visualize the association output from EPACTS

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not already logged in, please expand this section.''
<div class="mw-collapsible-content">
{{SeqShopLogin}}
</div>
</div>

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* [[GotCloud]] program
* Tutorial input files
* Setup an output directory
** It will leave your output directory from the previous tutorial in tact.
source /net/seqshop-server/home/mktrost/seqshop/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment
** If you want to view the detail of the setup, type
less /net/seqshop-server/home/mktrost/seqshop/setup.txt
and press 'q' to finish.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View setup.txt
<div class="mw-collapsible-content" style="width:800px">
export GC=/net/seqshop-server/home/mktrost/seqshop/gotcloud
export SS=/net/seqshop-server/home/mktrost/seqshop/example
export EPACTS=/net/seqshop-server/home/mktrost/seqshop/epacts
export OUT=~/out
mkdir -p ${OUT}
</div>
</div>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial builds on the alignment & snpcall tutorials, if you have not already, please first run those tutorials: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical|Alignment Tutorial]] & [[SeqShop: Variant Calling and Filtering for SNPs Practical|SNP Calling Tutorial]]

=== Download & Build EPACTS ===
If you do not already have EPACTS:
* cd to where you want EPACTS installed (you can change this to any directory you want)
mkdir -p ~/seqshop
cd ~/seqshop/
* download, decompress, and build the version of epacts that was tested with this tutorial:
wget http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.6.tar.gz
tar xvf EPACTS-3.2.6.tar.gz
cd EPACTS-3.2.6
./configure --prefix=$HOME/seqshop/epacts
make
make install
cd ../..

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for EPACTS:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export EPACTS=~/seqshop/epacts</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv EPACTS ~/seqshop/epacts</pre>
</div>
</div>
</ul>
</ul>

</div>
</div>

== Preparing Input Files ==

=== Input VCF file ===

We will use SNP genotypes from the SNP calling session, after LD-aware genotype refinement.
Check the contents of the VCF file using the following command.

zless ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz

=== Phenotype Information ===

Phenotype information is prepared in PED format commonly used in other GWAS software such as MERLIN or PLINK.

less ${SS}/assoc/seqshop.ped

The first several line should look like below.

<div class="mw-collapsible mw-collapsed" style="width:200px">
View Data
<div class="mw-collapsible-content" style="width:800px">
#FAM_ID IND_ID DAD_ID MOM_ID SEX PHENO
HG00551 HG00551 0 0 0 0
HG00553 HG00553 0 0 0 0
HG00554 HG00554 0 0 0 0
HG00637 HG00637 0 0 0 0
HG00638 HG00638 0 0 0 0
HG00640 HG00640 0 0 0 1
HG00641 HG00641 0 0 0 1
HG00734 HG00734 0 0 0 1
HG00736 HG00736 0 0 0 0
...
</div>
</div>

Binary phenotype can be encoded as 0-1 or 1-2. If the column contains more than two distinct values, it will automatically be recognized as quantitative values.

EPACTS allows PED file to have a header line. The header line should contain the description of each column. EPACTS also accepts a standard PED format where .ped file contains the phenotype data and .dat file contains the information about each column.

=== Installed version of EPACTS ===

EPACTS are installed in the server. If you want to install EPACTS by yourself, visit [[EPACTS]] page for more details

ls $EPACTS/bin

<div class="mw-collapsible mw-collapsed" style="width:400px">
View EPACTS executable files
<div class="mw-collapsible-content" style="width:800px">
anno epacts epacts-cis-extract epacts-group epacts-multi epacts.pm epstopdf test_run_epacts.sh
bgzip epacts-anno epacts-download epacts-make-group epacts-pca-plot epacts-single pEmmax vcfast
chaps epacts-cat epacts-enrich epacts-make-kin epacts-plot epacts-zoom tabix wGetOptions.pm
</div>
</div>

Note that some tools undocumented in [[EPACTS]] documentation is under development and may not work.

== Annotating Variants with EPACTS ==

There are multiple software tools that provides a function to annotate variants, such as Variant Effect Predictor (VEP) that is used in 1000 Genomes Project. While most annotation software provides very similar results to each other, their computational efficiency can substantially vary. The annotation software EPACTS provides is extremely fast and can provide genome-wide annotation results in orders of magnitude faster than other widely available annotation software.

In order to annotate variants with EPACTS, one can use <code>epacts-anno</code> module.

mkdir --p $OUT/assoc
$EPACTS/bin/epacts-anno --in $OUT/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --out $OUT/assoc/snps.anno.vcf.gz --ref $SS/ref22/human.g1k.v37.chr22.fa

Then you will see a series of messages before annotation finishes.

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the expected messages
<div class="mw-collapsible-content" style="width:800px">
/home/hmkang/seqshop/epacts/bin/anno -i /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz -r \
/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa -f refGene -g /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz \
-c /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt -o /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz --inputFormat vcf \
-p /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt
The following parameters are available. Ones with "[]" are in effect:

Available Options
Required Parameters :
-i [/net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz]
-o [/net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz]
Gene Annotation
Parameters : -g [/home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz]
-r [/home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa]
--inputFormat [vcf], --checkReference, -f [refGene]
-p [/home/hmkang/seqshop/epacts/share/EPACTS/priority.txt]
-c [/home/hmkang/seqshop/epacts/share/EPACTS/codon.txt]
-u [], -d [], --se [], --si [], --outputFormat []
Other Annotation Tools : --genomeScore [], --bed [], --tabix []
Load reference genome /home/hmkang/seqshop/ref22/human_g1k_v37.chr22.fa...
DONE: 1 chromosomes and 51304566 bases are loaded.
Load codon file /home/hmkang/seqshop/epacts/share/EPACTS/codon.txt...
DONE: codon file loaded.
Load priority file /home/hmkang/seqshop/epacts/share/EPACTS/priority.txt...
DONE: 24 priority annotation types loaded.
Load gene file /home/hmkang/seqshop/epacts/share/EPACTS/hg19_gencodeV14.txt.gz...
DONE: 92627 gene loaded.
DONE: Generated frequency of each annotype type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq ].
DONE: Generated frequency of each highest priority annotation type in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq ].
Ts/Tv ratio: 2.35733
Ts observed: 2718 times; Tv observed: 1153 times.
DONE: Generated frequency of each base change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq ].
DONE: Generated frequency of each codon change in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq ].
DONE: Generated frequency of indel length in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq ].
..............................................
... Anno(tation) ...
... Xiaowei Zhan, Goncalo Abecasis ...
... Speical Thanks: ...
... Hyun Ming Kang, Yanming Li ...
... zhanxw@umich.edu ...
... Sep 2011 ...
................................................

DONE: 3871 varaints are annotated.
DONE: Generated annotation output in [ /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz ].
Annotation succeed!
mv /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
/home/hmkang/seqshop/epacts/bin/bgzip -c /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp > /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
/home/hmkang/seqshop/epacts/bin/tabix -pvcf -f /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz
rm /net/seqshop-server/hmkang/out/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz.tmp
rm /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.log /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.top.anno.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.anno.frq /net/seqshop- server/hmkang/out/assoc/snps.anno.vcf.gz.base.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.codon.frq /net/seqshop-server/hmkang/out/assoc/snps.anno.vcf.gz.indel.frq
</div>
</div>

After running annotation, you can check the annotation results. Let's look at the APOL g1 risk allele we manually examined in the SNP calling section.

$GC/bin/tabix $OUT/assoc/snps.anno.vcf.gz 22:36661906 | head -1 | cut -f 1-8

<div class="mw-collapsible mw-collapsed" style="width:400px">
View the annotation results
<div class="mw-collapsible-content" style="width:800px">
22 36661906 . A G 18 PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;FIC=-0.0092;
SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;OBS=145,191,0,0,3,2,0,0;STR=-0.040;
STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;
MQ30=0.000;SVM=1.51214;BAVGPOST=0.998;BRSQ=0.941;LDAF=0.0161;AVGPOST=1.0000;RSQ=1.0000;ERATE=0.0019;THETA=0.0013;
ANNO=Nonsynonymous:APOL1;ANNOFULL=APOL1/ENST00000397278.3:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000426053.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base971/1143:Codon324/381:Exon5/5):Exon|
APOL1/ENST00000422706.1:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/6):Exon|
APOL1/ENST00000319136.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1073/1245:Codon358/415:Exon7/7):Exon|
APOL1/ENST00000347595.7:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base662/834:Codon221/278:Exon3/3):Exon|
APOL1/ENST00000397279.4:+:Nonsynonymous(AGC/Ser/S->GGC/Gly/G:Base1025/1197:Codon342/399:Exon6/7):Exon
</div>
</div>

* What is the function of this variant?
* How many different transcript does the variant overlap with?
* How can you represent the variant in terms of amino acid changes?

== Single Variant Association Analysis ==

Let's run a single-variant association analysis using a score test.

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/single --region 22:36000000-37000000 --test b.score --pheno PHENO --run 2

After running it, you will see EPACTS output files by looking at

ls $OUT/assoc

The top association results can be viewed by

head $OUT/assoc/single.epacts.top5000

<div class="mw-collapsible mw-collapsed" style="width:400px">
View top association results
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE NS.CASE NS.CTRL AF.CASE AF.CTRL
22 36995620 36995620 22:36995620_A/G 62 36 1 0.29032 5.6717e-09 5.8262 31 31 0.51613 0.064516
22 36993088 36993088 22:36993088_G/C 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36997871 36997871 22:36997871_G/T 62 30 1 0.24194 7.3258e-07 4.9525 31 31 0.43548 0.048387
22 36987368 36987368 22:36987368_G/A 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36987861 36987861 22:36987861_A/G 62 31 1 0.25 2.0898e-06 4.7445 31 31 0.43548 0.064516
22 36985499 36985499 22:36985499_C/T 62 29 1 0.23387 5.7389e-06 4.5358 31 31 0.40323 0.064516
22 36978260 36978260 22:36978260_G/T 62 28 1 0.22581 1.5051e-05 4.3279 31 31 0.3871 0.064516
22 36998907 36998907 22:36998907_C/T 62 61 1 0.49194 0.00015557 -3.782 31 31 0.30645 0.67742
22 36667082 36667082 22:36667082_T/G 62 28 1 0.22581 0.0003506 -3.5747 31 31 0.080645 0.37097
</div>
</div>

You can look also visualize the results by QQ-plot and Manhattan plot

<div class="mw-collapsible mw-collapsed" style="width:400px">
View QQ plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.qq.pdf&</pre>
[[File:Single.epacts.qq.png]]
</div>
</div>

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Manhattan plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.epacts.mh.pdf&</pre>
[[File:Single.epacts.mh.png|900px]]
</div>
</div>

Also, you can create a zoom plot focusing on the region of interest

$EPACTS/bin/epacts-zoom --vcf $OUT/assoc/snps.anno.vcf.gz --pos 22:36995620 --prefix $OUT/assoc/single

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Zoom Plots
<div class="mw-collapsible-content" style="width:800px">
:<pre>evince $OUT/assoc/single.zoom.22.36995620.pdf&</pre>
[[File:Single.zoom.22.36995620.png]]
</div>
</div>

If you want to run EMMAX, you first need to create a kinship matrix

$EPACTS/bin/epacts-make-kin --vcf $OUT/assoc/snps.anno.vcf.gz --min-maf 0.01 --out $OUT/assoc/snps.anno.kinf --run 2 --chr 22

And run EMMAX test specifying the kinship matrix

$EPACTS/bin/epacts-single --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/emmax --region 22:36000000-37000000 --test q.emmax --pheno PHENO --run 2 --kinf $OUT/assoc/snps.anno.kinf

Then the results may look similar to previous ones.

cat $OUT/assoc/emmax.epacts.top5000

== Run Groupwise Test ==

To run group-wise test such as gene-level burden test, you need to create a marker group file. If your VCF is already annotated, you can create a group file by running

$EPACTS/bin/epacts make-group --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/snps.anno.grp --nonsyn

The group file is simply a list of marker per group name, as shown below.

cat $OUT/assoc/snps.anno.grp
APOL1 22:36655735_G/A 22:36657740_G/A 22:36661330_G/A 22:36661566_G/A 22:36661646_G/A 22:36661891_G/A 22:36661906_A/G
APOL2 22:36623731_T/C 22:36623920_G/A 22:36629466_T/A 22:36633107_C/A
APOL3 22:36537763_C/T 22:36537798_G/A 22:36556768_G/A 22:36556823_G/T
APOL4 22:36587154_G/T 22:36587202_G/A 22:36587223_G/T 22:36587346_C/T 22:36587511_C/T 22:36587704_T/C 22:36587886_C/T 22:36593714_G/A 22:36597744_A/C 22:36598049_C/G 22:36598058_T/C 22:36598081_A/T
APOL5 22:36122356_G/A 22:36122380_T/A 22:36122930_C/T 22:36123083_C/T 22:36124860_C/G
FOXRED2 22:36900271_T/C 22:36900806_A/G
MYH9 22:36681163_G/C 22:36684354_T/C 22:36710183_T/C
Metazoa_SRP 22:36711990_C/G
RBFOX2 22:36424450_A/C
RP4-633O19__A.1 22:36792162_G/A

If you have your own annotation, you can create your own burden test unit by modifying this file.

If you want to run a collapsing burden test (CMC), run the following command

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.collapse --test b.collapse --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

You can view the results by examining the output file
cat $OUT/assoc/group.collapse.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

You can run SKAT-O test in a similar way, but with a special tag

$EPACTS/bin/epacts group --ped $SS/assoc/seqshop.ped --vcf $OUT/assoc/snps.anno.vcf.gz --out $OUT/assoc/group.skato --test skat --skat-o --groupf $OUT/assoc/snps.anno.grp --pheno PHENO --run 2

And view output files

cat $OUT/assoc/group.skato.epacts

<div class="mw-collapsible mw-collapsed" style="width:400px">
View Output file
<div class="mw-collapsible-content" style="width:800px">
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
22 36655735 36661906 22:36655735-36661906_APOL1 62 0.14516 7 4 0 0.42748 1
22 36623731 36633107 22:36623731-36633107_APOL2 62 0.064516 4 1 0 0.038657 NA
22 36537763 36556823 22:36537763-36556823_APOL3 62 0.080645 4 2 0 0.40634 0
22 36587154 36598081 22:36587154-36598081_APOL4 62 0.14516 12 4 0 0.67891 0
22 36122356 36124860 22:36122356-36124860_APOL5 62 0.1129 5 2 0 0.15055 0.3
22 36900271 36900806 22:36900271-36900806_FOXRED2 NA NA 2 0 0 NA NA
22 36681163 36710183 22:36681163-36710183_MYH9 62 0.032258 3 1 0 1 NA
22 36711990 36711990 22:36711990-36711990_Metazoa_SRP NA NA 1 0 0 NA NA
22 36424450 36424450 22:36424450-36424450_RBFOX2 62 0.032258 1 1 0 1 NA
22 36792162 36792162 22:36792162-36792162_RP4-633O19__A.1 NA NA 1 0 0 NA NA
</div>
</div>

EPACTS for DIAGRAM

2012-12-14T16:04:01Z

Clement Ma: /* 5.  Report results */

= Motivation and Rationale =

[[EPACTS|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 [http://www.stat.duke.edu/~scs/Courses/Stat376/Papers/GibbsFieldEst/BiasReductionMLE.pdf (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|minimac]] or I[[IMPUTE2|mpute2]].  Starting with minimac or impute2 output:

#Download and install EPACTS
#Prepare VCF file with genotypes / dosages
#Prepare PED file with phenotypes and covariates
#Run EPACTS association pipeline
#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.

*Please download the latest version of EPACTS here:  http://www.sph.umich.edu/csg/kang/epacts/download/
*Uncompress EPACTS package to the directory you would like to install and then type the following commands

>  tar xzvf EPACTS-3.0.0.tar.gz
> cd EPACTS-3.0.0
> ./configure --prefix [INSTALL DIRECTORY]
> make
> make install

*Perform a test run by running the following command

EPACTS-3.0.0/example/test_run_epacts.sh

=== For Local Users in CSG ===

For users in CSG, EPACTS can be found here:
<pre>/net/fantasia/home/hmkang/sw/epacts2.1/</pre>
=== 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:
<pre>$ EPACTS-3.0.0/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz

$ EPACTS-3.0.0/example/1000G_dummy_pheno.ped
</pre>
<br> Run the single variant score test on the example data using this command:
<pre>$ EPACTS-3.0.0/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 &
</pre>
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 [http://genome.sph.umich.edu/wiki/Test_EPACTS_for_DIAGRAM 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 pseudo VCF format.  Download the tool from [http://www.sph.umich.edu/csg/cfuchsb/dose2vcf_v0.5.gz here]. If you used rs numbers during imputation, you can find mapping tables ready for dose2vcf [http://www.sph.umich.edu/csg/cfuchsb/mapping_rs_ALL.GIANT.phase1_release_v3.20101123.tgz here (214 Mb) ]

<br>

To run the wrapper program, use the following command
<pre>> dose2vcf/dose2vcf --dose [ .dose file ] --info [ .info file ] --out [ output file prefix ]
</pre>
For example, for imputed chromosome 1 data from FUSION:
<pre>> dose2vcf/dose2vcf --dose FUSION.GWAS.1KG.imp.chr1.dose --info FUSION.GWAS.1KG.imp.chr1.info --out out/FUSION.GWAS.1KG.imp.chr1</pre>
The expected output file is:
<pre>out/FUSION.GWAS.1KG.imp.chr1.vcf

</pre>
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 [http://vcftools.sourceforge.net/docs.html#concat 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
<pre>bgzip input.vcf ## this command will produce input.vcf.gz
tabix -pvcf -f input.vcf.gz ## this command will produce input.vcf.gz.tbi
</pre>
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.<br>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 [http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml 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:

#Family ID
#Individual ID
#Paternal ID
#Maternal ID
#Sex (1=male; 2=female; other=unknown)
#Phenotype (1 = control; 2 = case)

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

#QT
#AGE<br>

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.
<pre>#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
</pre>
Alternatively, you can prepare a PED file without a header, and include a corresponding DAT file describing the column headers.
<pre>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

</pre>
The corresponding DAT file is:
<pre>A DISEASE
T QT
C AGE</pre>
Key:  A = binary trait; T = quantitative trait; C = covariate<br>

== 4.  Run EPACTS association pipeline ==

For detailed description of options, use:
<pre>EPACTS-3.0.0/epacts single -man
</pre>
<br>

=== Primary analyses (without BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' separate association analyses to be completed '''without adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#1._Typical_DIAGRAM_analysis_using_existing_association_pipeline A.  Typical DIAGRAM analysis using existing association pipeline]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#2._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression B.  Analysis of low frequency variants using Firth bias-corrected logistic regression]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

=== Analysis for QC [please submit as soon as the analysis is complete] ===

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.

<br>

{| width="1650" border="1" cellpadding="1" cellspacing="1" align="center"
|-
| align="center" | '''Association Analysis'''
| align="center" | '''Statistical Test'''
| align="center" | '''Subset of SNPs'''
| align="center" | '''Output File Type'''
| align="center" | '''Output Filename Format'''
|-
| [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#3._Analysis_of_chromosome_20_using_logistic_regression_score_test C.  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
|}

<br>

=== Secondary analyses (with BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' secondary analyses''' adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#4._Typical_DIAGRAM_analysis_using_existing_association_pipeline_.28with_BMI.29 D.  Typical DIAGRAM analysis using existing association pipeline (with BMI)]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#5._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression_.28with_BMI.29 E.  Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI)]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

'''Please send the 2 Primary analyses and the QC analysis when complete.'''

<br>

==== Filename key: ====

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

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

<br>

----

=== A. Typical DIAGRAM analysis using existing association pipeline<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software.   [[Image:1000Genomes_march2012_imputation_analysis_plan_08312012_v2.pdf]] (Updated Dec 14, 2012)

For frequently asked questions regarding the file format, please see:  [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications]

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

As preparation for the Firth test analysis, we encourage you to analyze the data using the Wald test first, since it is computationally much faster.  This will be a good way to check if your VCF and PED files for every chromosome are correctly formatted for EPACTS and resolve any problems you may have with your imputation or input files.
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== B. 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:
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== C. 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:
<pre>EPACTS-3.0.0/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
</pre>
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).

=== D. Typical DIAGRAM analysis using existing association pipeline (with BMI)<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software including BMI adjustment.

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

<br>

=== E. Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI) ===

Again use the Firth test on EPACTS for your analysis with BMI

== 5.  Report results<br> ==
For '''analysis 1''', please follow the following results file guidelines:   [[Image:1000Genomes_march2012_imputation_analysis_plan_08312012_v2.pdf]] (Updated Dec 14, 2012)

For frequently asked questions regarding the file format, please see:  [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications]

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

#'''Firth test (no BMI): ''' DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz
#'''Score test (chr 20, no BMI):  '''DIAGRAMv4_iSNPs_XXX_1000G_KKK_SCR_YYY_ZZZ.epacts.gz

<br>

The FTP hostname is:  '''ftp.broadinstitute.org'''.  Please place your files into to the /incoming/ directory.

Here's an example score test .epacts file
<pre>$ 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

</pre>
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 =

== EPACTS installation errors ==

== Errors when running EPACTS ==

=== Rscript execution error: No such file or directory ===
<pre>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

</pre>
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.

<br>

=== 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:

<br>
<pre>#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

</pre>
<br> 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".

== Results FIle Clarifications ==

=== 1. How do I code the INDEL variant names and alleles? ===

Please use the variant name and the allele name directly from IMPUTE or minimac. Please do NOT recode variant names or alleles. We will do this step in the analysis for consistency.

ACTION IF YOU HAVE UPLOADED YOUR FILE: If you have recoded your INDEL alleles, please tell us so we can remove your file and let us know when you can reupload with the original variable and allele names.

=== <br>2. The document asks for the number of homozygotes and heterozygotes in case and control. How do I get this from my data? Is this relevant for imputed data? ===

These numbers were relevant to genotyped data but not for imputed data. We didn't intend to ask for this. To retain the same file format between results already submitted and those to be submitted please retain the columns with a "." for the value.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip these columns.

=== 3. For the "Imputed" variable, what does imputed mean in the context of the data output from MACH and IMPUTE? ===

This is a hold over from the last round of analysis where we asked for results separately from genotyped SNPs and imputed SNPs and wanted to distinguish between the two. We will use r2_hat or info measures to estimate the accuracy of the genotypes. This column will be retained for consistency with files already submitted but should be filled in with "." or "1". It will not be used in the analysis.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip this column.

EPACTS for DIAGRAM

2012-12-14T16:03:10Z

Clement Ma: /* 5.  Report results */

= Motivation and Rationale =

[[EPACTS|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 [http://www.stat.duke.edu/~scs/Courses/Stat376/Papers/GibbsFieldEst/BiasReductionMLE.pdf (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|minimac]] or I[[IMPUTE2|mpute2]].  Starting with minimac or impute2 output:

#Download and install EPACTS
#Prepare VCF file with genotypes / dosages
#Prepare PED file with phenotypes and covariates
#Run EPACTS association pipeline
#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.

*Please download the latest version of EPACTS here:  http://www.sph.umich.edu/csg/kang/epacts/download/
*Uncompress EPACTS package to the directory you would like to install and then type the following commands

>  tar xzvf EPACTS-3.0.0.tar.gz
> cd EPACTS-3.0.0
> ./configure --prefix [INSTALL DIRECTORY]
> make
> make install

*Perform a test run by running the following command

EPACTS-3.0.0/example/test_run_epacts.sh

=== For Local Users in CSG ===

For users in CSG, EPACTS can be found here:
<pre>/net/fantasia/home/hmkang/sw/epacts2.1/</pre>
=== 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:
<pre>$ EPACTS-3.0.0/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz

$ EPACTS-3.0.0/example/1000G_dummy_pheno.ped
</pre>
<br> Run the single variant score test on the example data using this command:
<pre>$ EPACTS-3.0.0/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 &
</pre>
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 [http://genome.sph.umich.edu/wiki/Test_EPACTS_for_DIAGRAM 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 pseudo VCF format.  Download the tool from [http://www.sph.umich.edu/csg/cfuchsb/dose2vcf_v0.5.gz here]. If you used rs numbers during imputation, you can find mapping tables ready for dose2vcf [http://www.sph.umich.edu/csg/cfuchsb/mapping_rs_ALL.GIANT.phase1_release_v3.20101123.tgz here (214 Mb) ]

<br>

To run the wrapper program, use the following command
<pre>> dose2vcf/dose2vcf --dose [ .dose file ] --info [ .info file ] --out [ output file prefix ]
</pre>
For example, for imputed chromosome 1 data from FUSION:
<pre>> dose2vcf/dose2vcf --dose FUSION.GWAS.1KG.imp.chr1.dose --info FUSION.GWAS.1KG.imp.chr1.info --out out/FUSION.GWAS.1KG.imp.chr1</pre>
The expected output file is:
<pre>out/FUSION.GWAS.1KG.imp.chr1.vcf

</pre>
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 [http://vcftools.sourceforge.net/docs.html#concat 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
<pre>bgzip input.vcf ## this command will produce input.vcf.gz
tabix -pvcf -f input.vcf.gz ## this command will produce input.vcf.gz.tbi
</pre>
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.<br>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 [http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml 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:

#Family ID
#Individual ID
#Paternal ID
#Maternal ID
#Sex (1=male; 2=female; other=unknown)
#Phenotype (1 = control; 2 = case)

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

#QT
#AGE<br>

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.
<pre>#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
</pre>
Alternatively, you can prepare a PED file without a header, and include a corresponding DAT file describing the column headers.
<pre>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

</pre>
The corresponding DAT file is:
<pre>A DISEASE
T QT
C AGE</pre>
Key:  A = binary trait; T = quantitative trait; C = covariate<br>

== 4.  Run EPACTS association pipeline ==

For detailed description of options, use:
<pre>EPACTS-3.0.0/epacts single -man
</pre>
<br>

=== Primary analyses (without BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' separate association analyses to be completed '''without adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#1._Typical_DIAGRAM_analysis_using_existing_association_pipeline A.  Typical DIAGRAM analysis using existing association pipeline]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#2._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression B.  Analysis of low frequency variants using Firth bias-corrected logistic regression]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

=== Analysis for QC [please submit as soon as the analysis is complete] ===

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.

<br>

{| width="1650" border="1" cellpadding="1" cellspacing="1" align="center"
|-
| align="center" | '''Association Analysis'''
| align="center" | '''Statistical Test'''
| align="center" | '''Subset of SNPs'''
| align="center" | '''Output File Type'''
| align="center" | '''Output Filename Format'''
|-
| [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#3._Analysis_of_chromosome_20_using_logistic_regression_score_test C.  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
|}

<br>

=== Secondary analyses (with BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' secondary analyses''' adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#4._Typical_DIAGRAM_analysis_using_existing_association_pipeline_.28with_BMI.29 D.  Typical DIAGRAM analysis using existing association pipeline (with BMI)]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#5._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression_.28with_BMI.29 E.  Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI)]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

'''Please send the 2 Primary analyses and the QC analysis when complete.'''

<br>

==== Filename key: ====

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

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

<br>

----

=== A. Typical DIAGRAM analysis using existing association pipeline<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software.   [[Image:1000Genomes_march2012_imputation_analysis_plan_08312012_v2.pdf]] (Updated Dec 14, 2012)

For frequently asked questions regarding the file format, please see:  [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications]

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

As preparation for the Firth test analysis, we encourage you to analyze the data using the Wald test first, since it is computationally much faster.  This will be a good way to check if your VCF and PED files for every chromosome are correctly formatted for EPACTS and resolve any problems you may have with your imputation or input files.
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== B. 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:
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== C. 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:
<pre>EPACTS-3.0.0/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
</pre>
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).

=== D. Typical DIAGRAM analysis using existing association pipeline (with BMI)<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software including BMI adjustment.

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

<br>

=== E. Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI) ===

Again use the Firth test on EPACTS for your analysis with BMI

== 5.  Report results<br> ==
For '''analysis 1''', please follow the following results file guidelines:   [[Image:1000Genomes_march2012_imputation_analysis_plan_08312012_v2.pdf]] (Updated Dec 14, 2012)

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

#'''Firth test (no BMI): ''' DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz
#'''Score test (chr 20, no BMI):  '''DIAGRAMv4_iSNPs_XXX_1000G_KKK_SCR_YYY_ZZZ.epacts.gz

<br>

The FTP hostname is:  '''ftp.broadinstitute.org'''.  Please place your files into to the /incoming/ directory.

Here's an example score test .epacts file
<pre>$ 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

</pre>
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 =

== EPACTS installation errors ==

== Errors when running EPACTS ==

=== Rscript execution error: No such file or directory ===
<pre>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

</pre>
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.

<br>

=== 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:

<br>
<pre>#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

</pre>
<br> 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".

== Results FIle Clarifications ==

=== 1. How do I code the INDEL variant names and alleles? ===

Please use the variant name and the allele name directly from IMPUTE or minimac. Please do NOT recode variant names or alleles. We will do this step in the analysis for consistency.

ACTION IF YOU HAVE UPLOADED YOUR FILE: If you have recoded your INDEL alleles, please tell us so we can remove your file and let us know when you can reupload with the original variable and allele names.

=== <br>2. The document asks for the number of homozygotes and heterozygotes in case and control. How do I get this from my data? Is this relevant for imputed data? ===

These numbers were relevant to genotyped data but not for imputed data. We didn't intend to ask for this. To retain the same file format between results already submitted and those to be submitted please retain the columns with a "." for the value.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip these columns.

=== 3. For the "Imputed" variable, what does imputed mean in the context of the data output from MACH and IMPUTE? ===

This is a hold over from the last round of analysis where we asked for results separately from genotyped SNPs and imputed SNPs and wanted to distinguish between the two. We will use r2_hat or info measures to estimate the accuracy of the genotypes. This column will be retained for consistency with files already submitted but should be filled in with "." or "1". It will not be used in the analysis.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip this column.

EPACTS for DIAGRAM

2012-12-14T16:02:58Z

Clement Ma: /* A. Typical DIAGRAM analysis using existing association pipeline */

= Motivation and Rationale =

[[EPACTS|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 [http://www.stat.duke.edu/~scs/Courses/Stat376/Papers/GibbsFieldEst/BiasReductionMLE.pdf (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|minimac]] or I[[IMPUTE2|mpute2]].  Starting with minimac or impute2 output:

#Download and install EPACTS
#Prepare VCF file with genotypes / dosages
#Prepare PED file with phenotypes and covariates
#Run EPACTS association pipeline
#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.

*Please download the latest version of EPACTS here:  http://www.sph.umich.edu/csg/kang/epacts/download/
*Uncompress EPACTS package to the directory you would like to install and then type the following commands

>  tar xzvf EPACTS-3.0.0.tar.gz
> cd EPACTS-3.0.0
> ./configure --prefix [INSTALL DIRECTORY]
> make
> make install

*Perform a test run by running the following command

EPACTS-3.0.0/example/test_run_epacts.sh

=== For Local Users in CSG ===

For users in CSG, EPACTS can be found here:
<pre>/net/fantasia/home/hmkang/sw/epacts2.1/</pre>
=== 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:
<pre>$ EPACTS-3.0.0/example/1000G_exome_chr20_example_softFiltered.calls.vcf.gz

$ EPACTS-3.0.0/example/1000G_dummy_pheno.ped
</pre>
<br> Run the single variant score test on the example data using this command:
<pre>$ EPACTS-3.0.0/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 &
</pre>
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 [http://genome.sph.umich.edu/wiki/Test_EPACTS_for_DIAGRAM 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 pseudo VCF format.  Download the tool from [http://www.sph.umich.edu/csg/cfuchsb/dose2vcf_v0.5.gz here]. If you used rs numbers during imputation, you can find mapping tables ready for dose2vcf [http://www.sph.umich.edu/csg/cfuchsb/mapping_rs_ALL.GIANT.phase1_release_v3.20101123.tgz here (214 Mb) ]

<br>

To run the wrapper program, use the following command
<pre>> dose2vcf/dose2vcf --dose [ .dose file ] --info [ .info file ] --out [ output file prefix ]
</pre>
For example, for imputed chromosome 1 data from FUSION:
<pre>> dose2vcf/dose2vcf --dose FUSION.GWAS.1KG.imp.chr1.dose --info FUSION.GWAS.1KG.imp.chr1.info --out out/FUSION.GWAS.1KG.imp.chr1</pre>
The expected output file is:
<pre>out/FUSION.GWAS.1KG.imp.chr1.vcf

</pre>
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 [http://vcftools.sourceforge.net/docs.html#concat 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
<pre>bgzip input.vcf ## this command will produce input.vcf.gz
tabix -pvcf -f input.vcf.gz ## this command will produce input.vcf.gz.tbi
</pre>
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.<br>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 [http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml 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:

#Family ID
#Individual ID
#Paternal ID
#Maternal ID
#Sex (1=male; 2=female; other=unknown)
#Phenotype (1 = control; 2 = case)

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

#QT
#AGE<br>

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.
<pre>#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
</pre>
Alternatively, you can prepare a PED file without a header, and include a corresponding DAT file describing the column headers.
<pre>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

</pre>
The corresponding DAT file is:
<pre>A DISEASE
T QT
C AGE</pre>
Key:  A = binary trait; T = quantitative trait; C = covariate<br>

== 4.  Run EPACTS association pipeline ==

For detailed description of options, use:
<pre>EPACTS-3.0.0/epacts single -man
</pre>
<br>

=== Primary analyses (without BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' separate association analyses to be completed '''without adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#1._Typical_DIAGRAM_analysis_using_existing_association_pipeline A.  Typical DIAGRAM analysis using existing association pipeline]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#2._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression B.  Analysis of low frequency variants using Firth bias-corrected logistic regression]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

=== Analysis for QC [please submit as soon as the analysis is complete] ===

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.

<br>

{| width="1650" border="1" cellpadding="1" cellspacing="1" align="center"
|-
| align="center" | '''Association Analysis'''
| align="center" | '''Statistical Test'''
| align="center" | '''Subset of SNPs'''
| align="center" | '''Output File Type'''
| align="center" | '''Output Filename Format'''
|-
| [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#3._Analysis_of_chromosome_20_using_logistic_regression_score_test C.  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
|}

<br>

=== Secondary analyses (with BMI) [please submit as soon as the analysis is complete] ===

There are '''2''' secondary analyses''' adjusting for BMI'''.

{| width="1650" border="1" align="left" cellpadding="1" cellspacing="1"
|-
! scope="col" |
Association Analysis

! scope="col" |
Statistical Test

! scope="col" |
Subset of SNPs

! scope="col" |
Output File Type

! scope="col" |
Output Filename Format

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#4._Typical_DIAGRAM_analysis_using_existing_association_pipeline_.28with_BMI.29 D.  Typical DIAGRAM analysis using existing association pipeline (with BMI)]

<br>

|
Wald or likelihood ratio

|
All SNPs with

MAC >= 1

|
Custom file

based on DIAGRAM format

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_TTT_YYY_ZZZ.txt

<br>

|-
|
[http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#5._Analysis_of_low_frequency_variants_using_Firth_bias-corrected_logistic_regression_.28with_BMI.29 E.  Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI)]

<br>

|
Firth bias-corrected

|
SNPs with

200 >= MAC >= 1

|
EPACTS output file

|
DIAGRAMv4_iSNPs_XXX_adjBMI_1000G_KKK_FBC_YYY_ZZZ.epacts.gz

<br>

|}

<br>

<br>

<br>

<br>

<br>

<br>

<br>

'''Please send the 2 Primary analyses and the QC analysis when complete.'''

<br>

==== Filename key: ====

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

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

<br>

----

=== A. Typical DIAGRAM analysis using existing association pipeline<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software.   [[Image:1000Genomes_march2012_imputation_analysis_plan_08312012_v2.pdf]] (Updated Dec 14, 2012)

For frequently asked questions regarding the file format, please see:  [http://genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications genome.sph.umich.edu/wiki/EPACTS_for_DIAGRAM#Results_FIle_Clarifications]

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

As preparation for the Firth test analysis, we encourage you to analyze the data using the Wald test first, since it is computationally much faster.  This will be a good way to check if your VCF and PED files for every chromosome are correctly formatted for EPACTS and resolve any problems you may have with your imputation or input files.
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== B. 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:
<pre>EPACTS-3.0.0/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
</pre>
'''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")!

=== C. 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:
<pre>EPACTS-3.0.0/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
</pre>
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).

=== D. Typical DIAGRAM analysis using existing association pipeline (with BMI)<br> ===

This is the typical DIAGRAM analysis using your current association pipeline and software including BMI adjustment.

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

<br>

=== E. Analysis of low frequency variants using Firth bias-corrected logistic regression (with BMI) ===

Again use the Firth test on EPACTS for your analysis with BMI

== 5.  Report results<br> ==
For '''analysis 1''', please follow the following results file guidelines:   [[Image:1000Genomes march2012 imputation analysis plan 08312012.pdf]]

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

#'''Firth test (no BMI): ''' DIAGRAMv4_iSNPs_XXX_1000G_KKK_FBC_YYY_ZZZ.epacts.gz
#'''Score test (chr 20, no BMI):  '''DIAGRAMv4_iSNPs_XXX_1000G_KKK_SCR_YYY_ZZZ.epacts.gz

<br>

The FTP hostname is:  '''ftp.broadinstitute.org'''.  Please place your files into to the /incoming/ directory.

Here's an example score test .epacts file
<pre>$ 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

</pre>
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 =

== EPACTS installation errors ==

== Errors when running EPACTS ==

=== Rscript execution error: No such file or directory ===
<pre>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

</pre>
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.

<br>

=== 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:

<br>
<pre>#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

</pre>
<br> 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".

== Results FIle Clarifications ==

=== 1. How do I code the INDEL variant names and alleles? ===

Please use the variant name and the allele name directly from IMPUTE or minimac. Please do NOT recode variant names or alleles. We will do this step in the analysis for consistency.

ACTION IF YOU HAVE UPLOADED YOUR FILE: If you have recoded your INDEL alleles, please tell us so we can remove your file and let us know when you can reupload with the original variable and allele names.

=== <br>2. The document asks for the number of homozygotes and heterozygotes in case and control. How do I get this from my data? Is this relevant for imputed data? ===

These numbers were relevant to genotyped data but not for imputed data. We didn't intend to ask for this. To retain the same file format between results already submitted and those to be submitted please retain the columns with a "." for the value.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip these columns.

=== 3. For the "Imputed" variable, what does imputed mean in the context of the data output from MACH and IMPUTE? ===

This is a hold over from the last round of analysis where we asked for results separately from genotyped SNPs and imputed SNPs and wanted to distinguish between the two. We will use r2_hat or info measures to estimate the accuracy of the genotypes. This column will be retained for consistency with files already submitted but should be filled in with "." or "1". It will not be used in the analysis.

ACTION IF YOU HAVE UPLOADED YOUR FILE: No action. You do not need to redo the file. We will skip this column.