How to compare groups using WGS data? We have whole genome sequencing data for patients (not-cancer) (n=60) and for healthy controls (n=20). The sequencing centre has provided us with the best practice bioinformatics analyses including reads mapping (.BAM) and variant calling using GATK (.vcf) as well as annotation (annotated .vcf and .gVCF). 

What should be our next steps? We are interested to see if there are any differences (global and/or specific) between the groups.
I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.

An overview of Hail can be found here:

http://blog.cloudera.com/blog/2017/05/hail-scalable-genomics-analysis-with-spark/

A tutorial on association testing can be found here:

https://hail.is/hail/tutorial.html#Association-testing