Tutorial

1. Data input #

Example files are provided in the example directory: test.bed, test.bim, test.fam, and test_simu_pheno.txt.

Note: The example dataset is intentionally small and should only be used for testing program functionality. To obtain stable and reliable results, analyses should be performed on large-scale datasets.

PLINK Binary File #

PLINK binary files consist of three components:

• *.bed (binary genotype file): Stores genotype data in a compact binary format.
• *.bim (variant information file): Contains SNP identifiers, chromosome positions, and allele information.
• *.fam (sample information file): Includes individual IDs, family structure, sex, and phenotype data.

In fastGxE, missing genotype values for the focus SNP are imputed using the mean genotype. However, it is recommended to preprocess missing genotypes with Beagle or other imputation software for improved accuracy.

Data File Format #

The data file should contain individual IDs, confounding covariates, interacting environmental covariates, and phenotypic values. It must follow the structure below:

1. Header Line: A header row is required.
2. Delimiter: Data should be separated by blanks or tabs.
3. First Column: The first column must be individual IDs.
4. Matching Individuals:
   • Only individuals present in both the PLINK binary file and the data file will be used in the analysis.
   • Individuals missing in either the PLINK binary file or the data file will be automatically removed by fastGxE during the analysis.

An example data file

Associated options in fastGxE

NA Na na NAN NaN nan -NAN -NaN -nan N/A n/a

Any symbol listed above will be treated as a missing value. Symbols are separated by spaces, and users can add custom missing values to the list. If a variable is included in the model and contains a missing value in this field, the entire row will be discarded.

2. Calculate the Genomic Relationship Matrix (GRM) #

Single-Step Calculation #

For standard datasets, compute the GRM in a single step:

fastgxe --make-grm --bfile test --code-type 2 --out test

Parallel Computation for Large Datasets (e.g., Biobank Data) #

For large-scale datasets, partition the GRM into 100 parts, compute them in parallel, and then merge the results:

fastgxe --make-grm --bfile test --code-type 2 --npart 100 1 --out test
fastgxe --make-grm --bfile test --code-type 2 --npart 100 2 --out test
...
fastgxe --make-grm --bfile test --code-type 2 --npart 100 100 --out test

Once all parts are computed, merge them into a single GRM file:

fastgxe --process-grm --merge --grm test.agrm --npart 100

Sample Clustering Based on Relatedness #

To group samples such that any pair with an estimated relatedness greater than 0.05 belongs to the same subgroup:

fastgxe --process-grm --group --grm test.agrm --cut-value 0.05 --out test.agrm

3. Genome-Wide GxE Analysis #

Using the previously generated sparse GRM, along with a PLINK binary file and a data file, perform a GxE analysis as follows:

fastgxe --test-gxe --grm test --bfile test --data test_simu_pheno.txt --trait pheno \
--env-int E1:E40 --out test_gxe

• The column pheno represents the trait values.
• Columns from E1 to E40 are considered interacting environmental factors.

Parallel Computation for Large Datasets #

For large datasets, partition the SNPs into 10 parts and compute them in parallel. The results can then be merged using a custom Python or R script.

fastgxe --test-gxe --split-task 10 1 --grm test --bfile test --data test_simu_pheno.txt --trait pheno \
--env-int E1:E40 --out test_gxe 

fastgxe --test-gxe --split-task 10 2 --grm test --bfile test --data test_simu_pheno.txt --trait pheno \
--env-int E1:E40 --out test_gxe 

...

fastgxe --test-gxe --split-task 10 10 --grm test --bfile test --data test_simu_pheno.txt --trait pheno \
--env-int E1:E40 --out test_gxe

Once all tasks are completed, the results can be merged using a custom Python or R script.

4. GxE Heritability Estimation Using the Method of Moments #

To obtain an unbiased estimate of GxE heritability, use the method of moments:

Note: we need large sample size to get stable results. The example data have too small size.

fastgxe --mom --bfile test --data test_simu_pheno.txt --trait pheno \
--env-int E1:E40 --out test_mom

5. Identifying Environmental Factors Driving GxE #

Once GxE genomic loci are identified, the environmental factors driving the GxE effect of the lead SNP can be determined using mmSuSiE. mmSuSiE is written in python can be found at https://github.com/chaoning/mmsusie

First reformat the genomic relationship matrix with fastgxe

# this step will produce the 
fastgxe --process-grm --reformat --grm test.agrm --out-fmt 1 --out test.agrm

from mmsusie import MMSuSiE
model = MMSuSiE()

# Define input files
bed_file = "test"              # PLINK binary genotype file prefix (.bed/.bim/.fam)
grm_file = "test"              # GRM (genetic relationship matrix) file prefix
pheno_file = "test_simu_pheno.txt"  # Phenotype file with environmental covariates and trait

# Define a list of interacting environmental variables: E1, E2, ..., E40
env_int = [f"E{i}" for i in range(1, 41)]

# Perform LD pruning to identify leading SNPs for GxE signals
assoc_file = "test_gxe.res"    # Association results from fastGxE or similar method
df_leading = model.ld_pure(assoc_file, bed_file, ld_r2=0.1, snp="SNP", p="p_gxe", p_cutoff=5e-8)
print(df_leading)              # Print the pruned set of lead SNPs

# Specify parameters for mmSuSiE fine-mapping
snp_id = "rs550011"            # Lead SNP to analyze
trait = "pheno"                # Target phenotype column name
varcom_file = "test_gxe.var"   # Variance components file estimated from fastGxE --test-gxe
out_file = "test_mmsusie"      # Output file prefix for mmSuSiE results

# Run the mmSuSiE model to identify environmental factors driving the GxE effect
res_dct = model.run(pheno_file, trait, env_int, grm_file, bed_file, snp_id, varcom_file, out_file)

# Print the credible sets (CS) identified by mmSuSiE (index start from 0)
print(res_dct['cs'])

# Extract posterior inclusion probabilities (PIPs)
df_pip = res_dct["pip"]
print(df_pip)

# Display environmental factors with PIP > 0.5
print(df_pip[df_pip["pip"] > 0.5])