Step 2. Performing KnockOffGWAS

To perform KnockOffGWAS it is necessary that the data has been prepared as described in section Step 1. File Preparation.

You should have already have created a script called set_dirs.sh saved it in your analysis directory. This script should look something like:

# Run this file using "source set_dirs.sh" to set the following variable
DATA=/nobackup/proj/my_account/data

To run the analysis there is a script called run_knockoff_gwas.sh with the following parameters.

Command-Line Parameters for run_knockoff_gwas.sh

  1. Chromosome number.

  2. Path and filename prefix of the data.

  3. The phenotype name to give to the phenotype data.

  4. The false discovery rate (FDR).

  5. The directory name used to store the results. This directory will be created automatically by the pipeline.

As the analysis will take a long time is is best to run the analysis as an array job on an HPC machine. This can be done by creating a script in the hpc directory called my_analysis.sh (or something) and should look something like the following:

#!/bin/bash
#SBATCH --partition=default_free
#SBATCH --account=my_account
#SBATCH --cpus-per-task=1
#SBATCH --mem=100GB
#SBATCH --array=1-22                      # Run tasks for given chromosomes
#SBATCH --output=slurm_anal_%a.out

# Load modules

module load BCFtools/1.22-GCC-13.3.0
module load PLINK/1.9b_6.21-x86_64
module load R/4.5.1-gfbf-2024a

# Set dirs
source ./set_dirs.sh

date
echo "Running on $HOSTNAME PBC analysis"

../new_knockoffgwas_pipeline/run_knockoff_gwas.sh $SLURM_ARRAY_TASK_ID $DATA/mydata pbc 0.1 results_d25_w3_FDR20

date

As before, this will need to be updated for the requirements of the HPC machine that you are using. The script requires BCFTools, Plink version 1.9 and R, so these must be loaded.

The chromosome is set using the task job number given by $SLURM_ARRAY_TASK_ID.

The path and filename prefix of the data is given next by $DATA/mydata. The data should be formatted as described in the previous section; see Initial data preparation.

The phenotype name is given next, which is set here to pbc for the Primary Biliary Cholangitis (PBC) dataset.

The results directory is set to results.

In this case it is set to pbc for the Primary Biliary Cholangitis (PBC) dataset.

The False Discovery Rate (FDR) is set to 0.1 for the KnockOffGWAS analysis.

Run the analysis as an array job on the HPC with the following command:

sbatch hpc/my_analysis.sh

or whatever is appropriate for the HPC machine you are using.

Results for the analysis will be stored in the results with a statistics file for each chromosome 1-22 and seven different resolutions 0-6, for example: results_chr1_chr1_res0_stats.txt. There are also accompanying files for regions passing the FDR criteria, set to 10% above, for example: results_chr1_chr1_res0_discoveries.txt.

Goodness of fit plots (PNG files) are created and stored in results/knockoffs_gof_plots.

Removing temporary files

A number of files are produced during the data preprocessing and during analysis calculations. These files can be removed by running a script, which also restores the .fam file which needed to have its format altered. If you remove the temporary files and wish to rerun the analysis for any reason you will need to rerun the data preprocessing scripts also.

The script, run_remove_temporary_files.sh, can be ran as follows:

source set_dirs.sh
../new_knockoffgwas_pipeline/run_remove_temporary_files.sh $DATA/mydata

Where the parameter is the path and filename of the data used in the analysis.