Genomics and Machine Learning Uncover Diapause Mechanisms and Predict Voltinism in the Bark Beetle Ips typographus

================================================================================
Genomics and Machine Learning Uncover Diapause Mechanisms and Predict 
Voltinism in the Bark Beetle Ips typographus
================================================================================

DATASET DESCRIPTION
-------------------
This repository contains the code and analytical pipeline used to identify 
genetic polymorphisms associated with diapause phenotypes in the European 
spruce bark beetle, Ips typographus. The analysis combines genome-wide 
association study (GWAS) with machine learning (ML) to predict facultative 
versus obligate diapause phenotypes, with implications for understanding 
voltinism and outbreak risk.

================================================================================
STUDY OVERVIEW
================================================================================

Citation: 
Palmieri, L., Dowle, E.J., Nadachowska-Brzyska, K., Schopf, A., Dobart, N., 
Stauffer, C., Schuler, H., Ragland, G.J., & Schebeck, M. (2025). Genomics 
and Machine Learning Uncover Diapause Mechanisms and Predict Voltinism in 
the Bark Beetle Ips typographus. Molecular Ecology.

Authors: 
Luciano Palmieri

Contact: luciano.palmierirocha@unibz.it

code revised with Claude Sonnet 4.5

================================================================================
ABBREVIATIONS AND CODES
================================================================================

Sample Population Codes
------------------------
* LOW: Central European low-altitude population (Prinzersdorf, Austria; 
  48°22'N, 15°48'E; 350 m)
* HIGH: Central European high-altitude population (Gesäuse, Austria; 
  47°35'N, 14°38'E; 1500 m)
* NORTH: Northern European population (Vindeln, Sweden; 64°12'N, 19°43'E; 
  300 m)

Phenotype Classifications
-------------------------
* Facultative diapause: Individuals capable of producing multiple generations 
  per year under favorable conditions
* Obligate diapause: Individuals with univoltine life cycle regardless of 
  environmental conditions

Technical Abbreviations
-----------------------
* GWAS: Genome-Wide Association Study
* ML: Machine Learning
* BF: Bayes Factor (strength of association in GWAS)
* VCA: Variance Component Analysis
* VC: Variance Component
* SNP: Single Nucleotide Polymorphism
* ddRAD: Double-Digest Restriction-Associated DNA sequencing
* AUROC: Area Under the Receiver Operating Characteristic Curve
* GO: Gene Ontology
* JHE: Juvenile Hormone Esterase
* LG: Linkage Group (chromosome)

================================================================================
FILE DESCRIPTIONS
================================================================================

1. params-Ips_ipyrad.txt
-------------------------
Description: Parameter file for ipyrad v.0.9.95 assembly
Purpose: Configuration for ddRAD sequence assembly and SNP calling
Usage: Input for ipyrad pipeline to generate VCF files from raw sequencing data

2. phenotype_pruned_filtered75missing_LD.vcf
---------------------------------------------
Description: Filtered and LD-pruned VCF file for phenotyped individuals
Purpose: Final quality-controlled SNP dataset for GWAS and ML analyses
Content: 11,867 SNPs across 296 phenotyped individuals after filtering for:
  - Maximum 25% missing data per SNP
  - Linkage disequilibrium (r² < 0.2)
Populations included: LOW (obligate and facultative), HIGH (obligate and 
facultative), NORTH (obligate)
Format: Standard VCF format (can be opened with vcftools, PLINK, or R packages)

3. WILD.recode.vcf
------------------
Description: Filtered VCF file for wild (unphenotyped) individuals
Purpose: SNP dataset for phenotype prediction in natural populations
Content: Same 11,867 SNPs as phenotyped dataset, across 232 wild individuals
Populations included: LOW (n=79), HIGH (n=77), NORTH (n=76)
Format: Standard VCF format

4. pop_assign.txt
-----------------
Description: Population assignment file
Purpose: Maps sample IDs to population groups for Baypass analysis
Format: Two-column tab-separated file with no header:
  - Column 1: Sample ID (matching VCF sample names)
  - Column 2: Population label (either "Diapause" or "NoDiapause" for 
    phenotyped; "LOW", "HIGH", or "NORTH" for wild)
Usage: Required input for converting VCF to Baypass format

5. reshaper_baypass.py
----------------------
Description: Python script for VCF to Baypass format conversion
Purpose: Converts standard VCF files to Baypass genotype format (allele 
counts per population)
Usage: python reshaper_baypass.py input.vcf pop_assign.txt output.geno
Input: VCF file and population assignment file
Output: Baypass-compatible genotype matrix (.geno file)
Dependencies: Python 3.x, standard libraries
Source: https://gitlab.com/YDorant/Toolbox/-/blob/master/reshaper_baypass.py

6. VCA_and_PopGEN_script.r
---------------------------
Description: R script for population structure analysis
Purpose: Performs Variance Component Analysis (VCA) and STRUCTURE analysis
Input: Filtered VCF file from ipyrad
Output: 
  - VC scores for population structure correction
  - VC plots showing geographic clustering
  - STRUCTURE plots showing admixture proportions
Dependencies: vcfR, adegenet, rrBLUP, STRUCTURE software

7. Filtering_Liftoff_and_Baypass_GWAS.txt and 
   Filtering_Liftoff_and_Baypass_GWAS.bat
-----------------------------------------------
Description: Complete pipeline for VCF quality control, gene annotation, 
and GWAS analysis
Purpose: 
  - Filter SNPs for missing data and linkage disequilibrium
  - Map gene annotations using Liftoff
  - Run Baypass GWAS analysis
  - Calculate median Bayes Factors across runs
Input: Raw VCF file, reference genome, gene annotations
Output:
  - Filtered VCF files
  - Gene annotation mappings
  - Baypass association results with Bayes Factors
  - List of significant SNPs
Dependencies: vcftools, plink, Liftoff, Baypass, R

8. Hyperparameter_tuning.py
----------------------------
Description: Python script for LightGBM hyperparameter optimization
Purpose: Uses Optuna to find optimal model parameters via Bayesian optimization
Input: Filtered SNP genotype matrix with phenotype labels
Output: Optimized hyperparameters for LightGBM classifier
Dependencies: Python 3.x, optuna, lightgbm, pandas, numpy, scikit-learn

9. LightGBM_training_KFold.py
------------------------------
Description: Python script for training and evaluating LightGBM classifier
Purpose: 
  - Train gradient boosting model on phenotyped individuals
  - Perform 50-fold cross-validation
  - Rank SNPs by feature importance
  - Calculate AUROC and accuracy metrics
Input: 
  - SNP genotype matrix for phenotyped individuals
  - Top variance components (VC1, VC2, VC3)
  - Optimized hyperparameters from Hyperparameter_tuning.py
Output:
  - Trained model file (.pkl)
  - SNP importance rankings
  - Cross-validation performance metrics
Dependencies: Python 3.x, lightgbm, pandas, numpy, scikit-learn, joblib

10. Wild_phenotype_classifier.py
---------------------------------
Description: Python script for predicting diapause phenotypes in wild 
populations
Purpose: Apply trained LightGBM model to classify wild beetles
Input: 
  - SNP genotype matrix for wild individuals
  - Trained LightGBM model (.pkl)
  - VC scores for wild individuals
Output: 
  - Predicted phenotypes for each individual
  - Confidence scores
  - Population-level phenotype frequency estimates
Dependencies: Python 3.x, lightgbm, pandas, numpy, joblib

11. SNP_mapping_and_GO_enrichment_script.r
-------------------------------------------
Description: R script for functional annotation and GO enrichment analysis
Purpose: 
  - Map significant SNPs to genomic features
  - Identify genes harboring significant SNPs
  - Test for GO term enrichment
Input: 
  - Significant SNPs from GWAS (BF > 6) and ML analysis
  - Gene annotation file (GFF/GTF)
  - Gene-to-GO mapping file
  - Gene length file
Output:
  - SNP-to-gene mappings
  - GO enrichment results
  - Enriched biological pathways
Dependencies: GenomicRanges, goseq

================================================================================
ANALYTICAL PIPELINE (ORDER OF EXECUTION)
================================================================================

STEP 1: RAD Sequencing and SNP Calling
---------------------------------------
Command: ipyrad -p params-Ips_ipyrad.txt -s 1234567
Input: Raw ddRAD sequencing reads (available at NCBI SRA: PRJDB38037)
Output: VCF file with ~134,410 SNPs across 574 individuals

STEP 2: Population Structure Analysis
--------------------------------------
Command: Rscript VCA_and_PopGEN_script.r
Input: Raw VCF file from Step 1
Output: VC scores, population structure plots

STEP 3: VCF Filtering, Gene Annotation, and GWAS
-------------------------------------------------
Command (Windows): Filtering_Liftoff_and_Baypass_GWAS.bat
Command (Linux/Mac): bash Filtering_Liftoff_and_Baypass_GWAS.txt

Input: 
  - Filtered VCF from ipyrad
  - Reference genome and annotations
  - VC scores from Step 2

Output: 
  - Filtered VCF (11,867 SNPs after QC and LD pruning)
  - Gene annotations
  - Baypass results with Bayes Factors for each SNP

Key GWAS Results:
  - 185 SNPs with BF > 2 (weak association)
  - 48 SNPs with BF > 6 (moderate association)
  - 13 SNPs with BF > 10 (strong association)
  - Top SNP: Juvenile hormone esterase (JHE) locus (BF = 18.42)

STEP 4: Machine Learning - Hyperparameter Optimization
-------------------------------------------------------
Command: python Hyperparameter_tuning.py
Input: Phenotyped individuals SNP matrix + VC scores
Output: Optimal hyperparameters (learning rate, number of leaves, depth, etc.)

STEP 5: Machine Learning - Model Training and Validation
---------------------------------------------------------
Command: python LightGBM_training_KFold.py

Input: 
  - Phenotyped individuals (296 samples)
  - Top 25 ML-ranked SNPs + GWAS SNPs (BF > 6)
  - VC1, VC2, VC3 scores
  - Optimized hyperparameters

Output: 
  - Trained model: BEST_model_phenotyped_kfold_lgbmAUROC.pkl
  - SNP importance rankings
  - Performance: AUROC = 0.918, Mean Accuracy = 0.856

Top ML Results:
  - JHE locus ranked #1 (importance score = 2847.81)
  - Only 2 SNPs overlapped between GWAS BF>10 and ML top 25
  - Best model: ML top 25 + GWAS BF>6 + VCA

STEP 6: Phenotype Prediction in Wild Populations
-------------------------------------------------
Command: python Wild_phenotype_classifier.py

Input: 
  - Wild individuals SNP matrix (232 samples)
  - Trained model from Step 5
  - VC scores for wild individuals

Output: Predicted phenotype frequencies per population

Wild Population Predictions:
  - NORTH (n=76): 98% obligate diapause (univoltine)
  - HIGH (n=77): 46% obligate, 54% facultative
  - LOW (n=79): 41% obligate, 59% facultative (multivoltine potential)

STEP 7: Functional Annotation and GO Enrichment
------------------------------------------------
Command: Rscript SNP_mapping_and_GO_enrichment_script.r

Input: 
  - Significant SNPs from GWAS and ML
  - Gene annotations
  - Gene lengths
  - Gene-to-GO mappings

Output: 
  - SNP-gene overlaps
  - Enriched GO terms
  - Biological pathways implicated in diapause

Key Enriched Pathways (both GWAS and ML):
  - Chromosome organization (GO:0051276)
  - Endocytosis (GO:0006897)
  - Scavenger receptor activity (GO:0005044)

================================================================================
REQUIRED SOFTWARE AND VERSIONS
================================================================================

Bioinformatics Tools
--------------------
* ipyrad: v0.9.95 (RAD-seq assembly)
* vcftools: v0.1.16 (VCF filtering)
* plink: v1.90b6.21 (LD pruning)
* Liftoff: v1.6.3 (gene annotation transfer)
* Baypass: v2.41 (GWAS analysis)
* STRUCTURE: v2.3 (population structure)

R and R Packages
----------------
* R: v4.4.2
* vcfR: v1.7.0.91
* adegenet: v1.3.1
* rrBLUP: v4.6.3
* GenomicRanges (Bioconductor)
* goseq: v1.60.0
* pophelper: v2.3.1

Python and Python Packages
---------------------------
* Python: v3.8+
* lightgbm: v4.5.0
* optuna: v4.2.0
* pandas: v2.1.4
* numpy: v2.1.3
* scikit-learn: v1.5.2
* joblib: v1.4.2

================================================================================
DATA FILES REQUIRED
================================================================================

Input Data Files (NOT INCLUDED)
----------------
These files should be obtained separately or generated from raw sequencing data:

1. Raw sequencing reads: Available at NCBI SRA (BioProject: PRJDB38037)
2. Reference genome: Ips typographus genome assembly 
   (Nadachowska-Brzyska et al., 2025)
3. Gene annotations: Transferred from Powell et al. (2021) using Liftoff
4. Phenotype data: Binary labels (0/1) for obligate/facultative diapause 
   (determined experimentally)
5. Gene-to-GO mapping: Gene IDs with associated GO terms
6. Gene lengths: Gene IDs with exon-based lengths in base pairs

Data Files Included in This Repository
---------------------------------------
* phenotype_pruned_filtered75missing_LD.vcf: Filtered SNPs for 296 phenotyped 
  individuals
* WILD.recode.vcf: Filtered SNPs for 232 wild individuals
* pop_assign.txt: Population assignments for all samples

================================================================================
KEY PARAMETERS AND THRESHOLDS
================================================================================

VCF Filtering
-------------
* Missing data threshold: 25%
* LD pruning: r² = 0.2, window = 50 SNPs, step = 10 SNPs
* Final SNP count: 11,867 (after filtering from 134,410)

Baypass GWAS
------------
* MCMC runs: 5 independent runs with different seeds
* Burn-in: 2,500 iterations
* Pilot runs: 20 (npilot)
* Pilot length: 200 iterations
* Association strength interpretation:
  - BF > 2: Weak association
  - BF > 6: Moderate association
  - BF > 10: Strong association

LightGBM Machine Learning
--------------------------
* Cross-validation: 50-fold
* Hyperparameter optimization: 2,000 Optuna trials
* Binary classification threshold: 0.5
* Missing data handling: Native LightGBM imputation
* Feature encoding: 0 (homozygous ref), 2 (heterozygous), 3 (homozygous alt)

GO Enrichment
-------------
* Background: All genes in RAD-seq dataset
* Foreground: Genes with SNPs (BF > 6 for GWAS; non-zero importance for ML)
* Test: Wallenius' noncentral hypergeometric distribution
* Length bias correction: Included via goseq

================================================================================
EXPECTED OUTPUTS AND INTERPRETATION
================================================================================

GWAS Results
------------
* High BF SNPs: Strong candidates for diapause regulation
* JHE locus: Top-ranked SNP suggests juvenile hormone pathway involvement
* Moderate BF SNPs: Additional candidates for multi-locus architecture

Machine Learning Results
------------------------
* High importance SNPs: Most informative for phenotype prediction
* AUROC > 0.9: Excellent discriminative ability
* Accuracy ~86%: Robust phenotype classification

Wild Population Predictions
----------------------------
* High proportion of obligate diapause (>90%): Univoltine populations, lower 
  outbreak risk
* High proportion of facultative diapause (>50%): Multivoltine potential, 
  higher outbreak risk
* Intermediate proportions: Variable voltinism depending on environmental 
  conditions

================================================================================
CITATIONS
================================================================================

Key Methods References
----------------------
* ipyrad: Eaton & Overcast (2020) Bioinformatics 36:2592-2594
* vcftools: Danecek et al. (2011) Bioinformatics 27:2156-2158
* plink: Purcell et al. (2007) Am J Hum Genet 81:559-575
* Liftoff: Shumate & Salzberg (2021) Bioinformatics 37:1639-1643
* Baypass: Gautier (2015) Genetics 201:1555-1579
* LightGBM: Ke et al. (2017) Advances in Neural Information Processing 
  Systems 30:3149-3157
* goseq: Young et al. (2010) Genome Biology 11:R14

Reference Genome
----------------
* Powell et al. (2021) Communications Biology 4:1059
* Nadachowska-Brzyska et al. (2025) in preparation

================================================================================
HOW TO CITE THIS DATASET
================================================================================

Recommended citation format:

Palmieri, L., Dowle, E.J., Nadachowska-Brzyska, K., Schopf, A., Dobart, N., 
Stauffer, C., Schuler, H., Ragland, G.J., & Schebeck, M. (2025). Genomics 
and Machine Learning Uncover Diapause Mechanisms and Predict Voltinism in 
Ips typographus - Analysis Pipeline and Data [Data set]. Zenodo. 
https://doi.org/10.5281/zenodo.17543185

Last updated: November 2025