DOI: 10.5281/zenodo.17543184
title: Genomics and Machine Learning Uncover Diapause Mechanisms and Predict Voltinism in the Bark Beetle Ips typographus

description: ================================================================================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 assemblyPurpose: Configuration for ddRAD sequence assembly and SNP callingUsage: 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 individualsPurpose: Final quality-controlled SNP dataset for GWAS and ML analysesContent: 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) individualsPurpose: SNP dataset for phenotype prediction in natural populationsContent: Same 11,867 SNPs as phenotyped dataset, across 232 wild individualsPopulations included: LOW (n=79), HIGH (n=77), NORTH (n=76)Format: Standard VCF format

4. pop_assign.txt-----------------Description: Population assignment filePurpose: Maps sample IDs to population groups for Baypass analysisFormat: 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 conversionPurpose: Converts standard VCF files to Baypass genotype format (allele counts per population)Usage: python reshaper_baypass.py input.vcf pop_assign.txt output.genoInput: VCF file and population assignment fileOutput: Baypass-compatible genotype matrix (.geno file)Dependencies: Python 3.x, standard librariesSource: https://gitlab.com/YDorant/Toolbox/-/blob/master/reshaper_baypass.py

6. VCA_and_PopGEN_script.r---------------------------Description: R script for population structure analysisPurpose: Performs Variance Component Analysis (VCA) and STRUCTURE analysisInput: Filtered VCF file from ipyradOutput:   - VC scores for population structure correction  - VC plots showing geographic clustering  - STRUCTURE plots showing admixture proportionsDependencies: 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 analysisPurpose:   - Filter SNPs for missing data and linkage disequilibrium  - Map gene annotations using Liftoff  - Run Baypass GWAS analysis  - Calculate median Bayes Factors across runsInput: Raw VCF file, reference genome, gene annotationsOutput:  - Filtered VCF files  - Gene annotation mappings  - Baypass association results with Bayes Factors  - List of significant SNPsDependencies: vcftools, plink, Liftoff, Baypass, R

8. Hyperparameter_tuning.py----------------------------Description: Python script for LightGBM hyperparameter optimizationPurpose: Uses Optuna to find optimal model parameters via Bayesian optimizationInput: Filtered SNP genotype matrix with phenotype labelsOutput: Optimized hyperparameters for LightGBM classifierDependencies: Python 3.x, optuna, lightgbm, pandas, numpy, scikit-learn

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

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

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

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

STEP 1: RAD Sequencing and SNP Calling---------------------------------------Command: ipyrad -p params-Ips_ipyrad.txt -s 1234567Input: 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.rInput: Raw VCF file from Step 1Output: VC scores, population structure plots

STEP 3: VCF Filtering, Gene Annotation, and GWAS-------------------------------------------------Command (Windows): Filtering_Liftoff_and_Baypass_GWAS.batCommand (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.pyInput: Phenotyped individuals SNP matrix + VC scoresOutput: 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 Liftoff4. Phenotype data: Binary labels (0/1) for obligate/facultative diapause    (determined experimentally)5. Gene-to-GO mapping: Gene IDs with associated GO terms6. 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