This is a suit of GWAS workflows for CHIMGEN project built on NEXTFLOW framework.
I built the whole pipeline on NEXTFLOW framework. And NEXTFLOW can be installed by conda. To get started, you should clone this repo firstly.
git clone https://github.com/Jianhua-Wang/GWAS-pipeline.git
then set up the conda environment:
# create conda environment
cd bin
conda env create -f gwascondaenv.yml
# download qctool which is too big to upload to GitHub
wget -c https://www.well.ox.ac.uk/~gav/resources/qctool_v2.0.1-Ubuntu14.04-x86_64.tgz
tar -zxvf qctool_v2.0.1-Ubuntu14.04-x86_64.tgz
mv ./qctool_v2.0.1-Ubuntu14.04-x86_64/qctool ./
rm -rf qctool_v2.0.1-Ubuntu14.04-x86_64*
and use the commands below to activate/deactivate the environment:
conda activate gwas
conda deactivate
This usage of this pipeline is quite simple, for example:
nextflow run qc.nf -c qc.config
qc.nf
is the NEXTFLOW script for GWAS quality control and qc.config
is the configuration file that stores the parameters of the pipeline. NEXTFLOW will genere a work
folder in work directory for debuging and storing intermediate files. You can remove the Intermediate files by manually or using nextflow clean -f
.
This is a pipeline of genotype data quality assessment and control, including variant-level QC and individual-level QC. These steps are critical for the following imputation and association analysis, could reduce the false positive of GWAS results.
I referred to the Nature Protocols from Anderson et al. and a NEXTFLOW workflow from H3ABioNet. I didn't describe the details of every step in this tutorial, but it's very easy to understand the principles on the Nature Protocols.
PLINK binary .bed, .bim, and .fam file (10 sample)
There are various types of genotype formats, such as PLINK bed, gen, vcf, bgen, even pgen. You can find the description of them on PLINK's website. In this part, we mainly use PLINK bed format, which contains a bed file for genotype, a bim file for variant information, and a fam file for individual information, with the same prefix.
This pipeline does not only output the clean data but also generate tables of filtered variants and individuals and figures of data assessment. I put the tables and figures in a HTML report to read easily.
Output file | Description |
---|---|
test-nd-c-c.bim | bim file of clean data |
test-nd-c-c.bed | bed file of clean data |
test-nd-c-c.fam | fam file of clean data |
test.dups | duplicated variants |
test.badsex | sample with ambiguous sex |
test-nd-c.irem | removed sample due to missing rate |
test-nd-fail_IBD.txt | related sample |
test-nd-c-fail_het.txt | sample with extreme high or low heterozygous rate |
test-nd-c-c.irem | sample removed in phased 3 |
test-GWAS-QC_report.html | report of QC |
test_runtime.html | report of workflow running |
test_timeline.html | report of workflow running |
To run GWAS QC pipeline on the test data in input
folder, just activate conda environment and use:
conda activate gwas
nextflow run qc.nf -c qc.config
One can adjust the parameters such as prefix of input PLINK file, path of output file, and the cutoff of MAF modifying the qc.config
file:
params {
work_dir = "/$PWD"
input_dir = "${params.work_dir}/input"
output_dir = "${params.work_dir}/output/qc"
input_pat = "test"
pca_ref_dir = "./pca_ref/"
pca_ref_pat = "hapmap3_pca"
cut_maf = 0.0001
cut_mind = 0.03
cut_geno = 0.05
cut_hwe = 0.00001
pi_hat = 0.1875
f_lo_male = 0.8
f_hi_female = 0.2
times_of_meanhet = 3
}
or use the command line directly, like:
nextflow run qc.nf -c qc.config --output_dir qc1 --cut_maf 0.01
This command changed the path of output file. You can find a new folder named qc1
was created in output
folder and the clean data is with MAF >= 0.01.
⭐ The bin
folder is the default environmental variable of NEXTFLOW, so I just put the stable version of used tools in bin
folder instead of using a script to set up.
The stdout of NEXTFLOW looks like:
(gwas) ➜ GWAS-pipeline git:(master) ✗ nextflow run qc.nf -c qc.config
N E X T F L O W ~ version 0.30.1
Launching `qc.nf` [elated_easley] - revision: 16c33dbc96
[warm up] executor > local
[38/981164] Submitted process > identifyIndivDiscSexinfo (1)
[ea/0f8c22] Submitted process > getDuplicateMarkers (1)
[21/126daf] Submitted process > showHWEStats (1)
[f9/69193c] Submitted process > removeDuplicateSNPs (1)
[aa/ce6083] Submitted process > getInitMAF (1)
[2a/0c7329] Submitted process > removeQCPhase1 (1)
[99/6104a7] Submitted process > generateIndivMissingnessPlot (1)
[bb/9f881b] Submitted process > generateSnpMissingnessPlot (1)
[fe/33324b] Submitted process > showInitMAF (1)
[ff/38d030] Submitted process > calculateSampleHeterozygosity (1)
[a7/f6167b] Submitted process > pruneForIBD (1)
[12/a1cd0a] Submitted process > compPCA (1)
[b0/4c675a] Submitted process > getBadIndivsMissingHet (1)
[b7/c9ea9e] Submitted process > generateMissHetPlot (1)
[fa/1ef8f2] Submitted process > findRelatedIndiv (1)
[ee/217c0c] Submitted process > removeQCIndivs (1)
[56/6f17db] Submitted process > drawPCA (1)
[c5/ab2f2d] Submitted process > calculateMaf (1)
[7f/307d0a] Submitted process > generateHwePlot (1)
[2b/7eb74b] Submitted process > generateMafPlot (1)
[b3/fdb024] Submitted process > produceReports (1)
The output report is called //f/jianhua/jianhua_pipeline/GWAS-pipeline/output/qc/test-GWAS-QC_report.html
There are 21 steps in this QC pipeline and I wrote them in NEXTFLOW and lots of Python scripts, so it's not possible to explain each step clearly here. To understand the details in the pipeline, you need to be proficient in PLINK, NEXTFLOW and Python.
Here I use a simple example to demonstrate how the NEXTFLOW works.
process calculateMaf {
input:
file(plinks) from qc3A_ch
output:
file "${base}.frq" into maf_plot_ch
file "${base}.hwe" into hwe_scores_ch
script:
base = plinks[0].baseName
out = base.replace(".","_")
"""
plink --bfile $base --hardy --freq --out $out
"""
}
In NEXTFLOW, a single step is defined as a process
. Every process
has three sections: input, output, and script. The input file of every step is loaded in input section and the output file are pull out to other channels using wildcard. In script section, I used --hardy
and --freq
to calculate the HWE and MAF.
Genotype imputation is a statistical technique that is often used to increase the power and resolution of genetic association studies. Imputation methods work by using haplotype patterns in a reference panel, such as the 1000 Genomes Project (1000G) to predict unobserved genotypes in a study dataset.
I used the reference panel constructed by IMPUTE2. You can download the files using:
cd impute_ref
bash 01_prepare_reference.sh
The input files of imputation are the output files of QC pipeline which are PLINK bed format.
The original output of IMPUTE2 is .gen format which is readable and takes up lots of space. So I convert the .gen format to the .bgen format which is far small then .gen format and very fast to read using qctool.
(gwas) ➜ imputation git:(master) ✗ ll
total 180M
-rw-rw-r-- 2 jianhua jianhua 177M 9月 2 19:56 test-nd-c-c.bgen
-rw-rw-r-- 1 jianhua jianhua 2.9M 9月 2 19:56 test-nd-c-c_runtime.html
-rw-rw-r-- 1 jianhua jianhua 15K 9月 2 19:56 test-nd-c-c_timeline.html
Same as the usage of QC pipeline, it's very easy to run this pipeline:
nextflow run impute_with_onephased_prephasing.nf -c impute_with_onephased_prephasing.config
the parameters in configuration file are most fixed parameters:
params {
work_dir = "/$PWD"
output_dir = "${params.work_dir}/output/imputation"
input_dir = "${params.work_dir}/input"
input_pat = "test"
chromosomes_List = [21,22]
chromosomeSizesFile = "./hg19_chromosomes_size/hg19_chromosomes_size.txt"
map_dir = "./impute_ref/1000GP_Phase3/"
map_pattern = "genetic_map_chr%s_combined_b37.txt"
ref_panel_dir = "./impute_ref/1000GP_Phase3/"
ref_hap_pattern = "1000GP_Phase3_chr%s.hap.gz"
ref_leg_pattern = "1000GP_Phase3_chr%s.legend.gz"
ref_sample = "1000GP_Phase3.sample"
maxForks_shapeit = 10
maxForks_impute2 = 20
}
In the test example, I only imputed chr21 and chr22. You can replace [21,22] with 1..22 for autosomes imputation.
💥 I split file in every 5Mb bin for parallel running. You will get a error if there is no variant in this bin when using IMPUTE2. In order to ensure success, I have ignored such errors, the stdout looks like:
[1e/a15433] NOTE: Missing output file(s) `chr21-5000001-10000000.imputed` expected by process `impute2 (2)` -- Error is ignored
[e8/0b9125] NOTE: Missing output file(s) `chr21-10000001-15000000.imputed` expected by process `impute2 (3)` -- Error is ignored
[0d/448b04] NOTE: Missing output file(s) `chr22-50000001-51304566.imputed` expected by process `impute2 (21)` -- Error is ignored
As the graph shows, I split the input plink format into .gen by chromosomes, then impute under every 5,000,000 BP bin and merge the results of chromosomes, finally, convert the gen format files into bgen format.
This parallel strategy is easy to implement in NEXTFLOW, however, imputation will produce many large files.
As for association test, we reviewed lots of GWAS quantitative traits and summarize the association model they used. The most used and the earliest is generalized linear model, which is applied in PLINK and hail. Hail is a python module developed by Board institute. Different with GLM, mixed linear model involves the concept of random effect which is in estimated by genetic relationship matrix. And GCTA is one of the representative tools. As the imputation becoming mainstream, many GWAS use SNPTEST which was also developed by Marchini. It's based on a missing data likelihood theory, fully accounting for the uncertainty in genotypes.
- results of imputation (bgen format)
- information of sample: ID, phenotypes, covariates (sample format)
cat test-nd-c-c.sample
ID_1 ID_2 missing tiv age height sex
0 0 0 P C C D
1 1 0 1 21 170 1
2 2 0 2 22 171 2
3 3 0 3 23 172 1
4 4 0 4 24 173 2
5 5 0 5 25 174 2
6 6 0 6 26 175 1
7 7 0 7 27 176 1
8 8 0 8 28 177 2
9 9 0 9 29 178 1
10 10 0 10 30 179 2
Above is the sample file I used in the example. This format has a two-line head and the first three are always ID_1, ID_2 and missing. Column missing
is missing rate of every individual and I fill it with zero because SNPTEST will calculate by itself instead of using this column. tiv
is the interested phenotype while the second line is "P" for "Phenotype". And I used age, height, and sex as covariates. In second head, continuous covariates are represented by 'C' while discrete covariates are represented by 'D'.
Summary statistics and a PDF file with Manhattan plot and qq plot.
conda activate gwas
nextflow run snptest_frequentist.nf -c snptest_frequentist.config
conda deactivate
- In PCA of QC section, I used hapmap3 as the background. So the PCA is based on the overlapped variants between hapmap3 and input data. Sometimes it will report an error if the overlapped variants has inconsistent alleles, so make sure your variants are all forward strand, which is also important in imputation.
- In imputation section, I only imputed the variants in autosomes because the association analysis of chrX is complicated, although imputation on chrX is possible.
- You will got nothing on the manhattan plot and qq plot, and the P-value in results are all NA, cause SNPTEST requires at least 100 sample but we only have 10 in the test data. You will see the correct results with true data.
- Anderson, C. A., Pettersson, F. H., Clarke, G. M., Cardon, L. R., Morris, A. P., & Zondervan, K. T. (2010). Data quality control in genetic case-control association studies. Nature Protocols, 5(9), 1564–1573. https://doi.org/10.1038/nprot.2010.116
- Van Leeuwen, E. M., Kanterakis, A., Deelen, P., Kattenberg, M. V., Slagboom, P. E., De Bakker, P. I. W., … Hottenga, J. J. (2015). Population-specific genotype imputations using minimac or IMPUTE2. Nature Protocols, 10(9), 1285–1296. https://doi.org/10.1038/nprot.2015.077
- Clarke, G. M., Anderson, C. A., Pettersson, F. H., Cardon, L. R., Morris, A. P., & Zondervan, K. T. (2011). Basic statistical analysis in genetic case-control studies. Nature Protocols, 6(2), 121–133. https://doi.org/10.1038/nprot.2010.182