• Welcome
• Definitions and workflow
• Download
• Directory architecture
• Quality control of genetic data
• Input data required for RV-EXCALIBER
• Essential dependencies
• Retained dependency-generated files
• Definitions of script inputs
• Preparing the RV-EXCALIBER scripts
• Running the RV-EXCALIBER scripts
• Expected outputs from the RV-EXCALIBER scripts
• A note on result interpretation
• Contact information
• Citation

Welcome to the official GitHub repository for Rare Variant EXome CALIBration using External Repositories (RV-EXCALIBER), a comprehensive framework for conducting rare variant association testing by leveraging large exome sequencing repositories as controls.

The summary-level allele frequencies derived from publicly available exome sequencing databases, such as the genome Aggregation Database (gnomAD), represent those expected to be found the general population and can therefore be leveraged as a control distribution. However, using such databases in their raw form will lead to biased genetic associations.

RV-EXCALIBER works by instituting 1) a correction that accounts for the increased variance due to the presence of rare variants in linkage disequilibrium, 2) an individual-level correction factor (iCF) that accounts for global variations in population substructure and sequencing technology, and 3) a gene-level correction factor (gCF) that accounts for granular deviations in allele burden bias amongst gene sets.

gnomAD: an external whole-exome sequencing repository containing summary-level allele frequency information that will be used as a control dataset in RV-EXCALIBER

internal testing dataset: a whole-exome sequencing dataset containing individual-level genotypes that will be used in conjunction with the summary-level allele frequencies from gnomAD to conduct a gene-based rare variant association test using RV-EXCALIBER

internal ranking dataset: a whole-exome sequencing dataset containing individual-level genotypes that will be used to calibrate association signals from the gene-based rare variant association test that was conducted using your testing dataset and gnomAD (see definition of testing dataset above)

RV-EXCALIBER can be downloaded locally by executing a git clone command as follows (please ensure you have git installed prior to issuing the following commands):

# Enter a directory in which RV-EXCALIBER is to be downloaded

cd genetics

# Initialize a local git repository

git init

# Execute the git clone command

git clone https://github.com/GMELab/RV-EXCALIBER.git

Expected download time: 1-5 minutes

Space required on local system: 155 MB

# Enter the cloned RV-EXCALIBER directory

cd /genetics/RV-EXCALIBER

# View directory architecture

ls -l | awk '{ print $9 }'
gnomAD_211_exomes_hg19
gnomAD_211_exomes_hg19_coverage
gnomAD_211_exomes_hg19_filter
default_ranking_dataset
scripts

Contains re-formatted per-variant annotation files (per autosomal chromosome) containing refGene, gnomAD 211, and dbNSFP annotations.

Contains a .bed file containing "high coverage coding" autosomal regions from gnomAD. A "high coverage coding" region is defined as an individual exon where a mean read depth of 20X was achieved for at least 90% of individuals in gnomAD. This .bed file will be intersected against variant sites in your internal testing and/or ranking dataset and in gnomAD if the coverage input variable is set to hcc in the get_RVBurdenMatrix_internal_rvexcaliber.clean.sh and the get_RVBurdenMatrix_gnomAD_rvexcaliber.clean.sh scripts. Please see a more detailed explanation of input variables in the Definitions of script inputs section.

Contains a .txt file with the FILTER status of all gnomAD autosomal variant sites. The per-variant gnomAD FILTER status is used to eliminate variants that are mutually observed in both your internal case dataset and gnomAD which do not have a FILTER status of PASS in gnomAD.

Contains a .txt file with the summary associations obtained by using RV-EXCALIBER on n = 6,082 healthy controls from the Myocardial Infarction Genetics (MIGen) consortium. This will act as your internal ranking dataset if the internal_ranking_dataset input variable is set to NA in the get_SummaryAssociations_iCF_gCF_rvexcaliber.clean.sh script. Please see a more detailed explanation of input variables in the Definitions of script inputs section of the README.

Contains the RV-EXCALIBER shell and Rscripts.

The quality-control of genetic data is essential towards achieving robust downstream results. A thorough description of the variant-level and sample-level quality control procedures that we implemented on our exome sequencing datasets can be found in the Supplementary Appendix (Section 4-6) from our bioRxiv preprint: https://www.biorxiv.org/content/10.1101/2020.02.03.931519v1

The current release of RV-EXCALIBER has been designed to conduct rare variant association tests using whole-exome sequencing (WES) data. The only input required are the binary plink files (i.e. .bed, .bim, .fam) from your WES dataset(s) of choice. If your WES data is stored in Variant Call File (VCF) format, it can easily be converted to plink binary format using the --make-bed command within plink v1.9 as shown below:

plink=/genetics/tools/plink_v1.9/plink
plink --noweb \
      --vcf /genetics/inData/MIGen_BioImage_ExS_QC.bg.vcf.gz \
      --vcf-half-call missing \
      --keep-allele-order \
      --make-bed \
      --out /genetics/inData/MIGen_BioImage_ExS_QC

If your VCF is correctly formatted, it will be converted to plink binary format:

cd /genetics/inData
ls -l | awk '{ print $9 }'
MIGen_BioImage_ExS_QC.bg.vcf.gz
MIGen_BioImage_ExS_QC.bg.vcf.gz.tbi
MIGen_BioImage_ExS_QC.bed
MIGen_BioImage_ExS_QC.bim
MIGen_BioImage_ExS_QC.fam
MIGen_BioImage_ExS_QC.log
MIGen_BioImage_ExS_QC.nosex

RV-EXCALIBER makes use of 3 dependencies to prepare the rare variant burden matrices that are used as input for rare variant association testing. Please ensure to have these dependencies downloaded to directories with appropriate read/write privileges.

The more updated version of the following R packages are still compatible. Please select the appropriate CRAN mirror during the install.

RV-EXCALIBER requires users to execute 4 shell scripts which function to 0) annotate variants in the internal testing and/or ranking dataset with the refGene, gnomAD 211, and dbNSFP databases, 1) generate a rare variant burden matrix for your internal testing and/or ranking dataset, 2) generate a rare variant burden matrix for gnomAD, and 3) generate summary gene-based rare variant association statistics

This section functions to clearly define the necessary inputs for each RV-EXCALIBER shell script. For examples on how to run the following scripts, please refer to the Running the RV-EXCALIBER scripts section

Purpose: Download the refGene, gnomAD 211, and dbNSFP databases to ANNOVAR's /humandb directory, if not done so already

Purpose: Generate a rare variant burden matrix for the internal testing and/or ranking datasets Number of required input variables: 9

Purpose: If your internal testing and/or ranking rare variant burden matrices generated in the get_RVBurdenMatrix_internal_rvexcaliber.clean.sh script were split by chromosome, this script combines the resulting split rare variant burden matrices into a single matrix Number of required input variables: 7

Purpose: Generate a rare variant burden matrix for gnomAD that corresponds to you internal testing and/or ranking dataset Number of required input variables: 6

Purpose: Generate iCF and gCF-adjusted gene-based rare variant association statistics Number of required input variables: 8

Each of the 4 shell scripts that are to be executed by the user will require either the path to the unpacked rvexcaliber directory or the paths to the dependencies to be defined. Variable names to these paths will be present in the header information of each script

1 path to be defined by the user

Important: the above script executes ANNOVAR's annotate_variation.pl script to download the necessary annotation databases used by RV-EXCALIBER (refGene, gnomAD 211, and dbNSFP)

4 paths to be defined by the user

1 path to be defined by the user

1 path to be defined by the user

1 path to be defined by the user

Before running the RV-EXCALIBER scripts, please ensure each script has been prepared correctly by referring the Preparing the RV-EXCALIBER scripts section

# Define path to shell scripts
rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

cd ${rvexcaliber_shell_scripts}

bash get_ANNOVAR_annotations.clean.sh

9 required inputs (see Definitions of script inputs for their definitions)

# Define path to shell scripts

rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

# If your internal testing or ranking dataset ** IS NOT ** split by chromosome  (or split in any other fashion):

cd ${rvexcaliber_shell_scripts}

bash get_RVBurdenMatrix_internal_rvexcaliber.clean.sh /genetics/inData /genetics/outData MIGen_BioImage_ExS_QC 0.001,0.005 0.01 0.009,0.025 nfe Y 10

# Define path to shell scripts
rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

# If your internal testing or ranking dataset ** IS ** split by chromosome (or split in some other fashion):

cd ${rvexcaliber_shell_scripts}

for chr in {1..22}; do

    bash get_RVBurdenMatrix_internal_rvexcaliber.clean.sh /genetics/inData /genetics/outData ${chr}_MIGen_BioImage_ExS_QC 0.001,0.005 0.01 0.009,0.025 nfe hcc 10

done

7 required inputs (see Definitions of script inputs for their definitions)

Important: If your internal testing or ranking datasets is split by chromosome, or split in any other way, (Condition B above), you must combine the split rare variant burden matrices (that were generated in the get_RVBurdenMatrix_internal_rvexcaliber.clean.sh script) into a single matrix. This can be done with the get_combined_internal_matrix_rvexcaliber.clean.sh script as follows:

# Step 1: Generate concatenated list of the full paths to the rare variant burden matrix for each chromosome that was generated in the 'get_RVBurdenMatrix_internal_rvexcaliber.clean.sh' script

for chr in {1..22}; do

  for gnomAD_MAF_threshold in 0.001 0.005; do

    for MCAP_threshold in 0.009 0.025; do

      for eth in nfe; do

        for coverage in hcc; do

          echo /genetics/outData/${chr}_MIGen_BioImage_ExS_QC_RVBurdenMatrix_${gnomAD_MAF_threshold}_${MCAP_threshold}_${eth}_${coverage}.txt.gz >> /genetics/outData/MIGen_BioImage_ExS_QC_RVBurdenMatrix_filelist.txt

       done

      done

    done

  done

done


# Step 2: Execute the get_combined_internal_matrix_rvexcaliber.clean.sh script as follows:

bash get_combined_internal_matrix_rvexcaliber.clean.sh /genetics/outData MIGen_BioImage_ExS_QC 0.001,0.005 0.009,0.025 nfe hcc /genetics/outData/MIGen_BioImage_ExS_QC_RVBurdenMatrix_filelist.txt

6 required inputs (see Definitions of script inputs for their definitions)

# Define path to shell scripts:

rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

# If your internal testing or ranking ** IS NOT ** split by chromosome:

cd ${rvexcaliber_shell_scripts}

bash get_RVBurdenMatrix_gnomAD_rvexcaliber.clean.sh /genetics/outData MIGen_BioImage_ExS_QC 0.001,0.005 0.009,0.025 nfe hcc

8 required inputs (see the Definitions of script inputs section for their definitions)

# Define path to shell scripts:

rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

# If you ** DO NOT** have a dedicated ranking dataset, set the 'internal_ranking_dataset' variable to NA

cd ${rvexcaliber_shell_scripts}

bash get_SummaryAssociations_iCF_gCF_rvexcaliber.clean.sh /genetics/outData MIGen_BioImage_ExS_QC NA 0.001,0.005 0.009,0.025 nfe hcc fulladjust

# Define path to shell scripts:

rvexcaliber_shell_scripts=/genetics/RV-EXCALIBER/scripts/shell

# If you ** DO ** have a dedicated ranking dataset, set the 'internal_ranking_dataset' variable to the name of your internal ranking dataset, which you defined in the 'get_RVBurdenMatrix_gnomAD_rvexcaliber.clean.sh' script

cd ${rvexcaliber_shell_scripts}

bash get_SummaryAssociations_iCF_gCF_rvexcaliber.clean.sh /genetics/outData MIGen_BioImage_ExS_QC MIGen_OHS_ExS_QC 0.001,0.005 0.009,0.025 nfe hcc fulladjust

For the purposes of explaining the output files, we will assign the following naming and filter variables to:

internal_dataset = MIGen_BioImage_ExS_QC

gnomAD_MAF_threshold = 0.001

MCAP_threshold = 0.025

eth = nfe

coverage = hcc

adjust_type = fulladjust

For the definition of these variables (and other variables), please refer to the Definitions of script inputs section.

This script will generate 4 output files

This script will generate 1 output file

This script will generate 2 output files

This script will generate 10 output files

The genomic inflation factor (GIF) measures the calibration of the distribution of the observed P-values to a null uniform distribution of P-values. While GIF values closer to 1 are desired, the extent of the calibration is largely a function of sample size, variant filtering criteria, and the comparator gnomAD ethnicity.

Therefore, if your internal testing dataset has a small sample size (i.e. n < 500), is filtered against stringent gnomAD MAF and/or MCAP thresholds, or is assessed against a more genetically diverse ethnic group in gnomAD, then the RV-EXCALIBER rare variant association test will be inherently underpowered, which will result in deflation of the GIF.

Any queries pertaining to the RV-EXCALIBER scripts or methodological framework can be addressed to either: Ricky Lali ([email protected]) or Guillaume Paré ([email protected])

Lali, R., Chong, M., Omidi, A. et al. Calibrated rare variant genetic risk scores for complex disease prediction using large exome sequence repositories. Nat Commun 12, 5852 (2021). https://doi.org/10.1038/s41467-021-26114-0