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: [e-mail] 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