Protocol for genetic analysis of population-scale ultra-low-depth sequencing data.

REAGENT or RESOURCE SOURCE IDENTIFIER Deposited data Example code (GitHub repository) This paper; published works 1 , 2 https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis ; https://doi.org/10.5281/zenodo.14263301 Human reference genome hg38 (GRCh38) UCSC Genome Browser https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ GATK resource bundle GATK gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf gs://gcp-public-data--broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz Correspondence file linking mutation sites to rsIDs dbSNP https://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/chr_rpts/ LD scores for heritability and genetic correlation analysis Bulik-Sullivan et al. 4 , 5 https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores Partitioned heritability analysis reference files Bulik-Sullivan et al. 4 , 5 https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays Example summary statistics from Japanese ENcyclopedia of GEnetic associations Biobank Japan 3 http://jenger.riken.jp/result GTEx V8 Cross Tissue Imputation (UTMOST) Prediction Models PredictDB 6 https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ 1000 Genome Project Phase 3 from the MAGMA website de Leeuw et al. 7 https://cncr.nl/research/magma/ Anatomical Therapeutic Chemical (ATC) classification of the drug data WHO Collaborating Centre for Drug Statistics Methodology https://atcddd.fhi.no The drug gene set built from DGIdb v.5 Cannon et al. 8 https://www.dgidb.org/data/latest/interactions.tsv Software and algorithms R (version 4.2.1 or current version) CRAN https://www.r-project.org/ ; RRID: SCR_003005 R Studio (2023.12.1 Build 402) Posit https://posit.co/products/open-source/rstudio/ ; RRID: SCR_000432 Python (3.6.0 or current version) Python Software Foundation https://www.python.org/downloads/ ; RRID: SCR_008394 Python (2.7.0) Python Software Foundation https://www.python.org/downloads/ ; RRID: SCR_008394 Burrows-Wheeler Alignment Tool (BWA, 0.7.1 or current version) Li et al. 9 https://github.com/lh3/bwa ; RRID: SCR_010910 GATK (4.0.0 or current version) GATK 10 https://github.com/broadinstitute/gatk/releases ; RRID: SCR_001876 SAMtools (1.15 or current version) Danecek et al. 11 http://www.htslib.org/ ; RRID: SCR_002105 BaseVarC (1.0) Zilong Li https://github.com/Zilong-Li/BaseVarC BCFtools (1.13 or current version) Danecek et al. 11 https://samtools.github.io/bcftools/howtos/install.html ; RRID: SCR_005227 PLINK2 (2.00 alpha) Chang et al. 12 https://www.cog-genomics.org/plink/2.0/ ; RRID: SCR_001757 Fastp Chen et al. 13 https://github.com/OpenGene/fastp ; RRID: SCR_016962 LDSC (1.0.1) Bulik-Sullivan et al. 4 , 5 https://github.com/bulik/ldsc ; RRID: SCR_022801 MAGMA (1.10) de Leeuw et al. 7 https://cncr.nl/research/magma/ ; RRID: SCR_005757 STITCH (1.5.3 or current version) Davies et al. 14 https://github.com/rwdavies/STITCH TwoSampleMR (0.6.8) Davey et al. 15 https://mrcieu.github.io/TwoSampleMR/ ; RRID: SCR_019010 Ieugwasr (1.0.1) Gibran et al. 16 https://github.com/MRCIEU/ieugwasr RSQLite (2.3.7 or current version) Müller et al. 17 https://cran.r-project.org/web/packages/RSQLite/index.html stringr (1.5.1 or current version) Hadley Wickham 18 https://cran.r-project.org/web/packages/stringr/index.html ggplot2 (3.4.4 or current version) Hadley Wickham 19 https://CRAN.R-project.org/package=ggplot2 ; RRID: SCR_014601

Conceptualization, J.Z. and H.Z.; methodology, J.Z., L.L., X.L., Y.L., X.Z., Y.W., and H.Z.; writing – original draft, J.Z.; writing – review and editing, J.Z., H.Z., and L.L.; visualization, J.Z.; supervision, H.Z., M.L., and X.J.

This protocol outlines a workflow for processing ultra-low depth (<1×) sequencing data on a population scale (over 10,000 samples). In theory, this pipeline can be applied to any type of ultra-low depth sequencing data, not limited to data from non-invasive prenatal testing (NIPT). Specifically, to generate ultra-low depth sequence data from NIPT, 5 mL peripheral blood was collected from pregnant women using vacuum tubes, with EDTA samples processed within 8 h at 2°C–8°C and ambient temperature tubes within 96 h at 6°C–35°C. Plasma was separated via two centrifugations (1,600 ×  g and 16,000 ×  g ). Free DNA was extracted, libraries prepared through end repair, adapter ligation, and PCR amplification. The libraries were then pooled, circularized, and converted into DNA nanoballs. Single-end sequencing was performed on a BGISEQ-500 platform with an average depth of 0.1× and a read length of 35 bp. For further details, please refer to the published work. 1 , 2 To closely examine ultra-low-depth sequencing data, we simulated NIPT data by down-sampling two female samples (HG00501 and HG00525) from the 1000 Genomes Project Phase 3. Specifically, single-end FASTQ files were downloaded from the project’s FTP site ( https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/ ) and down-sampled to 0.1× coverage using Seqtk ( https://github.com/lh3/seqtk ). The simulated NIPT sequence data has been made available on Google Cloud ( https://drive.google.com/drive/folders/1T6nGk9dZK7Fk_xNV1Ql8oQnt627oJO7d?usp=sharing ). Figure 1 The tree directory structure for the protocol The tree directory structure for the protocol Timing: 1–2 days 1. Initialize the Working Directory. a. Download scripts from the GitHub repository, the documentation is available at: https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis . b. Execute the code 0.Setup.sh: cd /path_to/Ultra-Low-Depth-Sequencing-Genetic-Analysis-main/code sh 0.Setup.sh CRITICAL: Throughout this article, the prefix /path_to denotes the working directory where the protocol is executed. Be sure to replace this placeholder with the path to your local working directory wherever it appears. Note: This step establishes the majority of the working directories utilized throughout the protocol and organizes all files into a structured tree directory. Readers are encouraged to compare their script output with Figure 1 to verify consistency and accuracy. 2. Install R, RStudio, Python 3.x, Python 2.x in a Linux, MacOS, or Windows environment. Note: We recommend running this protocol on a Linux-based operating system. Although the procedures may theoretically be compatible with other operating systems (e.g., Windows, macOS), they have not been tested or verified for use outside of Linux. The Table 1 lists the supported operating systems for each software based on their official installation instructions. For users on Windows 10 or later, the Windows Subsystem for Linux (WSL) provides full-featured Linux capabilities within the Windows environment. For more information on WSL installation and usage, please refer to the official documentation: https://learn.microsoft.com/en-us/windows/wsl/install . 3. Ensure sufficient high-performance computing (HPC) cluster resources are allocated for the analysis, and verify that Java Runtime Environment, R, and Python 3.x are installed on the cluster. Note: All the required software can be installed locally. However, we recommend using a HPC environment because GWAS analysis of ultra-low-depth sequence data typically involves a large number of samples, requiring significant computational resources. If a user's local device has sufficient storage capacity and computational power, they can install the software and run the protocol locally. a. Ensure that all software listed in the key resources table is installed on the HPC cluster. b. Ensure that all R packages listed in the key resources table are installed on the HPC cluster. Table 1 Compatibility of software used in this protocol across different operating systems Software & packages Linux Windows MacOS BWA 0.7.1 √ ○ ○ GATK 4.0.0 √ ○ √ SAMtools 1.15 √ ○ ○ BCFtools 1.13 √ ○ ○ BaseVarC 1.0 ○ ○ ○ PLINK 2 √ √ √ fastp √ ○ ○ MAGMA 1.10 √ √ √ python 3.6.0 √ √ √ python 2.7.0 √ √ √ Example code (this protocol) ○ ○ ○ LDSC 1.0.1 ○ ○ ○ R 4.2.1 √ √ √ R Studio √ √ √ STITCH 1.5.3 ○ ○ ○ TwoSampleMR 0.6.8 ○ ○ ○ ieugwasr 1.0.1 ○ ○ ○ RSQLite 2.3.7 ○ ○ ○ stringr 1.5.1 ○ ○ ○ ggplot2 3.4.4 ○ ○ ○ √, Completely compatible, provides binary executable files; ○, Principle compatibility, provides source code that can be compiled and run. Timing: 1–2 days Initialize the Working Directory. a. Download scripts from the GitHub repository, the documentation is available at: https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis . b. Execute the code 0.Setup.sh: Download scripts from the GitHub repository, the documentation is available at: https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis . Execute the code 0.Setup.sh: cd /path_to/Ultra-Low-Depth-Sequencing-Genetic-Analysis-main/code sh 0.Setup.sh CRITICAL: Throughout this article, the prefix /path_to denotes the working directory where the protocol is executed. Be sure to replace this placeholder with the path to your local working directory wherever it appears. Note: This step establishes the majority of the working directories utilized throughout the protocol and organizes all files into a structured tree directory. Readers are encouraged to compare their script output with Figure 1 to verify consistency and accuracy. Install R, RStudio, Python 3.x, Python 2.x in a Linux, MacOS, or Windows environment. Note: We recommend running this protocol on a Linux-based operating system. Although the procedures may theoretically be compatible with other operating systems (e.g., Windows, macOS), they have not been tested or verified for use outside of Linux. The Table 1 lists the supported operating systems for each software based on their official installation instructions. For users on Windows 10 or later, the Windows Subsystem for Linux (WSL) provides full-featured Linux capabilities within the Windows environment. For more information on WSL installation and usage, please refer to the official documentation: https://learn.microsoft.com/en-us/windows/wsl/install . Ensure sufficient high-performance computing (HPC) cluster resources are allocated for the analysis, and verify that Java Runtime Environment, R, and Python 3.x are installed on the cluster. Note: All the required software can be installed locally. However, we recommend using a HPC environment because GWAS analysis of ultra-low-depth sequence data typically involves a large number of samples, requiring significant computational resources. If a user's local device has sufficient storage capacity and computational power, they can install the software and run the protocol locally. a. Ensure that all software listed in the key resources table is installed on the HPC cluster. b. Ensure that all R packages listed in the key resources table are installed on the HPC cluster. Note: All the required software can be installed locally. However, we recommend using a HPC environment because GWAS analysis of ultra-low-depth sequence data typically involves a large number of samples, requiring significant computational resources. If a user's local device has sufficient storage capacity and computational power, they can install the software and run the protocol locally. Ensure that all software listed in the key resources table is installed on the HPC cluster. Ensure that all R packages listed in the key resources table are installed on the HPC cluster. Compatibility of software used in this protocol across different operating systems √, Completely compatible, provides binary executable files; ○, Principle compatibility, provides source code that can be compiled and run. Timing: 1–2 days Note: The required dependency files will occupy approximately 40 GB of space. Please ensure you have adequate storage prepared in advance. 4. Download the reference genome sequence for sequencing analysis. a. We recommend to use the “Hard-masked” assembly sequence of human reference genome hg38 (GRCh38). b. Download the ‘hg38.fa.masked.gz’ file, the size is 487MB. The file and documentation are available at: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ . c. Move the downloaded file to the directory ‘/code/Supporting_data/hg38’. CRITICAL: Before using the downloaded reference genome sequence, it must be indexed with software tools such as BWA and SAMtools; otherwise, the analysis cannot proceed. For detailed instructions, please refer to the BWA documentation ( https://bio-bwa.sourceforge.net/bwa.shtml ) and the SAMtools documentation ( https://www.htslib.org/doc/samtools-index.html ). Note: We recommend using the 'Hard-masked' sequence (hg38.fa.masked.gz from the download page) due to the typically low sequencing depth (<1×) and single-ended nature of the data. Repeat regions are not ideal for these datasets, as they can hinder analysis. The masked sequence improves alignment performance. For more details, visit: https://gatk.broadinstitute.org/hc/en-us/articles/17295731870235-Masked-reference-genomes . 5. Download the GATK resource bundle for quality control in bam files. The file sizes are listed on the download page. Note: The bundle files we used in this section refer to the official examples used by GATK ( https://github.com/gatk-workflows/gatk4-data-processing/blob/master/processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json ). Altogether, we need to download three files along with their index files. a. Download dbSNP release 138 at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf. Note: dbSNP VCF files are available in various release versions. The file provided here is based on the official parameters of GATK. We recommend using the latest version of dbSNP, which can be accessed here: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/ . b. Download Mills_and_1000G_gold_standard.indels.hg38.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz . c. Download Homo_sapiens_assembly38.known_indels.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz . d. Move all of the downloaded files to the directory ‘/code/Supporting_data/gatk-bundle’. 6. Download the correspondence file that links mutation sites to rsIDs for heritability and genetic correlation analysis. a. You can download the data at: https://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/chr_rpts/ , the data sizes are listed on the webpage. b. Move the downloaded file to the directory ‘/code/Supporting_data/dbSNP’. 7. Download the LD scores required for LDSC heritability and genetic correlation analysis. a. You can download the LD scores data at https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores . The data size for East Asian LD scores is 12.3 MB. b. Move the downloaded file to the directory ‘/code/Supporting_data/LDSC’. 8. Download the LDSC reference files required for Partitioned Heritability analysis. a. The required files are listed below and can be downloaded from the URL provided: https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays . The file sizes are listed on the webpage. b. Download the baseline model LD scores, built using 1KGP Phase3 EAS samples, from `1000G_Phase3_EAS_baselineLD_v2.2_ldscores.tgz`. c. Download the regression weights files, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_weights_hm3_no_MHC.tgz`. d. Download the allele frequencies, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_plinkfiles.tgz`. e. Download the cell-type LD score reference file, built from cell-types and tissues in the GTEx database and Franke Lab, including `Multi_tissue_gene_expr_EAS_1000G_v3_ldscores.tar.gz` and `Multi_tissue_gene_expr.EAS.ldcts`. f. Move all of the downloaded files to the directory ‘/code/Supporting_data/LDSC’. 9. Download the summary statistics for Mendelian Randomization study. Here we used summary data from the Japanese ENcyclopedia of GEnetic associations 3 as example. a. File sizes and documentation are available at: http://jenger.riken.jp/en/result . b. Refer to Table 2 for the download links to the summary data. Note: Please note that Table 2 in this section provides summary data for 46 BBJ (BioBank Japan) datasets, which encompass all the BBJ data used in our published studies. 1 , 2 Readers may choose to download a subset of the summaries based on their specific needs. Table 2 The download links for BBJ summary data BBJ-Trait-summary URL Rheumatoid arthritis in Asian http://jenger.riken.jp:8080/download/Rheumatoid_Arthritis_Asian Open-angle glaucoma http://jenger.riken.jp:8080/download/Open-angle_glaucoma_autosome Type 2 Diabetes http://jenger.riken.jp:8080/download/Type_2_Diabetes Smoking initiation http://jenger.riken.jp:8080/download/Smoking_initiation_autosomes_all Smoking cessation http://jenger.riken.jp:8080/download/Smoking_cessation_autosomes_all Arrhythmia http://jenger.riken.jp:8080/download/Arrhythmia Asthma http://jenger.riken.jp:8080/download/Asthma Atopic dermatitis http://jenger.riken.jp:8080/download/Atopic_dermatitis Biliary tract cancer http://jenger.riken.jp:8080/download/Biliary_tract_cancer Cataract http://jenger.riken.jp:8080/download/Cataract Cerebral aneurysm http://jenger.riken.jp:8080/download/Cerebral_aneurysm Cervical cancer http://jenger.riken.jp:8080/download/Cervical_cancer Chronic hepatitis B http://jenger.riken.jp:8080/download/Chronic_hepatitis_B Chronic hepatitis C http://jenger.riken.jp:8080/download/Chronic_hepatitis_C Chronic obstructive pulmonary disease http://jenger.riken.jp:8080/download/Chronic_obstructive_pulmonary_disease Cirrhosis http://jenger.riken.jp:8080/download/Cirrhosis Congestive heart failure http://jenger.riken.jp:8080/download/Congestive_heart_failure Drug eruption http://jenger.riken.jp:8080/download/Drug_eruption Endometrial cancer http://jenger.riken.jp:8080/download/Endometrial_cancer Endometriosis http://jenger.riken.jp:8080/download/Endometriosis Epilepsy http://jenger.riken.jp:8080/download/Epilepsy Esophageal cancer http://jenger.riken.jp:8080/download/Esophageal_cancer Gastric cancer http://jenger.riken.jp:8080/download/Gastric_cancer Glaucoma http://jenger.riken.jp:8080/download/Glaucoma Graves disease http://jenger.riken.jp:8080/download/Graves_disease Hematological malignancy http://jenger.riken.jp:8080/download/Hematological_malignancy Interstitial lung disease http://jenger.riken.jp:8080/download/Interstitial_lung_disease Ischemic stroke http://jenger.riken.jp:8080/download/Ischemic_stroke Keloid http://jenger.riken.jp:8080/download/Keloid Lung cancer http://jenger.riken.jp:8080/download/Lung_cancer Nephrotic syndrome http://jenger.riken.jp:8080/download/Nephrotic_syndrome Osteoporosis http://jenger.riken.jp:8080/download/Osteoporosis Ovarian cancer http://jenger.riken.jp:8080/download/Ovarian_cancer Pancreatic cancer http://jenger.riken.jp:8080/download/Pancreatic_cancer Periodontal disease http://jenger.riken.jp:8080/download/Periodontal_disease Peripheral artery disease http://jenger.riken.jp:8080/download/Peripheral_artery_disease Pollinosis http://jenger.riken.jp:8080/download/Pollinosis Prostate cancer http://jenger.riken.jp:8080/download/Prostate_cancer Pulmonary tuberculosis http://jenger.riken.jp:8080/download/Pulmonary_tuberculosis Rheumatoid arthritis http://jenger.riken.jp:8080/download/Rheumatoid_arthritis Type 2 diabetes http://jenger.riken.jp:8080/download/Type_2_diabetes Urolithiasis http://jenger.riken.jp:8080/download/Urolithiasis Uterine fibroids http://jenger.riken.jp:8080/download/Uterine_fibroids hepatocellular carcinoma http://jenger.riken.jp:8080/download/hepatocellular_carcinoma Breast cancer http://jenger.riken.jp:8080/download/Breast_cancer Coronary artery disease http://jenger.riken.jp:8080/download/Coronary_artery_disease 10. Download the expression files for transcriptome-wide association study (TWAS) study. a. Here, we used the GTEx V8 Cross-Tissue Imputation (UTMOST) Prediction Models from PredictDB as an example. b. Download and documentation are available at: https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . c. Move all of the downloaded files to the directory ‘/code/Supporting_data/GTEx’. 11. Download and order the data files for drug-target analysis. a. Order anatomical therapeutic chemical (ATC) classification index Note: We need to obtain the ATC classification index from the WHO Collaborating Centre for Drug Statistics Methodology, which is available at: https://atcddd.fhi.no/atc_ddd_index_and_guidelines/order/ . b. Download the drug-gene interaction data from The Drug Gene Interaction Database: https://www.dgidb.org/data/latest/interactions.tsv . c. Download the MAGMA software dependencies: https://cncr.nl/research/magma/ . d. Move all of the downloaded files to the directory ‘/code/Supporting_data/magma’. Timing: 1–2 days Note: The required dependency files will occupy approximately 40 GB of space. Please ensure you have adequate storage prepared in advance. Download the reference genome sequence for sequencing analysis. a. We recommend to use the “Hard-masked” assembly sequence of human reference genome hg38 (GRCh38). b. Download the ‘hg38.fa.masked.gz’ file, the size is 487MB. The file and documentation are available at: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ . c. Move the downloaded file to the directory ‘/code/Supporting_data/hg38’. We recommend to use the “Hard-masked” assembly sequence of human reference genome hg38 (GRCh38). Download the ‘hg38.fa.masked.gz’ file, the size is 487MB. The file and documentation are available at: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ . Move the downloaded file to the directory ‘/code/Supporting_data/hg38’. CRITICAL: Before using the downloaded reference genome sequence, it must be indexed with software tools such as BWA and SAMtools; otherwise, the analysis cannot proceed. For detailed instructions, please refer to the BWA documentation ( https://bio-bwa.sourceforge.net/bwa.shtml ) and the SAMtools documentation ( https://www.htslib.org/doc/samtools-index.html ). Note: We recommend using the 'Hard-masked' sequence (hg38.fa.masked.gz from the download page) due to the typically low sequencing depth (<1×) and single-ended nature of the data. Repeat regions are not ideal for these datasets, as they can hinder analysis. The masked sequence improves alignment performance. For more details, visit: https://gatk.broadinstitute.org/hc/en-us/articles/17295731870235-Masked-reference-genomes . Download the GATK resource bundle for quality control in bam files. The file sizes are listed on the download page. Note: The bundle files we used in this section refer to the official examples used by GATK ( https://github.com/gatk-workflows/gatk4-data-processing/blob/master/processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json ). Altogether, we need to download three files along with their index files. a. Download dbSNP release 138 at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf. Note: dbSNP VCF files are available in various release versions. The file provided here is based on the official parameters of GATK. We recommend using the latest version of dbSNP, which can be accessed here: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/ . b. Download Mills_and_1000G_gold_standard.indels.hg38.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz . c. Download Homo_sapiens_assembly38.known_indels.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz . d. Move all of the downloaded files to the directory ‘/code/Supporting_data/gatk-bundle’. Note: The bundle files we used in this section refer to the official examples used by GATK ( https://github.com/gatk-workflows/gatk4-data-processing/blob/master/processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json ). Altogether, we need to download three files along with their index files. Download dbSNP release 138 at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf. Note: dbSNP VCF files are available in various release versions. The file provided here is based on the official parameters of GATK. We recommend using the latest version of dbSNP, which can be accessed here: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/ . Note: dbSNP VCF files are available in various release versions. The file provided here is based on the official parameters of GATK. We recommend using the latest version of dbSNP, which can be accessed here: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/ . Download Mills_and_1000G_gold_standard.indels.hg38.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz . Download Homo_sapiens_assembly38.known_indels.vcf.gz at gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz . Move all of the downloaded files to the directory ‘/code/Supporting_data/gatk-bundle’. Download the correspondence file that links mutation sites to rsIDs for heritability and genetic correlation analysis. a. You can download the data at: https://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/chr_rpts/ , the data sizes are listed on the webpage. b. Move the downloaded file to the directory ‘/code/Supporting_data/dbSNP’. You can download the data at: https://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/chr_rpts/ , the data sizes are listed on the webpage. Move the downloaded file to the directory ‘/code/Supporting_data/dbSNP’. Download the LD scores required for LDSC heritability and genetic correlation analysis. a. You can download the LD scores data at https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores . The data size for East Asian LD scores is 12.3 MB. b. Move the downloaded file to the directory ‘/code/Supporting_data/LDSC’. You can download the LD scores data at https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores . The data size for East Asian LD scores is 12.3 MB. Move the downloaded file to the directory ‘/code/Supporting_data/LDSC’. Download the LDSC reference files required for Partitioned Heritability analysis. a. The required files are listed below and can be downloaded from the URL provided: https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays . The file sizes are listed on the webpage. b. Download the baseline model LD scores, built using 1KGP Phase3 EAS samples, from `1000G_Phase3_EAS_baselineLD_v2.2_ldscores.tgz`. c. Download the regression weights files, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_weights_hm3_no_MHC.tgz`. d. Download the allele frequencies, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_plinkfiles.tgz`. e. Download the cell-type LD score reference file, built from cell-types and tissues in the GTEx database and Franke Lab, including `Multi_tissue_gene_expr_EAS_1000G_v3_ldscores.tar.gz` and `Multi_tissue_gene_expr.EAS.ldcts`. f. Move all of the downloaded files to the directory ‘/code/Supporting_data/LDSC’. The required files are listed below and can be downloaded from the URL provided: https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays . The file sizes are listed on the webpage. Download the baseline model LD scores, built using 1KGP Phase3 EAS samples, from `1000G_Phase3_EAS_baselineLD_v2.2_ldscores.tgz`. Download the regression weights files, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_weights_hm3_no_MHC.tgz`. Download the allele frequencies, built using 1KGP Phase3 EAS samples, `1000G_Phase3_EAS_plinkfiles.tgz`. Download the cell-type LD score reference file, built from cell-types and tissues in the GTEx database and Franke Lab, including `Multi_tissue_gene_expr_EAS_1000G_v3_ldscores.tar.gz` and `Multi_tissue_gene_expr.EAS.ldcts`. Move all of the downloaded files to the directory ‘/code/Supporting_data/LDSC’. Download the summary statistics for Mendelian Randomization study. Here we used summary data from the Japanese ENcyclopedia of GEnetic associations 3 as example. a. File sizes and documentation are available at: http://jenger.riken.jp/en/result . b. Refer to Table 2 for the download links to the summary data. Note: Please note that Table 2 in this section provides summary data for 46 BBJ (BioBank Japan) datasets, which encompass all the BBJ data used in our published studies. 1 , 2 Readers may choose to download a subset of the summaries based on their specific needs. Table 2 The download links for BBJ summary data BBJ-Trait-summary URL Rheumatoid arthritis in Asian http://jenger.riken.jp:8080/download/Rheumatoid_Arthritis_Asian Open-angle glaucoma http://jenger.riken.jp:8080/download/Open-angle_glaucoma_autosome Type 2 Diabetes http://jenger.riken.jp:8080/download/Type_2_Diabetes Smoking initiation http://jenger.riken.jp:8080/download/Smoking_initiation_autosomes_all Smoking cessation http://jenger.riken.jp:8080/download/Smoking_cessation_autosomes_all Arrhythmia http://jenger.riken.jp:8080/download/Arrhythmia Asthma http://jenger.riken.jp:8080/download/Asthma Atopic dermatitis http://jenger.riken.jp:8080/download/Atopic_dermatitis Biliary tract cancer http://jenger.riken.jp:8080/download/Biliary_tract_cancer Cataract http://jenger.riken.jp:8080/download/Cataract Cerebral aneurysm http://jenger.riken.jp:8080/download/Cerebral_aneurysm Cervical cancer http://jenger.riken.jp:8080/download/Cervical_cancer Chronic hepatitis B http://jenger.riken.jp:8080/download/Chronic_hepatitis_B Chronic hepatitis C http://jenger.riken.jp:8080/download/Chronic_hepatitis_C Chronic obstructive pulmonary disease http://jenger.riken.jp:8080/download/Chronic_obstructive_pulmonary_disease Cirrhosis http://jenger.riken.jp:8080/download/Cirrhosis Congestive heart failure http://jenger.riken.jp:8080/download/Congestive_heart_failure Drug eruption http://jenger.riken.jp:8080/download/Drug_eruption Endometrial cancer http://jenger.riken.jp:8080/download/Endometrial_cancer Endometriosis http://jenger.riken.jp:8080/download/Endometriosis Epilepsy http://jenger.riken.jp:8080/download/Epilepsy Esophageal cancer http://jenger.riken.jp:8080/download/Esophageal_cancer Gastric cancer http://jenger.riken.jp:8080/download/Gastric_cancer Glaucoma http://jenger.riken.jp:8080/download/Glaucoma Graves disease http://jenger.riken.jp:8080/download/Graves_disease Hematological malignancy http://jenger.riken.jp:8080/download/Hematological_malignancy Interstitial lung disease http://jenger.riken.jp:8080/download/Interstitial_lung_disease Ischemic stroke http://jenger.riken.jp:8080/download/Ischemic_stroke Keloid http://jenger.riken.jp:8080/download/Keloid Lung cancer http://jenger.riken.jp:8080/download/Lung_cancer Nephrotic syndrome http://jenger.riken.jp:8080/download/Nephrotic_syndrome Osteoporosis http://jenger.riken.jp:8080/download/Osteoporosis Ovarian cancer http://jenger.riken.jp:8080/download/Ovarian_cancer Pancreatic cancer http://jenger.riken.jp:8080/download/Pancreatic_cancer Periodontal disease http://jenger.riken.jp:8080/download/Periodontal_disease Peripheral artery disease http://jenger.riken.jp:8080/download/Peripheral_artery_disease Pollinosis http://jenger.riken.jp:8080/download/Pollinosis Prostate cancer http://jenger.riken.jp:8080/download/Prostate_cancer Pulmonary tuberculosis http://jenger.riken.jp:8080/download/Pulmonary_tuberculosis Rheumatoid arthritis http://jenger.riken.jp:8080/download/Rheumatoid_arthritis Type 2 diabetes http://jenger.riken.jp:8080/download/Type_2_diabetes Urolithiasis http://jenger.riken.jp:8080/download/Urolithiasis Uterine fibroids http://jenger.riken.jp:8080/download/Uterine_fibroids hepatocellular carcinoma http://jenger.riken.jp:8080/download/hepatocellular_carcinoma Breast cancer http://jenger.riken.jp:8080/download/Breast_cancer Coronary artery disease http://jenger.riken.jp:8080/download/Coronary_artery_disease File sizes and documentation are available at: http://jenger.riken.jp/en/result . Refer to Table 2 for the download links to the summary data. Note: Please note that Table 2 in this section provides summary data for 46 BBJ (BioBank Japan) datasets, which encompass all the BBJ data used in our published studies. 1 , 2 Readers may choose to download a subset of the summaries based on their specific needs. Table 2 The download links for BBJ summary data BBJ-Trait-summary URL Rheumatoid arthritis in Asian http://jenger.riken.jp:8080/download/Rheumatoid_Arthritis_Asian Open-angle glaucoma http://jenger.riken.jp:8080/download/Open-angle_glaucoma_autosome Type 2 Diabetes http://jenger.riken.jp:8080/download/Type_2_Diabetes Smoking initiation http://jenger.riken.jp:8080/download/Smoking_initiation_autosomes_all Smoking cessation http://jenger.riken.jp:8080/download/Smoking_cessation_autosomes_all Arrhythmia http://jenger.riken.jp:8080/download/Arrhythmia Asthma http://jenger.riken.jp:8080/download/Asthma Atopic dermatitis http://jenger.riken.jp:8080/download/Atopic_dermatitis Biliary tract cancer http://jenger.riken.jp:8080/download/Biliary_tract_cancer Cataract http://jenger.riken.jp:8080/download/Cataract Cerebral aneurysm http://jenger.riken.jp:8080/download/Cerebral_aneurysm Cervical cancer http://jenger.riken.jp:8080/download/Cervical_cancer Chronic hepatitis B http://jenger.riken.jp:8080/download/Chronic_hepatitis_B Chronic hepatitis C http://jenger.riken.jp:8080/download/Chronic_hepatitis_C Chronic obstructive pulmonary disease http://jenger.riken.jp:8080/download/Chronic_obstructive_pulmonary_disease Cirrhosis http://jenger.riken.jp:8080/download/Cirrhosis Congestive heart failure http://jenger.riken.jp:8080/download/Congestive_heart_failure Drug eruption http://jenger.riken.jp:8080/download/Drug_eruption Endometrial cancer http://jenger.riken.jp:8080/download/Endometrial_cancer Endometriosis http://jenger.riken.jp:8080/download/Endometriosis Epilepsy http://jenger.riken.jp:8080/download/Epilepsy Esophageal cancer http://jenger.riken.jp:8080/download/Esophageal_cancer Gastric cancer http://jenger.riken.jp:8080/download/Gastric_cancer Glaucoma http://jenger.riken.jp:8080/download/Glaucoma Graves disease http://jenger.riken.jp:8080/download/Graves_disease Hematological malignancy http://jenger.riken.jp:8080/download/Hematological_malignancy Interstitial lung disease http://jenger.riken.jp:8080/download/Interstitial_lung_disease Ischemic stroke http://jenger.riken.jp:8080/download/Ischemic_stroke Keloid http://jenger.riken.jp:8080/download/Keloid Lung cancer http://jenger.riken.jp:8080/download/Lung_cancer Nephrotic syndrome http://jenger.riken.jp:8080/download/Nephrotic_syndrome Osteoporosis http://jenger.riken.jp:8080/download/Osteoporosis Ovarian cancer http://jenger.riken.jp:8080/download/Ovarian_cancer Pancreatic cancer http://jenger.riken.jp:8080/download/Pancreatic_cancer Periodontal disease http://jenger.riken.jp:8080/download/Periodontal_disease Peripheral artery disease http://jenger.riken.jp:8080/download/Peripheral_artery_disease Pollinosis http://jenger.riken.jp:8080/download/Pollinosis Prostate cancer http://jenger.riken.jp:8080/download/Prostate_cancer Pulmonary tuberculosis http://jenger.riken.jp:8080/download/Pulmonary_tuberculosis Rheumatoid arthritis http://jenger.riken.jp:8080/download/Rheumatoid_arthritis Type 2 diabetes http://jenger.riken.jp:8080/download/Type_2_diabetes Urolithiasis http://jenger.riken.jp:8080/download/Urolithiasis Uterine fibroids http://jenger.riken.jp:8080/download/Uterine_fibroids hepatocellular carcinoma http://jenger.riken.jp:8080/download/hepatocellular_carcinoma Breast cancer http://jenger.riken.jp:8080/download/Breast_cancer Coronary artery disease http://jenger.riken.jp:8080/download/Coronary_artery_disease Note: Please note that Table 2 in this section provides summary data for 46 BBJ (BioBank Japan) datasets, which encompass all the BBJ data used in our published studies. 1 , 2 Readers may choose to download a subset of the summaries based on their specific needs. The download links for BBJ summary data Download the expression files for transcriptome-wide association study (TWAS) study. a. Here, we used the GTEx V8 Cross-Tissue Imputation (UTMOST) Prediction Models from PredictDB as an example. b. Download and documentation are available at: https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . c. Move all of the downloaded files to the directory ‘/code/Supporting_data/GTEx’. Here, we used the GTEx V8 Cross-Tissue Imputation (UTMOST) Prediction Models from PredictDB as an example. Download and documentation are available at: https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . Move all of the downloaded files to the directory ‘/code/Supporting_data/GTEx’. Download and order the data files for drug-target analysis. a. Order anatomical therapeutic chemical (ATC) classification index Note: We need to obtain the ATC classification index from the WHO Collaborating Centre for Drug Statistics Methodology, which is available at: https://atcddd.fhi.no/atc_ddd_index_and_guidelines/order/ . b. Download the drug-gene interaction data from The Drug Gene Interaction Database: https://www.dgidb.org/data/latest/interactions.tsv . c. Download the MAGMA software dependencies: https://cncr.nl/research/magma/ . d. Move all of the downloaded files to the directory ‘/code/Supporting_data/magma’. Order anatomical therapeutic chemical (ATC) classification index Note: We need to obtain the ATC classification index from the WHO Collaborating Centre for Drug Statistics Methodology, which is available at: https://atcddd.fhi.no/atc_ddd_index_and_guidelines/order/ . Note: We need to obtain the ATC classification index from the WHO Collaborating Centre for Drug Statistics Methodology, which is available at: https://atcddd.fhi.no/atc_ddd_index_and_guidelines/order/ . Download the drug-gene interaction data from The Drug Gene Interaction Database: https://www.dgidb.org/data/latest/interactions.tsv . Download the MAGMA software dependencies: https://cncr.nl/research/magma/ . Move all of the downloaded files to the directory ‘/code/Supporting_data/magma’.

The protocol provides practical tutorials for most analytical methodologies detailed in two published articles, 1 , 2 primarily focusing on the analysis of ultra-low depth sequencing data from large sample sizes. The process is divided into two parts: the first involves routine sequencing data analysis and genotype imputation, resulting in a VCF file containing genotype dosage data. The second part encompasses GWAS and post-GWAS tasks, which yield results suitable for visualizing Manhattan plots. Subsequent analyses in the post-GWAS phase predominantly rely on these results to estimate heritability for specific phenotypes and assess genetic correlations between pairs of phenotypes. The outcomes of Mendelian randomization analyses include a range of p -values and specific MR methodologies, allowing for the inference of causal relationships between phenotypes. Transcriptome-wide association studies produce a similar set of p -values, indicating associations between gene expression in various tissues and phenotypes. Finally, drug target analysis identifies potential drug targets for a given phenotype, as indicated by the area under the curve and associated p -values.

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Huanhuan Zhu ( [email protected] ). For technical or coding issues, please contact the technical contact, Jingyu Zeng ( [email protected] ). This study did not generate new unique reagents. The code is available on GitHub at: https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis .

A significant challenge associated with this protocol is the extensive computational time and resources required to process the large sample sizes typically involved, as indicated by previous studies with sample sizes around 30,000. 1 , 2 Readers intending to undertake such analyses should ensure they have adequate computational resources. Furthermore, the workflow employed in this study has not been evaluated on ultra-low depth data beyond non-invasive prenatal testing (NIPT), as noted in the current literature; therefore, no such analyses have been published. Researchers are encouraged to adjust parameters independently to meet their specific analytical needs. Lastly, it is important to note that this workflow is limited to the analysis of point mutations; structural variants and insertion-deletion mutations fall outside the scope of this protocol.

A code repository is available on GitHub ( https://github.com/Mirror-fish/Ultra-Low-Depth-Sequencing-Genetic-Analysis ) containing all the necessary code. Before getting started, please download the code. Although the code can be organized in a workflow structure, it is not presented in that format here for clarity. Readers are encouraged to modify the workflow as needed. Before executing any code file, it is essential to review and adjust the parameters to ensure they align with local conditions. Timing: ∼4 hours per sample Timing: ∼4 hours per sample The primary objective of this step is to perform quality control and alignment on the raw FASTQ sequencing data obtained. As a result, each sample will yield a bam file that has been aligned, deduplicated, and subjected to Base Quality Score Recalibration (BQSR). In this section, readers have to manually modify portions of the code according to their specific circumstances. 1. Modify the arguments of workflow_bwa.sh. a. Specify the software and the absolute paths of the required source data at the beginning of the script 1.1.Workflow_bwa.sh. b. Insert the path of the GATK bundle into the folder name specified in the section “ download dependencies for analysis ”. ### static path for tools and reference on local disk. ### Tools gatk=/path_to/gatk-4.0.4.0/gatk java=/path_to/jdk1.8.0_131/bin/java bwa=/path_to/bwa/0.7.16/bwa samtools=/path_to/samtools fastp=./path_to/fastp-0.23.2/fastp ### Reference Genome hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz ### GATK bundle gatk_bundle_dir=/path_to/Required_reference_data/gatk-bundle #Check the sequencing platform platform=ILLUMINA #where to perform the analysis, the folder's absolute path out_path=/path_to/1.Quality_control_and_alignment_of_the_FASTQ_files Note: In our published work, 1 , 2 we used single-end sequencing data. This protocol also includes code tailored for single-end sequencing. For different sequencing types, the BWA software employs distinct algorithms for read alignment. For more details, please refer to the BWA manual: https://bio-bwa.sourceforge.net/bwa.shtml . Note: Please verify the platform that generated the sequencing files, as different sequencing platforms require specific names in the 'platform' argument. See https://samtools.github.io/hts-specs/SAMv1.pdf for details. 2. Start the alignment and quality control. Note: In this section, we will handle the analysis scripts for all samples, as they will be aligned and quality-controlled using the same script. a. Input the different sample names and their corresponding FASTQ files and execute the script. sh /path_to/1.1.Workflow_bwa.sh Sample_name1 Fastq_file1 sh /path_to/1.1.Workflow_bwa.sh Sample_name2 Fastq_file2 ... Note: We typically avoid running this script directly due to its length. Instead, you can adjust the execution mode based on the machine you are using, such as splitting the script into multiple tasks and submitting them individually for execution. Modify the arguments of workflow_bwa.sh. a. Specify the software and the absolute paths of the required source data at the beginning of the script 1.1.Workflow_bwa.sh. b. Insert the path of the GATK bundle into the folder name specified in the section “ download dependencies for analysis ”. ### static path for tools and reference on local disk. ### Tools gatk=/path_to/gatk-4.0.4.0/gatk java=/path_to/jdk1.8.0_131/bin/java bwa=/path_to/bwa/0.7.16/bwa samtools=/path_to/samtools fastp=./path_to/fastp-0.23.2/fastp ### Reference Genome hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz ### GATK bundle gatk_bundle_dir=/path_to/Required_reference_data/gatk-bundle #Check the sequencing platform platform=ILLUMINA #where to perform the analysis, the folder's absolute path out_path=/path_to/1.Quality_control_and_alignment_of_the_FASTQ_files Note: In our published work, 1 , 2 we used single-end sequencing data. This protocol also includes code tailored for single-end sequencing. For different sequencing types, the BWA software employs distinct algorithms for read alignment. For more details, please refer to the BWA manual: https://bio-bwa.sourceforge.net/bwa.shtml . Note: Please verify the platform that generated the sequencing files, as different sequencing platforms require specific names in the 'platform' argument. See https://samtools.github.io/hts-specs/SAMv1.pdf for details. Specify the software and the absolute paths of the required source data at the beginning of the script 1.1.Workflow_bwa.sh. Insert the path of the GATK bundle into the folder name specified in the section “ download dependencies for analysis ”. ### static path for tools and reference on local disk. ### Tools gatk=/path_to/gatk-4.0.4.0/gatk java=/path_to/jdk1.8.0_131/bin/java bwa=/path_to/bwa/0.7.16/bwa samtools=/path_to/samtools fastp=./path_to/fastp-0.23.2/fastp ### Reference Genome hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz ### GATK bundle gatk_bundle_dir=/path_to/Required_reference_data/gatk-bundle #Check the sequencing platform platform=ILLUMINA #where to perform the analysis, the folder's absolute path out_path=/path_to/1.Quality_control_and_alignment_of_the_FASTQ_files Note: In our published work, 1 , 2 we used single-end sequencing data. This protocol also includes code tailored for single-end sequencing. For different sequencing types, the BWA software employs distinct algorithms for read alignment. For more details, please refer to the BWA manual: https://bio-bwa.sourceforge.net/bwa.shtml . Note: Please verify the platform that generated the sequencing files, as different sequencing platforms require specific names in the 'platform' argument. See https://samtools.github.io/hts-specs/SAMv1.pdf for details. ### static path for tools and reference on local disk. ### Tools gatk=/path_to/gatk-4.0.4.0/gatk java=/path_to/jdk1.8.0_131/bin/java bwa=/path_to/bwa/0.7.16/bwa samtools=/path_to/samtools fastp=./path_to/fastp-0.23.2/fastp ### Reference Genome hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz ### GATK bundle gatk_bundle_dir=/path_to/Required_reference_data/gatk-bundle #Check the sequencing platform platform=ILLUMINA #where to perform the analysis, the folder's absolute path out_path=/path_to/1.Quality_control_and_alignment_of_the_FASTQ_files Note: In our published work, 1 , 2 we used single-end sequencing data. This protocol also includes code tailored for single-end sequencing. For different sequencing types, the BWA software employs distinct algorithms for read alignment. For more details, please refer to the BWA manual: https://bio-bwa.sourceforge.net/bwa.shtml . Note: Please verify the platform that generated the sequencing files, as different sequencing platforms require specific names in the 'platform' argument. See https://samtools.github.io/hts-specs/SAMv1.pdf for details. Start the alignment and quality control. Note: In this section, we will handle the analysis scripts for all samples, as they will be aligned and quality-controlled using the same script. a. Input the different sample names and their corresponding FASTQ files and execute the script. sh /path_to/1.1.Workflow_bwa.sh Sample_name1 Fastq_file1 sh /path_to/1.1.Workflow_bwa.sh Sample_name2 Fastq_file2 ... Note: We typically avoid running this script directly due to its length. Instead, you can adjust the execution mode based on the machine you are using, such as splitting the script into multiple tasks and submitting them individually for execution. Note: In this section, we will handle the analysis scripts for all samples, as they will be aligned and quality-controlled using the same script. Input the different sample names and their corresponding FASTQ files and execute the script. sh /path_to/1.1.Workflow_bwa.sh Sample_name1 Fastq_file1 sh /path_to/1.1.Workflow_bwa.sh Sample_name2 Fastq_file2 ... Note: We typically avoid running this script directly due to its length. Instead, you can adjust the execution mode based on the machine you are using, such as splitting the script into multiple tasks and submitting them individually for execution. sh /path_to/1.1.Workflow_bwa.sh Sample_name1 Fastq_file1 sh /path_to/1.1.Workflow_bwa.sh Sample_name2 Fastq_file2 ... Note: We typically avoid running this script directly due to its length. Instead, you can adjust the execution mode based on the machine you are using, such as splitting the script into multiple tasks and submitting them individually for execution. Timing: 2–3 weeks (depends on the sample size) Timing: 2–3 weeks (depends on the sample size) In this step, genotype imputation is performed using the BAM file obtained from the previous step. Initially, the entire genome is divided into uniformly sized segments. Subsequently, genotype imputation is carried out on each segment individually. Ultimately, each segment yields a corresponding '.vcf' file. These files are then merged into a single comprehensive '.vcf' file using BCFtools, thus completing the process. Note: Generally, the larger the sample size, the more time-consuming the process becomes. In our published work, 1 , 2 we used approximately 30,000 samples and divided the genome into segments of 5,000,000 bp, resulting in about 620 tasks for genotype imputation. The analysis, conducted on 4 CPUs, took roughly 12 h per task. Since our cluster tasks run in parallel, and factoring in queuing time, the entire genotype imputation process (including error correction and reruns, as discussed in the troubleshooting section) took approximately 3 weeks to complete. 3. Modify the arguments of 2.1.Create_bash_STITCH.py. Note: Given the large sample size used in STITCH, failing to segment chromosome intervals in advance may lead to insufficient computing resources. Therefore, we provide a script that divides chromosomes into smaller chunks. Example input files for this part are provided. You can find additional examples in the folder ‘2.Genotype_imputation_using_STITCH’. a. Adjust the parameters in the script 2.1.Create_bash_STITCH.py before executing it. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/2.Genotype_imputation_using_STITCH/bash' #sample names and bam paths, need to match each other sample_namelist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' bamlist='/path_to/2.Genotype_imputation_using_STITCH/bamid.list' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #the path to the STITCH code STITCH_script ='/path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh' Note: The first parameter, chunk_size, represents the size of the segments into which we want to divide the genome. Larger chunks will require more computing resources and time. It is advisable to test different chunk sizes before starting to determine the optimal parameters. Note: We also need to provide the paths to two files: sample_namelist and bamlist. The former contains the names of all samples used in this round of the STITCH operation, while the latter lists the BAM file paths for the corresponding samples. Please ensure that the contents of both files correspond one-to-one. Note: The parameter CHR_length refers to a file that specifies the lengths of each chromosome (e.g., chr1, 248956422 bp). We have included a hg38 version of this file in our code package. If you are using a different version of the reference genome, please replace this file with the corresponding chromosome length file. The STITCH_script refers to 2.2.run_STITCH.sh, which we will edit in the next step. b. Once all the parameters are set, run the code using Python (version 3.x): Note: After completing the previous step, you should find the bash and output folders in the STITCH working directory. In the bash folder, there will be a file titled Chr#.sh, and its contents should resemble the following: #Example of the outcome bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_1_5000000 chr1 1 5000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_5000001_10000000 chr1 5000001 10000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_10000001_15000000 chr1 10000001 15000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_15000001_20000000 chr1 15000001 20000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_20000001_25000000 chr1 20000001 25000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_25000001_30000000 chr1 25000001 30000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_30000001_35000000 chr1 30000001 35000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_35000001_40000000 chr1 35000001 40000000 Note: Please note that you should not execute the code directly at this stage. First, some adjustments need to be made to the code within the 2.2.run_STITCH.sh file. It is also important to carefully review the output files for any errors. 4. Modify the arguments of 2.2.run_STITCH.sh. a. Specify the parameters for 2.2.run_STITCH.sh and execute the script. #Specify the arguments before execute the code related #path for bcftools bcftools=/path_to/bcftools ##sample names and bam paths, the sample file as 2.1.Create_bash_STITCH.py sample_namelist=/path_to/2.Genotype_imputation_using_STITCH/bampath.list bamlist=/path_to/2.Genotype_imputation_using_STITCH/bamid.list #output folder outdir=/path_to/2.Genotype_imputation_using_STITCH/output #Reference Genome, should be the same as 01.workflow_bwa.sh hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz #Before starting STITCH, please make sure you have correctly installed STITCH in your environment #If correctly installed, you can find the STITCH.R in the directory you install STITCH STITCH_R=/path_to/STITCH.R Note: The first parameter (bcftools) should be the path to the previously installed BCFtools. Additionally, please ensure that the paths for sample_namelist and bamlist match those used in 2.1.Create_bash_STITCH.py. Note: The parameter hg38 refers to the path for the reference genome; please ensure it matches the reference genome used in the ‘Quality Control and Alignment of the FASTQ Files’. Note: The last parameter STITCH_R , which should point to the folder where STITCH was previously installed. Please fill in the corresponding file path. b. Specify the parameters (‘-- K’ and ‘-- nGen’) for the imputation execution. ###Specify the arguments before execute! time /path_to/Rscript $STITCH_R \ --outputdir ${outdir}/${chr}/${chr_region} \ --bamlist $bamlist \ --sampleNames_file $sample_namelist \ --reference $hg38 \ --K 10 --nGen 16000 --nCores 1 \ --regionStart $start --regionEnd $end --chr $chr \ --buffer 250000 Note: According to the STITCH official documentation, “K is the number of ancestral haplotypes in the model. A larger K improves imputation accuracy for larger samples and higher coverages, but it also increases computation time, and accuracy may decrease with lower coverage“. 14 Therefore, it is recommended to test different ‘--K’ values to identify the most suitable one before beginning the analysis. The ‘--nGen’ value can be calculated using the formula: 4 ∗ Sample size / K, so this parameter can be determined once the sample size and ‘--K’ value are set. Note: The parameter ‘--nCores’ refers to the number of CPU cores allocated for each running job; you can adjust this based on your specific situation. Note: The ‘--buffer’ parameter is primarily used to account for potential connections between chromosomal segments. During imputation of a specified chromosomal region, STITCH extends the region by the length defined by ‘--buffer’ at both ends to preserve any potential linkage disequilibrium (LD) information. The buffer size should be set based on the LD characteristics of the data being used. Alternatives: STITCH can perform imputation without a reference panel; however, incorporating one reference panel can influence the number and quality of loci in the final results. For instructions on how to add a reference panel, please consult the official STITCH documentation ( https://github.com/rwdavies/STITCH?tab=readme-ov-file#about-reference-panels- ). 5. Start genotype imputation using STITCH. a. Run the ‘Chr#.sh’ script generated by 2.1.Create_bash_STITCH.py. b. Check the analysis results for each chromosome chunk in the output folder. Note: We typically avoid running these scripts directly due to their length. Instead, you can adjust the execution mode based on your machine, such as by splitting the script into multiple tasks and submitting them for individual execution. 6. Merge all genotype imputation results using BCFtools. a. Use the ls command to locate all completed files: ls /path_to/2.Genotype_imputation_using_STITCH/output/∗/∗.vcf.gz > /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list b. Specify the arguments in 2.3.Merge.sh and then execute the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list \ -Oz \ --threads 4 \ -o /path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz Note: The ‘--threads’ parameter allows you to adjust the number of parallel threads in BCFtools, which can enhance the execution speed of this step. The optimal setting may vary depending on the capabilities of the machine you are using. CRITICAL: The ‘.vcf’ file generated after this step serves as the basis for various analyses, including GWAS. It is important to note that this ‘.vcf’ file differs from standard ‘.vcf’ files in that the genotypes are represented using genotype dosage (DS). This distinction should be considered when conducting further analyses. Note: Generally, the larger the sample size, the more time-consuming the process becomes. In our published work, 1 , 2 we used approximately 30,000 samples and divided the genome into segments of 5,000,000 bp, resulting in about 620 tasks for genotype imputation. The analysis, conducted on 4 CPUs, took roughly 12 h per task. Since our cluster tasks run in parallel, and factoring in queuing time, the entire genotype imputation process (including error correction and reruns, as discussed in the troubleshooting section) took approximately 3 weeks to complete. Modify the arguments of 2.1.Create_bash_STITCH.py. Note: Given the large sample size used in STITCH, failing to segment chromosome intervals in advance may lead to insufficient computing resources. Therefore, we provide a script that divides chromosomes into smaller chunks. Example input files for this part are provided. You can find additional examples in the folder ‘2.Genotype_imputation_using_STITCH’. a. Adjust the parameters in the script 2.1.Create_bash_STITCH.py before executing it. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/2.Genotype_imputation_using_STITCH/bash' #sample names and bam paths, need to match each other sample_namelist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' bamlist='/path_to/2.Genotype_imputation_using_STITCH/bamid.list' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #the path to the STITCH code STITCH_script ='/path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh' Note: The first parameter, chunk_size, represents the size of the segments into which we want to divide the genome. Larger chunks will require more computing resources and time. It is advisable to test different chunk sizes before starting to determine the optimal parameters. Note: We also need to provide the paths to two files: sample_namelist and bamlist. The former contains the names of all samples used in this round of the STITCH operation, while the latter lists the BAM file paths for the corresponding samples. Please ensure that the contents of both files correspond one-to-one. Note: The parameter CHR_length refers to a file that specifies the lengths of each chromosome (e.g., chr1, 248956422 bp). We have included a hg38 version of this file in our code package. If you are using a different version of the reference genome, please replace this file with the corresponding chromosome length file. The STITCH_script refers to 2.2.run_STITCH.sh, which we will edit in the next step. b. Once all the parameters are set, run the code using Python (version 3.x): Note: After completing the previous step, you should find the bash and output folders in the STITCH working directory. In the bash folder, there will be a file titled Chr#.sh, and its contents should resemble the following: #Example of the outcome bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_1_5000000 chr1 1 5000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_5000001_10000000 chr1 5000001 10000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_10000001_15000000 chr1 10000001 15000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_15000001_20000000 chr1 15000001 20000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_20000001_25000000 chr1 20000001 25000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_25000001_30000000 chr1 25000001 30000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_30000001_35000000 chr1 30000001 35000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_35000001_40000000 chr1 35000001 40000000 Note: Please note that you should not execute the code directly at this stage. First, some adjustments need to be made to the code within the 2.2.run_STITCH.sh file. It is also important to carefully review the output files for any errors. Note: Given the large sample size used in STITCH, failing to segment chromosome intervals in advance may lead to insufficient computing resources. Therefore, we provide a script that divides chromosomes into smaller chunks. Example input files for this part are provided. You can find additional examples in the folder ‘2.Genotype_imputation_using_STITCH’. Adjust the parameters in the script 2.1.Create_bash_STITCH.py before executing it. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/2.Genotype_imputation_using_STITCH/bash' #sample names and bam paths, need to match each other sample_namelist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' bamlist='/path_to/2.Genotype_imputation_using_STITCH/bamid.list' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #the path to the STITCH code STITCH_script ='/path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh' Note: The first parameter, chunk_size, represents the size of the segments into which we want to divide the genome. Larger chunks will require more computing resources and time. It is advisable to test different chunk sizes before starting to determine the optimal parameters. Note: We also need to provide the paths to two files: sample_namelist and bamlist. The former contains the names of all samples used in this round of the STITCH operation, while the latter lists the BAM file paths for the corresponding samples. Please ensure that the contents of both files correspond one-to-one. Note: The parameter CHR_length refers to a file that specifies the lengths of each chromosome (e.g., chr1, 248956422 bp). We have included a hg38 version of this file in our code package. If you are using a different version of the reference genome, please replace this file with the corresponding chromosome length file. The STITCH_script refers to 2.2.run_STITCH.sh, which we will edit in the next step. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/2.Genotype_imputation_using_STITCH/bash' #sample names and bam paths, need to match each other sample_namelist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' bamlist='/path_to/2.Genotype_imputation_using_STITCH/bamid.list' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #the path to the STITCH code STITCH_script ='/path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh' Note: The first parameter, chunk_size, represents the size of the segments into which we want to divide the genome. Larger chunks will require more computing resources and time. It is advisable to test different chunk sizes before starting to determine the optimal parameters. Note: We also need to provide the paths to two files: sample_namelist and bamlist. The former contains the names of all samples used in this round of the STITCH operation, while the latter lists the BAM file paths for the corresponding samples. Please ensure that the contents of both files correspond one-to-one. Note: The parameter CHR_length refers to a file that specifies the lengths of each chromosome (e.g., chr1, 248956422 bp). We have included a hg38 version of this file in our code package. If you are using a different version of the reference genome, please replace this file with the corresponding chromosome length file. The STITCH_script refers to 2.2.run_STITCH.sh, which we will edit in the next step. Once all the parameters are set, run the code using Python (version 3.x): Note: After completing the previous step, you should find the bash and output folders in the STITCH working directory. In the bash folder, there will be a file titled Chr#.sh, and its contents should resemble the following: #Example of the outcome bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_1_5000000 chr1 1 5000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_5000001_10000000 chr1 5000001 10000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_10000001_15000000 chr1 10000001 15000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_15000001_20000000 chr1 15000001 20000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_20000001_25000000 chr1 20000001 25000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_25000001_30000000 chr1 25000001 30000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_30000001_35000000 chr1 30000001 35000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_35000001_40000000 chr1 35000001 40000000 Note: Please note that you should not execute the code directly at this stage. First, some adjustments need to be made to the code within the 2.2.run_STITCH.sh file. It is also important to carefully review the output files for any errors. Note: After completing the previous step, you should find the bash and output folders in the STITCH working directory. In the bash folder, there will be a file titled Chr#.sh, and its contents should resemble the following: #Example of the outcome bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_1_5000000 chr1 1 5000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_5000001_10000000 chr1 5000001 10000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_10000001_15000000 chr1 10000001 15000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_15000001_20000000 chr1 15000001 20000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_20000001_25000000 chr1 20000001 25000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_25000001_30000000 chr1 25000001 30000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_30000001_35000000 chr1 30000001 35000000 bash /Path_to/2.Genotype_imputation_using_STITCH/2.2.run_STITCH.sh chr1_35000001_40000000 chr1 35000001 40000000 Note: Please note that you should not execute the code directly at this stage. First, some adjustments need to be made to the code within the 2.2.run_STITCH.sh file. It is also important to carefully review the output files for any errors. Modify the arguments of 2.2.run_STITCH.sh. a. Specify the parameters for 2.2.run_STITCH.sh and execute the script. #Specify the arguments before execute the code related #path for bcftools bcftools=/path_to/bcftools ##sample names and bam paths, the sample file as 2.1.Create_bash_STITCH.py sample_namelist=/path_to/2.Genotype_imputation_using_STITCH/bampath.list bamlist=/path_to/2.Genotype_imputation_using_STITCH/bamid.list #output folder outdir=/path_to/2.Genotype_imputation_using_STITCH/output #Reference Genome, should be the same as 01.workflow_bwa.sh hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz #Before starting STITCH, please make sure you have correctly installed STITCH in your environment #If correctly installed, you can find the STITCH.R in the directory you install STITCH STITCH_R=/path_to/STITCH.R Note: The first parameter (bcftools) should be the path to the previously installed BCFtools. Additionally, please ensure that the paths for sample_namelist and bamlist match those used in 2.1.Create_bash_STITCH.py. Note: The parameter hg38 refers to the path for the reference genome; please ensure it matches the reference genome used in the ‘Quality Control and Alignment of the FASTQ Files’. Note: The last parameter STITCH_R , which should point to the folder where STITCH was previously installed. Please fill in the corresponding file path. b. Specify the parameters (‘-- K’ and ‘-- nGen’) for the imputation execution. ###Specify the arguments before execute! time /path_to/Rscript $STITCH_R \ --outputdir ${outdir}/${chr}/${chr_region} \ --bamlist $bamlist \ --sampleNames_file $sample_namelist \ --reference $hg38 \ --K 10 --nGen 16000 --nCores 1 \ --regionStart $start --regionEnd $end --chr $chr \ --buffer 250000 Note: According to the STITCH official documentation, “K is the number of ancestral haplotypes in the model. A larger K improves imputation accuracy for larger samples and higher coverages, but it also increases computation time, and accuracy may decrease with lower coverage“. 14 Therefore, it is recommended to test different ‘--K’ values to identify the most suitable one before beginning the analysis. The ‘--nGen’ value can be calculated using the formula: 4 ∗ Sample size / K, so this parameter can be determined once the sample size and ‘--K’ value are set. Note: The parameter ‘--nCores’ refers to the number of CPU cores allocated for each running job; you can adjust this based on your specific situation. Note: The ‘--buffer’ parameter is primarily used to account for potential connections between chromosomal segments. During imputation of a specified chromosomal region, STITCH extends the region by the length defined by ‘--buffer’ at both ends to preserve any potential linkage disequilibrium (LD) information. The buffer size should be set based on the LD characteristics of the data being used. Alternatives: STITCH can perform imputation without a reference panel; however, incorporating one reference panel can influence the number and quality of loci in the final results. For instructions on how to add a reference panel, please consult the official STITCH documentation ( https://github.com/rwdavies/STITCH?tab=readme-ov-file#about-reference-panels- ). Specify the parameters for 2.2.run_STITCH.sh and execute the script. #Specify the arguments before execute the code related #path for bcftools bcftools=/path_to/bcftools ##sample names and bam paths, the sample file as 2.1.Create_bash_STITCH.py sample_namelist=/path_to/2.Genotype_imputation_using_STITCH/bampath.list bamlist=/path_to/2.Genotype_imputation_using_STITCH/bamid.list #output folder outdir=/path_to/2.Genotype_imputation_using_STITCH/output #Reference Genome, should be the same as 01.workflow_bwa.sh hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz #Before starting STITCH, please make sure you have correctly installed STITCH in your environment #If correctly installed, you can find the STITCH.R in the directory you install STITCH STITCH_R=/path_to/STITCH.R Note: The first parameter (bcftools) should be the path to the previously installed BCFtools. Additionally, please ensure that the paths for sample_namelist and bamlist match those used in 2.1.Create_bash_STITCH.py. Note: The parameter hg38 refers to the path for the reference genome; please ensure it matches the reference genome used in the ‘Quality Control and Alignment of the FASTQ Files’. Note: The last parameter STITCH_R , which should point to the folder where STITCH was previously installed. Please fill in the corresponding file path. #Specify the arguments before execute the code related #path for bcftools bcftools=/path_to/bcftools ##sample names and bam paths, the sample file as 2.1.Create_bash_STITCH.py sample_namelist=/path_to/2.Genotype_imputation_using_STITCH/bampath.list bamlist=/path_to/2.Genotype_imputation_using_STITCH/bamid.list #output folder outdir=/path_to/2.Genotype_imputation_using_STITCH/output #Reference Genome, should be the same as 01.workflow_bwa.sh hg38=/path_to/Required_reference_data/hg38/hg38.fa.masked.gz #Before starting STITCH, please make sure you have correctly installed STITCH in your environment #If correctly installed, you can find the STITCH.R in the directory you install STITCH STITCH_R=/path_to/STITCH.R Note: The first parameter (bcftools) should be the path to the previously installed BCFtools. Additionally, please ensure that the paths for sample_namelist and bamlist match those used in 2.1.Create_bash_STITCH.py. Note: The parameter hg38 refers to the path for the reference genome; please ensure it matches the reference genome used in the ‘Quality Control and Alignment of the FASTQ Files’. Note: The last parameter STITCH_R , which should point to the folder where STITCH was previously installed. Please fill in the corresponding file path. Specify the parameters (‘-- K’ and ‘-- nGen’) for the imputation execution. ###Specify the arguments before execute! time /path_to/Rscript $STITCH_R \ --outputdir ${outdir}/${chr}/${chr_region} \ --bamlist $bamlist \ --sampleNames_file $sample_namelist \ --reference $hg38 \ --K 10 --nGen 16000 --nCores 1 \ --regionStart $start --regionEnd $end --chr $chr \ --buffer 250000 Note: According to the STITCH official documentation, “K is the number of ancestral haplotypes in the model. A larger K improves imputation accuracy for larger samples and higher coverages, but it also increases computation time, and accuracy may decrease with lower coverage“. 14 Therefore, it is recommended to test different ‘--K’ values to identify the most suitable one before beginning the analysis. The ‘--nGen’ value can be calculated using the formula: 4 ∗ Sample size / K, so this parameter can be determined once the sample size and ‘--K’ value are set. Note: The parameter ‘--nCores’ refers to the number of CPU cores allocated for each running job; you can adjust this based on your specific situation. Note: The ‘--buffer’ parameter is primarily used to account for potential connections between chromosomal segments. During imputation of a specified chromosomal region, STITCH extends the region by the length defined by ‘--buffer’ at both ends to preserve any potential linkage disequilibrium (LD) information. The buffer size should be set based on the LD characteristics of the data being used. Alternatives: STITCH can perform imputation without a reference panel; however, incorporating one reference panel can influence the number and quality of loci in the final results. For instructions on how to add a reference panel, please consult the official STITCH documentation ( https://github.com/rwdavies/STITCH?tab=readme-ov-file#about-reference-panels- ). ###Specify the arguments before execute! time /path_to/Rscript $STITCH_R \ --outputdir ${outdir}/${chr}/${chr_region} \ --bamlist $bamlist \ --sampleNames_file $sample_namelist \ --reference $hg38 \ --K 10 --nGen 16000 --nCores 1 \ --regionStart $start --regionEnd $end --chr $chr \ --buffer 250000 Note: According to the STITCH official documentation, “K is the number of ancestral haplotypes in the model. A larger K improves imputation accuracy for larger samples and higher coverages, but it also increases computation time, and accuracy may decrease with lower coverage“. 14 Therefore, it is recommended to test different ‘--K’ values to identify the most suitable one before beginning the analysis. The ‘--nGen’ value can be calculated using the formula: 4 ∗ Sample size / K, so this parameter can be determined once the sample size and ‘--K’ value are set. Note: The parameter ‘--nCores’ refers to the number of CPU cores allocated for each running job; you can adjust this based on your specific situation. Note: The ‘--buffer’ parameter is primarily used to account for potential connections between chromosomal segments. During imputation of a specified chromosomal region, STITCH extends the region by the length defined by ‘--buffer’ at both ends to preserve any potential linkage disequilibrium (LD) information. The buffer size should be set based on the LD characteristics of the data being used. Alternatives: STITCH can perform imputation without a reference panel; however, incorporating one reference panel can influence the number and quality of loci in the final results. For instructions on how to add a reference panel, please consult the official STITCH documentation ( https://github.com/rwdavies/STITCH?tab=readme-ov-file#about-reference-panels- ). Start genotype imputation using STITCH. a. Run the ‘Chr#.sh’ script generated by 2.1.Create_bash_STITCH.py. b. Check the analysis results for each chromosome chunk in the output folder. Run the ‘Chr#.sh’ script generated by 2.1.Create_bash_STITCH.py. Check the analysis results for each chromosome chunk in the output folder. Note: We typically avoid running these scripts directly due to their length. Instead, you can adjust the execution mode based on your machine, such as by splitting the script into multiple tasks and submitting them for individual execution. Merge all genotype imputation results using BCFtools. a. Use the ls command to locate all completed files: ls /path_to/2.Genotype_imputation_using_STITCH/output/∗/∗.vcf.gz > /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list b. Specify the arguments in 2.3.Merge.sh and then execute the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list \ -Oz \ --threads 4 \ -o /path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz Note: The ‘--threads’ parameter allows you to adjust the number of parallel threads in BCFtools, which can enhance the execution speed of this step. The optimal setting may vary depending on the capabilities of the machine you are using. CRITICAL: The ‘.vcf’ file generated after this step serves as the basis for various analyses, including GWAS. It is important to note that this ‘.vcf’ file differs from standard ‘.vcf’ files in that the genotypes are represented using genotype dosage (DS). This distinction should be considered when conducting further analyses. Use the ls command to locate all completed files: ls /path_to/2.Genotype_imputation_using_STITCH/output/∗/∗.vcf.gz > /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list ls /path_to/2.Genotype_imputation_using_STITCH/output/∗/∗.vcf.gz > /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list Specify the arguments in 2.3.Merge.sh and then execute the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list \ -Oz \ --threads 4 \ -o /path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz Note: The ‘--threads’ parameter allows you to adjust the number of parallel threads in BCFtools, which can enhance the execution speed of this step. The optimal setting may vary depending on the capabilities of the machine you are using. CRITICAL: The ‘.vcf’ file generated after this step serves as the basis for various analyses, including GWAS. It is important to note that this ‘.vcf’ file differs from standard ‘.vcf’ files in that the genotypes are represented using genotype dosage (DS). This distinction should be considered when conducting further analyses. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/2.Genotype_imputation_using_STITCH/completed_vcf.list \ -Oz \ --threads 4 \ -o /path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz Note: The ‘--threads’ parameter allows you to adjust the number of parallel threads in BCFtools, which can enhance the execution speed of this step. The optimal setting may vary depending on the capabilities of the machine you are using. CRITICAL: The ‘.vcf’ file generated after this step serves as the basis for various analyses, including GWAS. It is important to note that this ‘.vcf’ file differs from standard ‘.vcf’ files in that the genotypes are represented using genotype dosage (DS). This distinction should be considered when conducting further analyses. Timing: 1–2 days Timing: 1–2 days During this step, genome-wide association analysis (GWAS) will be conducted on the previously acquired population variant data (i.e., '.vcf' files). This section will begin by preparing covariates for the data: employing a dedicated single nucleotide variant detection software (BaseVarC) to perform single nucleotide variant detection on all '.bam' files. Subsequently, the data will be consolidated for the computation of a Principal Component Analysis (PCA). After completing the PCA, the results from the genotype imputation will be converted using PLINK2 to the appropriate format. Finally, the GWAS will be carried out using the designated phenotype file and covariate file. Visualization codes are provided for data summary. 7. Conduct segmented single-base mutation detection using BaseVarC and merge the results for principal component analysis (PCA). Note: The primary goal of this step is to perform variant detection across all samples, resulting in the creation of a ‘.vcf’ file that includes data from all samples. However, this ‘.vcf’ file is distinct from the one generated through genotype imputation. Its purpose is to provide a more comprehensive set of original variant sites for PCA analysis before genotype imputation. For subsequent GWAS analyses, the ‘.vcf’ file produced from genotype imputation will serve as the input for population variation. a. Modify the arguments of 3.1.BaseVarC_prep.py. Note: Due to the sample size considerations in STITCH, segmenting the processing can help conserve computing resources and time. Here, we provide a script that divides chromosomes into smaller chunks and generates scripts for BaseVarC. i. Specify the parameters in 3.1.BaseVarC_prep.py. ii. Execute the script to segment the chromosome regions. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/3.Genome-wide_association_study_using_PLINK2' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #sample names and bam paths, need to match each other bamlist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' #Reference Genome, should be the same as 01.workflow_bwa.sh ref_genome='/path_to/Required_reference_data/hg38/hg38.fa.masked.gz' #the path to path of BaseVarC Basevar='/path_to/BaseVarC' Note: It is recommended to ensure consistency with the parameters used previously for each setting here. b. Run the following command: Python 3.1.BaseVarC_prep.py. c. Execute the script generated. Note: In the directory ‘3.Genome-wide_association_study_using_PLINK2/bash’, you will find the relevant script named ‘CMD_basevar_chrx.sh’. You can refer to the methods previously used for running STITCH imputation to execute this script. d. Merge all BaseVarC results using BCFtools. You can use the ls command to locate all completed files: ls /path_to/3.Genome-wide_association_study_using_PLINK2/output/∗/∗.vcf.gz > /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list Alternatives: You can use the VCF list generated by 3.1.BaseVarC_prep.py. e. Specify the arguments in 3.2.Merge.sh, and then run the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list \ -Oz \ -o /path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz f. Modify the arguments of 3.3.Plink_load_forBaseVarC.sh. i. Enter the absolute path for PLINK2. ii. Fill in the corresponding parameters for the ‘.vcf’ file generated after executing 3.2.Merge.sh. iii. Run the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 3.2.Merge.sh merge_out=/path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz $plink2 \ --pca 5 \ --vcf $merge_out \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5 Note: In this step, we utilize the population single-base variation ‘.vcf’ file generated by BaseVarC for PCA calculations, rather than the ‘.vcf’ result file produced by STITCH. Note: In this step, we did not convert the ‘.vcf’ file to ‘.pgen’ format using PLINK2 before performing PCA calculations; therefore, rerunning this script will not save time. If you need to use this ‘.vcf’ data multiple times for PLINK2 operations, you can add the ‘--make-pgen’ option in the script and then replace the ‘--vcf’ parameter with ‘--pfile’ when recalculating PCA. Alternatives: In this protocol, the computations performed by PLINK2 are limited to 5 principal components (PCs). If the analysis involves a larger sample size or a more complex population structure, it is advisable to adjust the parameters for PCA. You can change the number of PCs calculated by using the ‘--pca’ option. 8. Perform GWAS using PLINK2. Note: In this step, we will conduct GWAS analysis using the PCA results obtained in the previous step as covariates, ultimately generating the data required to create the Manhattan plot. a. Load the pfile used in GWAS using PLINK2. i. Enter the corresponding parameters for the ‘.vcf’ file generated after completing 2.3.Merge.sh, along with the absolute path for PLINK2. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 2.3.Merge.sh merge_out=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz $plink2 \ --freq \ --hardy \ --make-pgen 'vzs' \ --set-all-var-ids chr@:# \ --vcf $merge_out dosage=DS \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH Note: In this step, we convert the ‘.vcf’ file generated by STITCH into a 'pgen’ file compatible with PLINK for GWAS. The generated ‘pgen’ file will serve as the genotype input for subsequent GWAS analyses. The ‘--make-pgen’ argument includes the ‘vzs’ parameter, which compresses the output genotype file to save space, making it particularly useful for large sample GWAS. b. Calculate the kinship among samples. i. Specify the arguments in ‘3.5.1.kinship_analysis.sh’ and execute the code: Note: In GWAS, statistical analyses based on general linear models require that samples be independent to ensure the reliability of the results. Samples with kinship relationships violate this assumption of independence and should be excluded from the dataset. In this step, we use PLINK2 to calculate the kinship among samples, enabling the exclusion of related samples in subsequent analyses. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --make-king triangle bin\ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship Note: The result of the ‘--make-king’ calculation is a kinship matrix file with the suffix ‘.king.bin’, which contains the pairwise kinship relationships between each sample. The kinship coefficient for duplicate samples and identical twins is 0.5. For first-degree relationships, such as parent-child or full siblings, the kinship coefficient is 0.25, and for second-degree relationships, it is 0.125, and so on. Alternatives: To exclude samples with first-degree relationships, add the ‘--king-cutoff /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship 0.177’ parameter when loading the pfile. Alternatively, you can use your preferred algorithm to identify the samples that need to be excluded based on the kinship matrix, and then remove these samples from the dataset using the ‘--remove’ flag. Note: 0.177 is the geometric mean of 0.25 and 0.125. When selecting a cutoff threshold, the geometric mean between two degrees of relationship is commonly used. Alternatives: In the context of GWAS, it's essential to note that the input samples for the generalized linear model must be independent of each other. However, when utilizing software designed for mixed linear models, such as GEMMA, 20 it is possible to incorporate a kinship matrix as a parameter for GWAS analysis. Interested readers are encouraged to conduct their analyses using relevant software. It is important to remember that the sample genotypes output by STITCH imputation are represented in the form of genotype dosages. Therefore, users should select software capable of handling genotype dosage data for their GWAS analysis. Although PLINK2 offers a hard-call mode that can convert genotype dosages into the traditional genotypic forms of 0, 1, and 2, this approach leads to the loss of some genetic information, which is not recommended. c. Prepare the phenotype file and covariates for GWAS analysis. CRITICAL: Before starting the GWAS, we need to prepare the phenotype and covariate files. Please refer to phenotype_example.txt and PCA_out_example.txt provided in our GitHub repository to edit the corresponding files. There are two files to edit: 1. Phenotype File: Each row represents a sample, with the first column containing the sample name and each subsequent column representing a different phenotype. For binary variables, 1 and 2 denote different groups, while 0 or NA indicates missing data. 2. Covariate File: This file follows a similar structure to the phenotype file. Each row corresponds to a sample, with the first column for the sample name and each subsequent column for a covariate. In GWAS research, we typically use the first five principal components (PCs) as covariates, along with other factors such as gender and age. Please set the number and type of covariates based on the characteristics of the phenotype relevant to your GWAS analysis. You can edit the ‘.eigenvac’ file generated by 3.3.Plink_load_forBaseVarC.sh by adding extra columns of covariates (e.g., age and BMI). Note: For PLINK2 input formats, please check the official document at: https://www.cog-genomics.org/plink/2.0/input#pheno . d. Perform GWAS analysis. i. Specify the phenotype name and output options in 3.5.2.run_GWAS.sh. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --glm \ --maf 0.05 \ --hwe 1e-6 \ --geno 0.1 dosage \ --pheno ./phenotype_genaral/pheno.csv \ --pheno-name pheno1 \ --pheno-quantile-normalize \ --covar /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5.eigenvec \ --covar-variance-standardize \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/GWAS/STITCH.pheno_1 Note: In the previous step, we consolidated all phenotypes into a single file. Here, you only need to update the content of the ‘ --pheno-name’ parameter to reflect the phenotype used in each round of GWAS analysis. The quality control parameters ( --maf, --hwe, --geno ) have been set using the most widely accepted methods available. Please review your own data and make any necessary adjustments. CRITICAL: Please note that different phenotypes may require distinct covariates, so be sure to adjust them according to your specific phenotype. Note: The final two lines of 3.5.2.run_GWAS.sh process the output file from the GWAS analysis, focusing specifically on retaining the additive effect component. We will ultimately use this result to generate a Manhattan plot. for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.linear`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.logistic.hybrid`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done Note: As illustrated, the GWAS result file contains entries for single loci tested against multiple covariates. It is necessary to retain only the entries labeled 'ADD' in the 'TEST' column ( Figure 2 ), which represent the outcomes post-adjustment for covariates. The final two lines of code in 3.5.run_GWAS.sh are executing this specific task. Figure 2 The ‘ADD’ lines of GWAS output 9. Visualize the GWAS result. a. Use the R package qqman to visualize the GWAS results: library(qqman) gwas_result <- read.table('STITCH.pheno_1.pheno1.glm.linear.add') manhattan(gwas_result, chr = `#CHROM`, bp = 'POS', p = 'P', snp = 'ID') qq(gwas_result$P) Note: In our code package, we provide a compressed file containing the GWAS results ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, you can use the data to perform visualization tests. Conduct segmented single-base mutation detection using BaseVarC and merge the results for principal component analysis (PCA). Note: The primary goal of this step is to perform variant detection across all samples, resulting in the creation of a ‘.vcf’ file that includes data from all samples. However, this ‘.vcf’ file is distinct from the one generated through genotype imputation. Its purpose is to provide a more comprehensive set of original variant sites for PCA analysis before genotype imputation. For subsequent GWAS analyses, the ‘.vcf’ file produced from genotype imputation will serve as the input for population variation. a. Modify the arguments of 3.1.BaseVarC_prep.py. Note: Due to the sample size considerations in STITCH, segmenting the processing can help conserve computing resources and time. Here, we provide a script that divides chromosomes into smaller chunks and generates scripts for BaseVarC. i. Specify the parameters in 3.1.BaseVarC_prep.py. ii. Execute the script to segment the chromosome regions. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/3.Genome-wide_association_study_using_PLINK2' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #sample names and bam paths, need to match each other bamlist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' #Reference Genome, should be the same as 01.workflow_bwa.sh ref_genome='/path_to/Required_reference_data/hg38/hg38.fa.masked.gz' #the path to path of BaseVarC Basevar='/path_to/BaseVarC' Note: It is recommended to ensure consistency with the parameters used previously for each setting here. b. Run the following command: Python 3.1.BaseVarC_prep.py. c. Execute the script generated. Note: In the directory ‘3.Genome-wide_association_study_using_PLINK2/bash’, you will find the relevant script named ‘CMD_basevar_chrx.sh’. You can refer to the methods previously used for running STITCH imputation to execute this script. d. Merge all BaseVarC results using BCFtools. You can use the ls command to locate all completed files: ls /path_to/3.Genome-wide_association_study_using_PLINK2/output/∗/∗.vcf.gz > /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list Alternatives: You can use the VCF list generated by 3.1.BaseVarC_prep.py. e. Specify the arguments in 3.2.Merge.sh, and then run the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list \ -Oz \ -o /path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz f. Modify the arguments of 3.3.Plink_load_forBaseVarC.sh. i. Enter the absolute path for PLINK2. ii. Fill in the corresponding parameters for the ‘.vcf’ file generated after executing 3.2.Merge.sh. iii. Run the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 3.2.Merge.sh merge_out=/path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz $plink2 \ --pca 5 \ --vcf $merge_out \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5 Note: In this step, we utilize the population single-base variation ‘.vcf’ file generated by BaseVarC for PCA calculations, rather than the ‘.vcf’ result file produced by STITCH. Note: In this step, we did not convert the ‘.vcf’ file to ‘.pgen’ format using PLINK2 before performing PCA calculations; therefore, rerunning this script will not save time. If you need to use this ‘.vcf’ data multiple times for PLINK2 operations, you can add the ‘--make-pgen’ option in the script and then replace the ‘--vcf’ parameter with ‘--pfile’ when recalculating PCA. Alternatives: In this protocol, the computations performed by PLINK2 are limited to 5 principal components (PCs). If the analysis involves a larger sample size or a more complex population structure, it is advisable to adjust the parameters for PCA. You can change the number of PCs calculated by using the ‘--pca’ option. Note: The primary goal of this step is to perform variant detection across all samples, resulting in the creation of a ‘.vcf’ file that includes data from all samples. However, this ‘.vcf’ file is distinct from the one generated through genotype imputation. Its purpose is to provide a more comprehensive set of original variant sites for PCA analysis before genotype imputation. For subsequent GWAS analyses, the ‘.vcf’ file produced from genotype imputation will serve as the input for population variation. Modify the arguments of 3.1.BaseVarC_prep.py. Note: Due to the sample size considerations in STITCH, segmenting the processing can help conserve computing resources and time. Here, we provide a script that divides chromosomes into smaller chunks and generates scripts for BaseVarC. i. Specify the parameters in 3.1.BaseVarC_prep.py. ii. Execute the script to segment the chromosome regions. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/3.Genome-wide_association_study_using_PLINK2' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #sample names and bam paths, need to match each other bamlist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' #Reference Genome, should be the same as 01.workflow_bwa.sh ref_genome='/path_to/Required_reference_data/hg38/hg38.fa.masked.gz' #the path to path of BaseVarC Basevar='/path_to/BaseVarC' Note: It is recommended to ensure consistency with the parameters used previously for each setting here. Note: Due to the sample size considerations in STITCH, segmenting the processing can help conserve computing resources and time. Here, we provide a script that divides chromosomes into smaller chunks and generates scripts for BaseVarC. Specify the parameters in 3.1.BaseVarC_prep.py. Execute the script to segment the chromosome regions. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/3.Genome-wide_association_study_using_PLINK2' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #sample names and bam paths, need to match each other bamlist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' #Reference Genome, should be the same as 01.workflow_bwa.sh ref_genome='/path_to/Required_reference_data/hg38/hg38.fa.masked.gz' #the path to path of BaseVarC Basevar='/path_to/BaseVarC' Note: It is recommended to ensure consistency with the parameters used previously for each setting here. #Specify the arguments before excute this code #this is the chunk size, the smaller the size, lower compute resoure required chunk_size = int('5000000') #the folder of working outdir='/path_to/3.Genome-wide_association_study_using_PLINK2' #the length of each chromosomes, example file was in GRCh38, please change if running on other versions of reference genome CHR_length='/path_to/2.Genotype_imputation_using_STITCH/CHR.len' #sample names and bam paths, need to match each other bamlist='/path_to/2.Genotype_imputation_using_STITCH/bampath.list' #Reference Genome, should be the same as 01.workflow_bwa.sh ref_genome='/path_to/Required_reference_data/hg38/hg38.fa.masked.gz' #the path to path of BaseVarC Basevar='/path_to/BaseVarC' Note: It is recommended to ensure consistency with the parameters used previously for each setting here. Run the following command: Python 3.1.BaseVarC_prep.py. Execute the script generated. Note: In the directory ‘3.Genome-wide_association_study_using_PLINK2/bash’, you will find the relevant script named ‘CMD_basevar_chrx.sh’. You can refer to the methods previously used for running STITCH imputation to execute this script. Note: In the directory ‘3.Genome-wide_association_study_using_PLINK2/bash’, you will find the relevant script named ‘CMD_basevar_chrx.sh’. You can refer to the methods previously used for running STITCH imputation to execute this script. Merge all BaseVarC results using BCFtools. You can use the ls command to locate all completed files: ls /path_to/3.Genome-wide_association_study_using_PLINK2/output/∗/∗.vcf.gz > /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list Alternatives: You can use the VCF list generated by 3.1.BaseVarC_prep.py. ls /path_to/3.Genome-wide_association_study_using_PLINK2/output/∗/∗.vcf.gz > /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list Alternatives: You can use the VCF list generated by 3.1.BaseVarC_prep.py. Specify the arguments in 3.2.Merge.sh, and then run the code. #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list \ -Oz \ -o /path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz #Specify the arguments before excute the code related #path for bcftools bcftools=/path_to/bcftools $bcftools concat \ -f /path_to/3.Genome-wide_association_study_using_PLINK2/completed_BaseVarC_vcf.list \ -Oz \ -o /path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz Modify the arguments of 3.3.Plink_load_forBaseVarC.sh. i. Enter the absolute path for PLINK2. ii. Fill in the corresponding parameters for the ‘.vcf’ file generated after executing 3.2.Merge.sh. iii. Run the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 3.2.Merge.sh merge_out=/path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz $plink2 \ --pca 5 \ --vcf $merge_out \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5 Note: In this step, we utilize the population single-base variation ‘.vcf’ file generated by BaseVarC for PCA calculations, rather than the ‘.vcf’ result file produced by STITCH. Note: In this step, we did not convert the ‘.vcf’ file to ‘.pgen’ format using PLINK2 before performing PCA calculations; therefore, rerunning this script will not save time. If you need to use this ‘.vcf’ data multiple times for PLINK2 operations, you can add the ‘--make-pgen’ option in the script and then replace the ‘--vcf’ parameter with ‘--pfile’ when recalculating PCA. Alternatives: In this protocol, the computations performed by PLINK2 are limited to 5 principal components (PCs). If the analysis involves a larger sample size or a more complex population structure, it is advisable to adjust the parameters for PCA. You can change the number of PCs calculated by using the ‘--pca’ option. Enter the absolute path for PLINK2. Fill in the corresponding parameters for the ‘.vcf’ file generated after executing 3.2.Merge.sh. Run the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 3.2.Merge.sh merge_out=/path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz $plink2 \ --pca 5 \ --vcf $merge_out \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5 Note: In this step, we utilize the population single-base variation ‘.vcf’ file generated by BaseVarC for PCA calculations, rather than the ‘.vcf’ result file produced by STITCH. Note: In this step, we did not convert the ‘.vcf’ file to ‘.pgen’ format using PLINK2 before performing PCA calculations; therefore, rerunning this script will not save time. If you need to use this ‘.vcf’ data multiple times for PLINK2 operations, you can add the ‘--make-pgen’ option in the script and then replace the ‘--vcf’ parameter with ‘--pfile’ when recalculating PCA. Alternatives: In this protocol, the computations performed by PLINK2 are limited to 5 principal components (PCs). If the analysis involves a larger sample size or a more complex population structure, it is advisable to adjust the parameters for PCA. You can change the number of PCs calculated by using the ‘--pca’ option. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 3.2.Merge.sh merge_out=/path_to/3.Genome-wide_association_study_using_PLINK2/final/BaseVarC.vcf.gz $plink2 \ --pca 5 \ --vcf $merge_out \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5 Note: In this step, we utilize the population single-base variation ‘.vcf’ file generated by BaseVarC for PCA calculations, rather than the ‘.vcf’ result file produced by STITCH. Note: In this step, we did not convert the ‘.vcf’ file to ‘.pgen’ format using PLINK2 before performing PCA calculations; therefore, rerunning this script will not save time. If you need to use this ‘.vcf’ data multiple times for PLINK2 operations, you can add the ‘--make-pgen’ option in the script and then replace the ‘--vcf’ parameter with ‘--pfile’ when recalculating PCA. Alternatives: In this protocol, the computations performed by PLINK2 are limited to 5 principal components (PCs). If the analysis involves a larger sample size or a more complex population structure, it is advisable to adjust the parameters for PCA. You can change the number of PCs calculated by using the ‘--pca’ option. Perform GWAS using PLINK2. Note: In this step, we will conduct GWAS analysis using the PCA results obtained in the previous step as covariates, ultimately generating the data required to create the Manhattan plot. a. Load the pfile used in GWAS using PLINK2. i. Enter the corresponding parameters for the ‘.vcf’ file generated after completing 2.3.Merge.sh, along with the absolute path for PLINK2. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 2.3.Merge.sh merge_out=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz $plink2 \ --freq \ --hardy \ --make-pgen 'vzs' \ --set-all-var-ids chr@:# \ --vcf $merge_out dosage=DS \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH Note: In this step, we convert the ‘.vcf’ file generated by STITCH into a 'pgen’ file compatible with PLINK for GWAS. The generated ‘pgen’ file will serve as the genotype input for subsequent GWAS analyses. The ‘--make-pgen’ argument includes the ‘vzs’ parameter, which compresses the output genotype file to save space, making it particularly useful for large sample GWAS. b. Calculate the kinship among samples. i. Specify the arguments in ‘3.5.1.kinship_analysis.sh’ and execute the code: Note: In GWAS, statistical analyses based on general linear models require that samples be independent to ensure the reliability of the results. Samples with kinship relationships violate this assumption of independence and should be excluded from the dataset. In this step, we use PLINK2 to calculate the kinship among samples, enabling the exclusion of related samples in subsequent analyses. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --make-king triangle bin\ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship Note: The result of the ‘--make-king’ calculation is a kinship matrix file with the suffix ‘.king.bin’, which contains the pairwise kinship relationships between each sample. The kinship coefficient for duplicate samples and identical twins is 0.5. For first-degree relationships, such as parent-child or full siblings, the kinship coefficient is 0.25, and for second-degree relationships, it is 0.125, and so on. Alternatives: To exclude samples with first-degree relationships, add the ‘--king-cutoff /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship 0.177’ parameter when loading the pfile. Alternatively, you can use your preferred algorithm to identify the samples that need to be excluded based on the kinship matrix, and then remove these samples from the dataset using the ‘--remove’ flag. Note: 0.177 is the geometric mean of 0.25 and 0.125. When selecting a cutoff threshold, the geometric mean between two degrees of relationship is commonly used. Alternatives: In the context of GWAS, it's essential to note that the input samples for the generalized linear model must be independent of each other. However, when utilizing software designed for mixed linear models, such as GEMMA, 20 it is possible to incorporate a kinship matrix as a parameter for GWAS analysis. Interested readers are encouraged to conduct their analyses using relevant software. It is important to remember that the sample genotypes output by STITCH imputation are represented in the form of genotype dosages. Therefore, users should select software capable of handling genotype dosage data for their GWAS analysis. Although PLINK2 offers a hard-call mode that can convert genotype dosages into the traditional genotypic forms of 0, 1, and 2, this approach leads to the loss of some genetic information, which is not recommended. c. Prepare the phenotype file and covariates for GWAS analysis. CRITICAL: Before starting the GWAS, we need to prepare the phenotype and covariate files. Please refer to phenotype_example.txt and PCA_out_example.txt provided in our GitHub repository to edit the corresponding files. There are two files to edit: 1. Phenotype File: Each row represents a sample, with the first column containing the sample name and each subsequent column representing a different phenotype. For binary variables, 1 and 2 denote different groups, while 0 or NA indicates missing data. 2. Covariate File: This file follows a similar structure to the phenotype file. Each row corresponds to a sample, with the first column for the sample name and each subsequent column for a covariate. In GWAS research, we typically use the first five principal components (PCs) as covariates, along with other factors such as gender and age. Please set the number and type of covariates based on the characteristics of the phenotype relevant to your GWAS analysis. You can edit the ‘.eigenvac’ file generated by 3.3.Plink_load_forBaseVarC.sh by adding extra columns of covariates (e.g., age and BMI). Note: For PLINK2 input formats, please check the official document at: https://www.cog-genomics.org/plink/2.0/input#pheno . d. Perform GWAS analysis. i. Specify the phenotype name and output options in 3.5.2.run_GWAS.sh. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --glm \ --maf 0.05 \ --hwe 1e-6 \ --geno 0.1 dosage \ --pheno ./phenotype_genaral/pheno.csv \ --pheno-name pheno1 \ --pheno-quantile-normalize \ --covar /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5.eigenvec \ --covar-variance-standardize \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/GWAS/STITCH.pheno_1 Note: In the previous step, we consolidated all phenotypes into a single file. Here, you only need to update the content of the ‘ --pheno-name’ parameter to reflect the phenotype used in each round of GWAS analysis. The quality control parameters ( --maf, --hwe, --geno ) have been set using the most widely accepted methods available. Please review your own data and make any necessary adjustments. CRITICAL: Please note that different phenotypes may require distinct covariates, so be sure to adjust them according to your specific phenotype. Note: The final two lines of 3.5.2.run_GWAS.sh process the output file from the GWAS analysis, focusing specifically on retaining the additive effect component. We will ultimately use this result to generate a Manhattan plot. for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.linear`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.logistic.hybrid`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done Note: As illustrated, the GWAS result file contains entries for single loci tested against multiple covariates. It is necessary to retain only the entries labeled 'ADD' in the 'TEST' column ( Figure 2 ), which represent the outcomes post-adjustment for covariates. The final two lines of code in 3.5.run_GWAS.sh are executing this specific task. Figure 2 The ‘ADD’ lines of GWAS output Note: In this step, we will conduct GWAS analysis using the PCA results obtained in the previous step as covariates, ultimately generating the data required to create the Manhattan plot. Load the pfile used in GWAS using PLINK2. i. Enter the corresponding parameters for the ‘.vcf’ file generated after completing 2.3.Merge.sh, along with the absolute path for PLINK2. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 2.3.Merge.sh merge_out=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz $plink2 \ --freq \ --hardy \ --make-pgen 'vzs' \ --set-all-var-ids chr@:# \ --vcf $merge_out dosage=DS \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH Note: In this step, we convert the ‘.vcf’ file generated by STITCH into a 'pgen’ file compatible with PLINK for GWAS. The generated ‘pgen’ file will serve as the genotype input for subsequent GWAS analyses. The ‘--make-pgen’ argument includes the ‘vzs’ parameter, which compresses the output genotype file to save space, making it particularly useful for large sample GWAS. Enter the corresponding parameters for the ‘.vcf’ file generated after completing 2.3.Merge.sh, along with the absolute path for PLINK2. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 2.3.Merge.sh merge_out=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz $plink2 \ --freq \ --hardy \ --make-pgen 'vzs' \ --set-all-var-ids chr@:# \ --vcf $merge_out dosage=DS \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH Note: In this step, we convert the ‘.vcf’ file generated by STITCH into a 'pgen’ file compatible with PLINK for GWAS. The generated ‘pgen’ file will serve as the genotype input for subsequent GWAS analyses. The ‘--make-pgen’ argument includes the ‘vzs’ parameter, which compresses the output genotype file to save space, making it particularly useful for large sample GWAS. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 #the output vcf file of 2.3.Merge.sh merge_out=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz $plink2 \ --freq \ --hardy \ --make-pgen 'vzs' \ --set-all-var-ids chr@:# \ --vcf $merge_out dosage=DS \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH Note: In this step, we convert the ‘.vcf’ file generated by STITCH into a 'pgen’ file compatible with PLINK for GWAS. The generated ‘pgen’ file will serve as the genotype input for subsequent GWAS analyses. The ‘--make-pgen’ argument includes the ‘vzs’ parameter, which compresses the output genotype file to save space, making it particularly useful for large sample GWAS. Calculate the kinship among samples. i. Specify the arguments in ‘3.5.1.kinship_analysis.sh’ and execute the code: Note: In GWAS, statistical analyses based on general linear models require that samples be independent to ensure the reliability of the results. Samples with kinship relationships violate this assumption of independence and should be excluded from the dataset. In this step, we use PLINK2 to calculate the kinship among samples, enabling the exclusion of related samples in subsequent analyses. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --make-king triangle bin\ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship Note: The result of the ‘--make-king’ calculation is a kinship matrix file with the suffix ‘.king.bin’, which contains the pairwise kinship relationships between each sample. The kinship coefficient for duplicate samples and identical twins is 0.5. For first-degree relationships, such as parent-child or full siblings, the kinship coefficient is 0.25, and for second-degree relationships, it is 0.125, and so on. Alternatives: To exclude samples with first-degree relationships, add the ‘--king-cutoff /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship 0.177’ parameter when loading the pfile. Alternatively, you can use your preferred algorithm to identify the samples that need to be excluded based on the kinship matrix, and then remove these samples from the dataset using the ‘--remove’ flag. Note: 0.177 is the geometric mean of 0.25 and 0.125. When selecting a cutoff threshold, the geometric mean between two degrees of relationship is commonly used. Alternatives: In the context of GWAS, it's essential to note that the input samples for the generalized linear model must be independent of each other. However, when utilizing software designed for mixed linear models, such as GEMMA, 20 it is possible to incorporate a kinship matrix as a parameter for GWAS analysis. Interested readers are encouraged to conduct their analyses using relevant software. It is important to remember that the sample genotypes output by STITCH imputation are represented in the form of genotype dosages. Therefore, users should select software capable of handling genotype dosage data for their GWAS analysis. Although PLINK2 offers a hard-call mode that can convert genotype dosages into the traditional genotypic forms of 0, 1, and 2, this approach leads to the loss of some genetic information, which is not recommended. Specify the arguments in ‘3.5.1.kinship_analysis.sh’ and execute the code: Note: In GWAS, statistical analyses based on general linear models require that samples be independent to ensure the reliability of the results. Samples with kinship relationships violate this assumption of independence and should be excluded from the dataset. In this step, we use PLINK2 to calculate the kinship among samples, enabling the exclusion of related samples in subsequent analyses. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --make-king triangle bin\ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship Note: The result of the ‘--make-king’ calculation is a kinship matrix file with the suffix ‘.king.bin’, which contains the pairwise kinship relationships between each sample. The kinship coefficient for duplicate samples and identical twins is 0.5. For first-degree relationships, such as parent-child or full siblings, the kinship coefficient is 0.25, and for second-degree relationships, it is 0.125, and so on. Alternatives: To exclude samples with first-degree relationships, add the ‘--king-cutoff /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship 0.177’ parameter when loading the pfile. Alternatively, you can use your preferred algorithm to identify the samples that need to be excluded based on the kinship matrix, and then remove these samples from the dataset using the ‘--remove’ flag. Note: 0.177 is the geometric mean of 0.25 and 0.125. When selecting a cutoff threshold, the geometric mean between two degrees of relationship is commonly used. Alternatives: In the context of GWAS, it's essential to note that the input samples for the generalized linear model must be independent of each other. However, when utilizing software designed for mixed linear models, such as GEMMA, 20 it is possible to incorporate a kinship matrix as a parameter for GWAS analysis. Interested readers are encouraged to conduct their analyses using relevant software. It is important to remember that the sample genotypes output by STITCH imputation are represented in the form of genotype dosages. Therefore, users should select software capable of handling genotype dosage data for their GWAS analysis. Although PLINK2 offers a hard-call mode that can convert genotype dosages into the traditional genotypic forms of 0, 1, and 2, this approach leads to the loss of some genetic information, which is not recommended. Note: In GWAS, statistical analyses based on general linear models require that samples be independent to ensure the reliability of the results. Samples with kinship relationships violate this assumption of independence and should be excluded from the dataset. In this step, we use PLINK2 to calculate the kinship among samples, enabling the exclusion of related samples in subsequent analyses. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --make-king triangle bin\ --out /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship Note: The result of the ‘--make-king’ calculation is a kinship matrix file with the suffix ‘.king.bin’, which contains the pairwise kinship relationships between each sample. The kinship coefficient for duplicate samples and identical twins is 0.5. For first-degree relationships, such as parent-child or full siblings, the kinship coefficient is 0.25, and for second-degree relationships, it is 0.125, and so on. Alternatives: To exclude samples with first-degree relationships, add the ‘--king-cutoff /path_to/3.Genome-wide_association_study_using_PLINK2/final/kinship 0.177’ parameter when loading the pfile. Alternatively, you can use your preferred algorithm to identify the samples that need to be excluded based on the kinship matrix, and then remove these samples from the dataset using the ‘--remove’ flag. Note: 0.177 is the geometric mean of 0.25 and 0.125. When selecting a cutoff threshold, the geometric mean between two degrees of relationship is commonly used. Alternatives: In the context of GWAS, it's essential to note that the input samples for the generalized linear model must be independent of each other. However, when utilizing software designed for mixed linear models, such as GEMMA, 20 it is possible to incorporate a kinship matrix as a parameter for GWAS analysis. Interested readers are encouraged to conduct their analyses using relevant software. It is important to remember that the sample genotypes output by STITCH imputation are represented in the form of genotype dosages. Therefore, users should select software capable of handling genotype dosage data for their GWAS analysis. Although PLINK2 offers a hard-call mode that can convert genotype dosages into the traditional genotypic forms of 0, 1, and 2, this approach leads to the loss of some genetic information, which is not recommended. Prepare the phenotype file and covariates for GWAS analysis. CRITICAL: Before starting the GWAS, we need to prepare the phenotype and covariate files. Please refer to phenotype_example.txt and PCA_out_example.txt provided in our GitHub repository to edit the corresponding files. There are two files to edit: 1. Phenotype File: Each row represents a sample, with the first column containing the sample name and each subsequent column representing a different phenotype. For binary variables, 1 and 2 denote different groups, while 0 or NA indicates missing data. 2. Covariate File: This file follows a similar structure to the phenotype file. Each row corresponds to a sample, with the first column for the sample name and each subsequent column for a covariate. In GWAS research, we typically use the first five principal components (PCs) as covariates, along with other factors such as gender and age. Please set the number and type of covariates based on the characteristics of the phenotype relevant to your GWAS analysis. You can edit the ‘.eigenvac’ file generated by 3.3.Plink_load_forBaseVarC.sh by adding extra columns of covariates (e.g., age and BMI). Note: For PLINK2 input formats, please check the official document at: https://www.cog-genomics.org/plink/2.0/input#pheno . CRITICAL: Before starting the GWAS, we need to prepare the phenotype and covariate files. Please refer to phenotype_example.txt and PCA_out_example.txt provided in our GitHub repository to edit the corresponding files. There are two files to edit: 1. Phenotype File: Each row represents a sample, with the first column containing the sample name and each subsequent column representing a different phenotype. For binary variables, 1 and 2 denote different groups, while 0 or NA indicates missing data. 2. Covariate File: This file follows a similar structure to the phenotype file. Each row corresponds to a sample, with the first column for the sample name and each subsequent column for a covariate. In GWAS research, we typically use the first five principal components (PCs) as covariates, along with other factors such as gender and age. Please set the number and type of covariates based on the characteristics of the phenotype relevant to your GWAS analysis. You can edit the ‘.eigenvac’ file generated by 3.3.Plink_load_forBaseVarC.sh by adding extra columns of covariates (e.g., age and BMI). Note: For PLINK2 input formats, please check the official document at: https://www.cog-genomics.org/plink/2.0/input#pheno . Perform GWAS analysis. i. Specify the phenotype name and output options in 3.5.2.run_GWAS.sh. ii. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --glm \ --maf 0.05 \ --hwe 1e-6 \ --geno 0.1 dosage \ --pheno ./phenotype_genaral/pheno.csv \ --pheno-name pheno1 \ --pheno-quantile-normalize \ --covar /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5.eigenvec \ --covar-variance-standardize \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/GWAS/STITCH.pheno_1 Note: In the previous step, we consolidated all phenotypes into a single file. Here, you only need to update the content of the ‘ --pheno-name’ parameter to reflect the phenotype used in each round of GWAS analysis. The quality control parameters ( --maf, --hwe, --geno ) have been set using the most widely accepted methods available. Please review your own data and make any necessary adjustments. CRITICAL: Please note that different phenotypes may require distinct covariates, so be sure to adjust them according to your specific phenotype. Note: The final two lines of 3.5.2.run_GWAS.sh process the output file from the GWAS analysis, focusing specifically on retaining the additive effect component. We will ultimately use this result to generate a Manhattan plot. for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.linear`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.logistic.hybrid`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done Note: As illustrated, the GWAS result file contains entries for single loci tested against multiple covariates. It is necessary to retain only the entries labeled 'ADD' in the 'TEST' column ( Figure 2 ), which represent the outcomes post-adjustment for covariates. The final two lines of code in 3.5.run_GWAS.sh are executing this specific task. Figure 2 The ‘ADD’ lines of GWAS output Specify the phenotype name and output options in 3.5.2.run_GWAS.sh. Execute the script. #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --glm \ --maf 0.05 \ --hwe 1e-6 \ --geno 0.1 dosage \ --pheno ./phenotype_genaral/pheno.csv \ --pheno-name pheno1 \ --pheno-quantile-normalize \ --covar /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5.eigenvec \ --covar-variance-standardize \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/GWAS/STITCH.pheno_1 Note: In the previous step, we consolidated all phenotypes into a single file. Here, you only need to update the content of the ‘ --pheno-name’ parameter to reflect the phenotype used in each round of GWAS analysis. The quality control parameters ( --maf, --hwe, --geno ) have been set using the most widely accepted methods available. Please review your own data and make any necessary adjustments. CRITICAL: Please note that different phenotypes may require distinct covariates, so be sure to adjust them according to your specific phenotype. Note: The final two lines of 3.5.2.run_GWAS.sh process the output file from the GWAS analysis, focusing specifically on retaining the additive effect component. We will ultimately use this result to generate a Manhattan plot. for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.linear`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.logistic.hybrid`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done Note: As illustrated, the GWAS result file contains entries for single loci tested against multiple covariates. It is necessary to retain only the entries labeled 'ADD' in the 'TEST' column ( Figure 2 ), which represent the outcomes post-adjustment for covariates. The final two lines of code in 3.5.run_GWAS.sh are executing this specific task. Figure 2 The ‘ADD’ lines of GWAS output #Specify the arguments before execute the code #path for plink2 plink2=/path_to/plink2 $plink2 \ --pfile /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH 'vzs' \ --read-freq /path_to/3.Genome-wide_association_study_using_PLINK2/final/STITCH.afreq \ --glm \ --maf 0.05 \ --hwe 1e-6 \ --geno 0.1 dosage \ --pheno ./phenotype_genaral/pheno.csv \ --pheno-name pheno1 \ --pheno-quantile-normalize \ --covar /path_to/3.Genome-wide_association_study_using_PLINK2/final/GWAS_pca_5.eigenvec \ --covar-variance-standardize \ --out /path_to/3.Genome-wide_association_study_using_PLINK2/GWAS/STITCH.pheno_1 Note: In the previous step, we consolidated all phenotypes into a single file. Here, you only need to update the content of the ‘ --pheno-name’ parameter to reflect the phenotype used in each round of GWAS analysis. The quality control parameters ( --maf, --hwe, --geno ) have been set using the most widely accepted methods available. Please review your own data and make any necessary adjustments. CRITICAL: Please note that different phenotypes may require distinct covariates, so be sure to adjust them according to your specific phenotype. Note: The final two lines of 3.5.2.run_GWAS.sh process the output file from the GWAS analysis, focusing specifically on retaining the additive effect component. We will ultimately use this result to generate a Manhattan plot. for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.linear`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done for i in `ls /path_to/3.Genome-wide_association_study_using_PLINK2/final/∗.logistic.hybrid`; do head $i -n1 >> $i.add && grep 'ADD' $i >> $i.add ; done Note: As illustrated, the GWAS result file contains entries for single loci tested against multiple covariates. It is necessary to retain only the entries labeled 'ADD' in the 'TEST' column ( Figure 2 ), which represent the outcomes post-adjustment for covariates. The final two lines of code in 3.5.run_GWAS.sh are executing this specific task. The ‘ADD’ lines of GWAS output Visualize the GWAS result. a. Use the R package qqman to visualize the GWAS results: library(qqman) gwas_result <- read.table('STITCH.pheno_1.pheno1.glm.linear.add') manhattan(gwas_result, chr = `#CHROM`, bp = 'POS', p = 'P', snp = 'ID') qq(gwas_result$P) Note: In our code package, we provide a compressed file containing the GWAS results ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, you can use the data to perform visualization tests. Use the R package qqman to visualize the GWAS results: library(qqman) gwas_result <- read.table('STITCH.pheno_1.pheno1.glm.linear.add') manhattan(gwas_result, chr = `#CHROM`, bp = 'POS', p = 'P', snp = 'ID') qq(gwas_result$P) Note: In our code package, we provide a compressed file containing the GWAS results ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, you can use the data to perform visualization tests. library(qqman) gwas_result <- read.table('STITCH.pheno_1.pheno1.glm.linear.add') manhattan(gwas_result, chr = `#CHROM`, bp = 'POS', p = 'P', snp = 'ID') qq(gwas_result$P) Note: In our code package, we provide a compressed file containing the GWAS results ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, you can use the data to perform visualization tests. Timing: 1–2 days Timing: 1–2 days This section primarily focuses on conducting heritability and genetic correlation analyses on various phenotypes using the summary statistics generated from the previous step’s GWAS. The process begins by converting the summary statistics into a format compatible with the software LDSC. Subsequently, LDSC is utilized to perform heritability analyses. Upon obtaining these results, the heritability outcomes of the phenotypes of interest are paired and input into LDSC for genetic correlation analysis. The final output includes data on the heritability of the phenotypes and the genetic correlation between pairs of phenotypes. Note: This analysis can be performed using the sample data provided in our code package, located at ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, the data will be ready for use. 10. Reference file preparation. Note: To start, it is crucial to convert the chromosomal positions (chr and pos) into their corresponding rsIDs. This process requires obtaining data files that map genetic variations to rsIDs across various chromosomes. These files can be sourced from the downloads provided by dbSNP and are listed in the key resources table for easy access. a. Adjust chrpos_2_rs.py. i. Specify the input and output parameters of chrpos_2_rs.py, which are outlined at the beginning of the code: # directory to the download files input_path = '/path_to/Required_reference_data/dbSNP/' outfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt' CRITICAL: Please note that this code will consolidate all downloaded dbSNP files (content from various chromosomes) into a single output file. To streamline this process, it is recommended to place all downloaded chromosome dbSNP files in a single directory beforehand. b. Update the parameters in the script 4.1.File_preparations.sh. i. Specifying the path to the PLINK2 software and the path to the genotype imputation results file generated in the “ genotype imputation using STITCH ” step. Python=/path_to/Python3 #plink afreq and chr_pos reference files #path for plink2 plink2=/path_to/plink2 #the input vcf file produced by 'Genotype imputation using STITCH' stitch=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz freq_out=/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out ii. Once these updates are made, execute the 4.1.File_preparations.sh script to complete the file preparation process: sh 4.1.File_preparations.sh CRITICAL: Please note that this step requires Python version 3. Be sure to avoid confusing the two versions. 11. Data conversion. a. Adjust glm_2_ss.py. i. Update the reference files specified at the beginning of the glm_2_ss.py script. Ensure that the primary inputs for this stage, derived from the previous ‘Reference file preparations’ step, are correctly assigned their respective file paths before moving on to the next phase. #outputs of the previous step index = pd.read_csv('/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt',sep='\t') freqfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out.afreq' b. Modify the 4.2.Conversion.sh script: #Python version3 Python=/path_to/Python3 #The python codes provided in github glm_2_ss=/path_to/4.Heritability_and_genetic_correlation_analysis/glm_2_ss.py ss_2_ldsc=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_ldsc.py ss_2_snplist=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_snplist.py c. Run the code, the external parameters include GWAS result data, phenotype names, and the designated output directory: sh 4.2.Conversion.sh /path_to/pheno1.add pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out CRITICAL: Please note that this step requires Python version 3. Be sure not to confuse the two versions. Note: Each GWAS result needs this transformation to calculate heritability and genetic correlation. We recommend converting all GWAS summary statistics in this step to facilitate subsequent genetic correlation analyses. Note: In this step, the output files will have the extensions ‘.ss’, ‘.snplist’, and ‘.ldsc’. The files with ‘.snplist’ and ‘.ldsc’ extensions will be used in subsequent steps. It is acceptable to set the same output folder for all phenotype files, as long as each phenotype name remains unique. 12. Heritability analysis. a. Specify the arguments in 4.3.Heritability_analysis.sh. i. Input the previously obtained LD reference files for the parameters ‘--ref-ld-chr’ and ‘--w-ld-chr’ : #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation munge_sumstats=/path_to/ldsc-master/munge_sumstats.py ldsc=/path_to/ldsc-master/ldsc.py ii. After that, run LDSC with: sh 4.3.Heritability_analysis.sh pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out /path_to/4.Heritability_and_genetic_correlation_analysis/heritability CRITICAL: Please note that all Python instances within the script 4.3.Heritability_analysis.sh utilize Python 2. It is essential to adjust the version accordingly and ensure that the correct path for the LDSC software functions is specified. 13. Genetic correlation analysis. a. Specify the arguments in 4.4.Genetic_correlation.sh: #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py b. Then run the code: sh 4.4.Genetic_correlation.sh pheno_name1 pheno_name2 /path_to/4.Heritability_and_genetic_correlation_analysis/heritability /path_to/4.Heritability_and_genetic_correlation_analysis/correlation CRITICAL: Please note that this step requires Python version 2. Make sure not to confuse it with Python version 3. Note: This analysis can be performed using the sample data provided in our code package, located at ‘/path_to/Supporting_data/example_GWAS_outcome/GWAS_quantitative_20.P_LCR.glm.linear.add.gz’. After decompressing the file, the data will be ready for use. Reference file preparation. Note: To start, it is crucial to convert the chromosomal positions (chr and pos) into their corresponding rsIDs. This process requires obtaining data files that map genetic variations to rsIDs across various chromosomes. These files can be sourced from the downloads provided by dbSNP and are listed in the key resources table for easy access. a. Adjust chrpos_2_rs.py. i. Specify the input and output parameters of chrpos_2_rs.py, which are outlined at the beginning of the code: # directory to the download files input_path = '/path_to/Required_reference_data/dbSNP/' outfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt' CRITICAL: Please note that this code will consolidate all downloaded dbSNP files (content from various chromosomes) into a single output file. To streamline this process, it is recommended to place all downloaded chromosome dbSNP files in a single directory beforehand. b. Update the parameters in the script 4.1.File_preparations.sh. i. Specifying the path to the PLINK2 software and the path to the genotype imputation results file generated in the “ genotype imputation using STITCH ” step. Python=/path_to/Python3 #plink afreq and chr_pos reference files #path for plink2 plink2=/path_to/plink2 #the input vcf file produced by 'Genotype imputation using STITCH' stitch=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz freq_out=/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out ii. Once these updates are made, execute the 4.1.File_preparations.sh script to complete the file preparation process: sh 4.1.File_preparations.sh CRITICAL: Please note that this step requires Python version 3. Be sure to avoid confusing the two versions. Note: To start, it is crucial to convert the chromosomal positions (chr and pos) into their corresponding rsIDs. This process requires obtaining data files that map genetic variations to rsIDs across various chromosomes. These files can be sourced from the downloads provided by dbSNP and are listed in the key resources table for easy access. Adjust chrpos_2_rs.py. i. Specify the input and output parameters of chrpos_2_rs.py, which are outlined at the beginning of the code: # directory to the download files input_path = '/path_to/Required_reference_data/dbSNP/' outfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt' CRITICAL: Please note that this code will consolidate all downloaded dbSNP files (content from various chromosomes) into a single output file. To streamline this process, it is recommended to place all downloaded chromosome dbSNP files in a single directory beforehand. Specify the input and output parameters of chrpos_2_rs.py, which are outlined at the beginning of the code: # directory to the download files input_path = '/path_to/Required_reference_data/dbSNP/' outfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt' CRITICAL: Please note that this code will consolidate all downloaded dbSNP files (content from various chromosomes) into a single output file. To streamline this process, it is recommended to place all downloaded chromosome dbSNP files in a single directory beforehand. # directory to the download files input_path = '/path_to/Required_reference_data/dbSNP/' outfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt' CRITICAL: Please note that this code will consolidate all downloaded dbSNP files (content from various chromosomes) into a single output file. To streamline this process, it is recommended to place all downloaded chromosome dbSNP files in a single directory beforehand. Update the parameters in the script 4.1.File_preparations.sh. i. Specifying the path to the PLINK2 software and the path to the genotype imputation results file generated in the “ genotype imputation using STITCH ” step. Python=/path_to/Python3 #plink afreq and chr_pos reference files #path for plink2 plink2=/path_to/plink2 #the input vcf file produced by 'Genotype imputation using STITCH' stitch=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz freq_out=/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out ii. Once these updates are made, execute the 4.1.File_preparations.sh script to complete the file preparation process: sh 4.1.File_preparations.sh CRITICAL: Please note that this step requires Python version 3. Be sure to avoid confusing the two versions. Specifying the path to the PLINK2 software and the path to the genotype imputation results file generated in the “ genotype imputation using STITCH ” step. Python=/path_to/Python3 #plink afreq and chr_pos reference files #path for plink2 plink2=/path_to/plink2 #the input vcf file produced by 'Genotype imputation using STITCH' stitch=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz freq_out=/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out Python=/path_to/Python3 #plink afreq and chr_pos reference files #path for plink2 plink2=/path_to/plink2 #the input vcf file produced by 'Genotype imputation using STITCH' stitch=/path_to/2.Genotype_imputation_using_STITCH/final/STITCH.vcf.gz freq_out=/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out Once these updates are made, execute the 4.1.File_preparations.sh script to complete the file preparation process: sh 4.1.File_preparations.sh CRITICAL: Please note that this step requires Python version 3. Be sure to avoid confusing the two versions. sh 4.1.File_preparations.sh CRITICAL: Please note that this step requires Python version 3. Be sure to avoid confusing the two versions. Data conversion. a. Adjust glm_2_ss.py. i. Update the reference files specified at the beginning of the glm_2_ss.py script. Ensure that the primary inputs for this stage, derived from the previous ‘Reference file preparations’ step, are correctly assigned their respective file paths before moving on to the next phase. #outputs of the previous step index = pd.read_csv('/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt',sep='\t') freqfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out.afreq' b. Modify the 4.2.Conversion.sh script: #Python version3 Python=/path_to/Python3 #The python codes provided in github glm_2_ss=/path_to/4.Heritability_and_genetic_correlation_analysis/glm_2_ss.py ss_2_ldsc=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_ldsc.py ss_2_snplist=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_snplist.py c. Run the code, the external parameters include GWAS result data, phenotype names, and the designated output directory: sh 4.2.Conversion.sh /path_to/pheno1.add pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out CRITICAL: Please note that this step requires Python version 3. Be sure not to confuse the two versions. Note: Each GWAS result needs this transformation to calculate heritability and genetic correlation. We recommend converting all GWAS summary statistics in this step to facilitate subsequent genetic correlation analyses. Note: In this step, the output files will have the extensions ‘.ss’, ‘.snplist’, and ‘.ldsc’. The files with ‘.snplist’ and ‘.ldsc’ extensions will be used in subsequent steps. It is acceptable to set the same output folder for all phenotype files, as long as each phenotype name remains unique. Adjust glm_2_ss.py. i. Update the reference files specified at the beginning of the glm_2_ss.py script. Ensure that the primary inputs for this stage, derived from the previous ‘Reference file preparations’ step, are correctly assigned their respective file paths before moving on to the next phase. #outputs of the previous step index = pd.read_csv('/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt',sep='\t') freqfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out.afreq' Update the reference files specified at the beginning of the glm_2_ss.py script. Ensure that the primary inputs for this stage, derived from the previous ‘Reference file preparations’ step, are correctly assigned their respective file paths before moving on to the next phase. #outputs of the previous step index = pd.read_csv('/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt',sep='\t') freqfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out.afreq' #outputs of the previous step index = pd.read_csv('/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/chrpos.rs.txt',sep='\t') freqfile = '/path_to/4.Heritability_and_genetic_correlation_analysis/LDSC_ref/STITCH_out.afreq' Modify the 4.2.Conversion.sh script: #Python version3 Python=/path_to/Python3 #The python codes provided in github glm_2_ss=/path_to/4.Heritability_and_genetic_correlation_analysis/glm_2_ss.py ss_2_ldsc=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_ldsc.py ss_2_snplist=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_snplist.py #Python version3 Python=/path_to/Python3 #The python codes provided in github glm_2_ss=/path_to/4.Heritability_and_genetic_correlation_analysis/glm_2_ss.py ss_2_ldsc=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_ldsc.py ss_2_snplist=/path_to/4.Heritability_and_genetic_correlation_analysis/ss_2_snplist.py Run the code, the external parameters include GWAS result data, phenotype names, and the designated output directory: sh 4.2.Conversion.sh /path_to/pheno1.add pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out CRITICAL: Please note that this step requires Python version 3. Be sure not to confuse the two versions. Note: Each GWAS result needs this transformation to calculate heritability and genetic correlation. We recommend converting all GWAS summary statistics in this step to facilitate subsequent genetic correlation analyses. Note: In this step, the output files will have the extensions ‘.ss’, ‘.snplist’, and ‘.ldsc’. The files with ‘.snplist’ and ‘.ldsc’ extensions will be used in subsequent steps. It is acceptable to set the same output folder for all phenotype files, as long as each phenotype name remains unique. sh 4.2.Conversion.sh /path_to/pheno1.add pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out CRITICAL: Please note that this step requires Python version 3. Be sure not to confuse the two versions. Note: Each GWAS result needs this transformation to calculate heritability and genetic correlation. We recommend converting all GWAS summary statistics in this step to facilitate subsequent genetic correlation analyses. Note: In this step, the output files will have the extensions ‘.ss’, ‘.snplist’, and ‘.ldsc’. The files with ‘.snplist’ and ‘.ldsc’ extensions will be used in subsequent steps. It is acceptable to set the same output folder for all phenotype files, as long as each phenotype name remains unique. Heritability analysis. a. Specify the arguments in 4.3.Heritability_analysis.sh. i. Input the previously obtained LD reference files for the parameters ‘--ref-ld-chr’ and ‘--w-ld-chr’ : #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation munge_sumstats=/path_to/ldsc-master/munge_sumstats.py ldsc=/path_to/ldsc-master/ldsc.py ii. After that, run LDSC with: sh 4.3.Heritability_analysis.sh pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out /path_to/4.Heritability_and_genetic_correlation_analysis/heritability CRITICAL: Please note that all Python instances within the script 4.3.Heritability_analysis.sh utilize Python 2. It is essential to adjust the version accordingly and ensure that the correct path for the LDSC software functions is specified. Specify the arguments in 4.3.Heritability_analysis.sh. i. Input the previously obtained LD reference files for the parameters ‘--ref-ld-chr’ and ‘--w-ld-chr’ : #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation munge_sumstats=/path_to/ldsc-master/munge_sumstats.py ldsc=/path_to/ldsc-master/ldsc.py ii. After that, run LDSC with: sh 4.3.Heritability_analysis.sh pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out /path_to/4.Heritability_and_genetic_correlation_analysis/heritability CRITICAL: Please note that all Python instances within the script 4.3.Heritability_analysis.sh utilize Python 2. It is essential to adjust the version accordingly and ensure that the correct path for the LDSC software functions is specified. Input the previously obtained LD reference files for the parameters ‘--ref-ld-chr’ and ‘--w-ld-chr’ : #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation munge_sumstats=/path_to/ldsc-master/munge_sumstats.py ldsc=/path_to/ldsc-master/ldsc.py #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation munge_sumstats=/path_to/ldsc-master/munge_sumstats.py ldsc=/path_to/ldsc-master/ldsc.py After that, run LDSC with: sh 4.3.Heritability_analysis.sh pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out /path_to/4.Heritability_and_genetic_correlation_analysis/heritability CRITICAL: Please note that all Python instances within the script 4.3.Heritability_analysis.sh utilize Python 2. It is essential to adjust the version accordingly and ensure that the correct path for the LDSC software functions is specified. sh 4.3.Heritability_analysis.sh pheno1 /path_to/4.Heritability_and_genetic_correlation_analysis/out /path_to/4.Heritability_and_genetic_correlation_analysis/heritability CRITICAL: Please note that all Python instances within the script 4.3.Heritability_analysis.sh utilize Python 2. It is essential to adjust the version accordingly and ensure that the correct path for the LDSC software functions is specified. Genetic correlation analysis. a. Specify the arguments in 4.4.Genetic_correlation.sh: #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py b. Then run the code: sh 4.4.Genetic_correlation.sh pheno_name1 pheno_name2 /path_to/4.Heritability_and_genetic_correlation_analysis/heritability /path_to/4.Heritability_and_genetic_correlation_analysis/correlation CRITICAL: Please note that this step requires Python version 2. Make sure not to confuse it with Python version 3. Specify the arguments in 4.4.Genetic_correlation.sh: #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py #Python version2 Python=/path_to/Python2 #You can check https://github.com/bulik/ldsc?tab=readme-ov-file#where-can-i-get-ld-scores #for LD reference score files LD_ref_path=/path_to/Supporting_data/LDSC/ #ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py Then run the code: sh 4.4.Genetic_correlation.sh pheno_name1 pheno_name2 /path_to/4.Heritability_and_genetic_correlation_analysis/heritability /path_to/4.Heritability_and_genetic_correlation_analysis/correlation CRITICAL: Please note that this step requires Python version 2. Make sure not to confuse it with Python version 3. sh 4.4.Genetic_correlation.sh pheno_name1 pheno_name2 /path_to/4.Heritability_and_genetic_correlation_analysis/heritability /path_to/4.Heritability_and_genetic_correlation_analysis/correlation CRITICAL: Please note that this step requires Python version 2. Make sure not to confuse it with Python version 3. Timing: 1 day Timing: 1 day In this section, we focus on Partitioned Heritability analysis, which uses the LDSC software along with partitioned annotation reference files from various cell and tissue types derived from different databases. This analysis decomposes heritability across different cell/tissue types to assess each tissue’s genetic contribution to the heritability of the phenotype of interest. First, we need to download the cell-type reference files provided by the official LDSC software. Then, we will use the summary statistics in the required LDSC format as the input to calculate partitioned heritability for each cell or tissue type. The final output will include coefficients and p -values for each tissue type. 14. Reference file preparation. a. Make sure that all required reference files have been successfully downloaded and placed in the correct folder. 15. Partitioned heritability calculation. Note: When calculating partitioned heritability, it is essential to use summary statistics files that meet LDSC requirements as input. Specifically, these are the ‘.sumstats.gz’ files generated in the previous calculations of heritability and genetic correlations. a. Modify the 5.1.Partitioned_heritability.sh before performing the calculation: # Partitioned_heritability.sh #Python version2 Python=/path_to/Python2 #for LD reference score files LD_PH_ref_path=/path_to/Required_reference_data/LDSC/ # ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py b. Execute the code with: sh /path_to/5.1.Partitioned_heritability.sh Traits.sumstats.gz Traits.ph2 Note: The result of the calculation is a file with the suffix ‘cell_type_results.txt’. Each line in the file represents a cell or tissue type, with four columns corresponding to the tissue name, coefficient, standard error of the coefficient, and the coefficient's p-value. Note: The reference files listed above were constructed based on the East Asian population. If your study focuses on other populations, please download the corresponding reference files for those populations. Reference file preparation. a. Make sure that all required reference files have been successfully downloaded and placed in the correct folder. Make sure that all required reference files have been successfully downloaded and placed in the correct folder. Partitioned heritability calculation. Note: When calculating partitioned heritability, it is essential to use summary statistics files that meet LDSC requirements as input. Specifically, these are the ‘.sumstats.gz’ files generated in the previous calculations of heritability and genetic correlations. a. Modify the 5.1.Partitioned_heritability.sh before performing the calculation: # Partitioned_heritability.sh #Python version2 Python=/path_to/Python2 #for LD reference score files LD_PH_ref_path=/path_to/Required_reference_data/LDSC/ # ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py b. Execute the code with: sh /path_to/5.1.Partitioned_heritability.sh Traits.sumstats.gz Traits.ph2 Note: The result of the calculation is a file with the suffix ‘cell_type_results.txt’. Each line in the file represents a cell or tissue type, with four columns corresponding to the tissue name, coefficient, standard error of the coefficient, and the coefficient's p-value. Note: The reference files listed above were constructed based on the East Asian population. If your study focuses on other populations, please download the corresponding reference files for those populations. Note: When calculating partitioned heritability, it is essential to use summary statistics files that meet LDSC requirements as input. Specifically, these are the ‘.sumstats.gz’ files generated in the previous calculations of heritability and genetic correlations. Modify the 5.1.Partitioned_heritability.sh before performing the calculation: # Partitioned_heritability.sh #Python version2 Python=/path_to/Python2 #for LD reference score files LD_PH_ref_path=/path_to/Required_reference_data/LDSC/ # ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py # Partitioned_heritability.sh #Python version2 Python=/path_to/Python2 #for LD reference score files LD_PH_ref_path=/path_to/Required_reference_data/LDSC/ # ldsc codes, can be found after LDSC installation ldsc=/path_to/ldsc-master/ldsc.py Execute the code with: sh /path_to/5.1.Partitioned_heritability.sh Traits.sumstats.gz Traits.ph2 Note: The result of the calculation is a file with the suffix ‘cell_type_results.txt’. Each line in the file represents a cell or tissue type, with four columns corresponding to the tissue name, coefficient, standard error of the coefficient, and the coefficient's p-value. Note: The reference files listed above were constructed based on the East Asian population. If your study focuses on other populations, please download the corresponding reference files for those populations. sh /path_to/5.1.Partitioned_heritability.sh Traits.sumstats.gz Traits.ph2 Note: The result of the calculation is a file with the suffix ‘cell_type_results.txt’. Each line in the file represents a cell or tissue type, with four columns corresponding to the tissue name, coefficient, standard error of the coefficient, and the coefficient's p-value. Note: The reference files listed above were constructed based on the East Asian population. If your study focuses on other populations, please download the corresponding reference files for those populations. Timing: <1 week Timing: <1 week In our published articles, 1 , 2 we used the summary statistics from our GWAS results as the exposure and the 46 GWAS results from the BBJ 3 as the outcome. In this section, we will use the content of our article as a case study to illustrate how to conduct two-sample Mendelian randomization (MR) analyses between phenotypes in two different groups. The results can ultimately be presented in the form of tables or heat maps. In this section, we provide a GWAS result on Google Drive ( https://drive.google.com/file/d/1R2hrZwiTlz-kVC4F_zf1M8sz0mMcsC1L/view?usp=sharing ), which can serve as the exposure data. You may convert it into the appropriate format following the instructions from the previous step for testing purposes. Additionally, we provide download links for the 46 outcome datasets used in our published studies. You can select the data you need by referring to the corresponding descriptions. 16. Prepare GWAS summary statistics files as exposure. Note: In this section, the summary statistics files we will use were generated during the ‘Heritability and Genetic Correlation Analysis’. If you have already completed the heritability analysis, you can consolidate all of the ‘.ss’ files into a single summary file, as demonstrated in ‘example_exposure_summaries.txt’ using: ls /path_to/4.Heritability_and_genetic_correlation_analysis/out/∗.ss > exposure_summaries.txt Alternatives: If you did not conduct the heritability or genetic correlation analysis, please refer to the first step of the ‘Heritability and Genetic Correlation Analysis’ to convert the necessary files. 17. Download external summary statistic files as the outcome. Note: In our published work, 1 we used our GWAS results for pregnancy phenotypes as exposure variables and external GWAS results from BioBank Japan (BBJ) 3 as outcome variables to conduct Mendelian randomization analysis. It is important to note that the GWAS summary statistics from BBJ are not the only option available. Users can download and use any publicly available GWAS summary statistics or their own data, depending on the specific needs of the study. In this section, we use BBJ's GWAS summary data as an example of outcome data. To obtain this data, please download them from http://jenger.riken.jp/result . Instead of converting the downloaded results, we will provide the case and control numbers externally. a. Find the case and control numbers on the BBJ webpage. b. Summarize the paths of the downloaded files along with the corresponding case and control numbers into a single summary file, as shown in example_outcome_summaries.txt: ID Study                 Case  Ctrl 1 /path_to/outcome/phenocode-Atrial_Fibrillation.tsv.gz    8180  28612 2 /path_to/outcome/phenocode-Rheumatoid_Arthritis_Asian.tsv.gz  3636  15554 8 /path_to/outcome/phenocode-Open-angle_glaucoma_autosome.tsv.gz 3980  18815 11 /path_to/outcome/phenocode-Colorectal_Cancer.tsv.gz   6692  27178 14 /path_to/outcome/phenocode-Type_2_Diabetes.tsv.gz    36614 155150 Alternatives: The selection of outcome data is not limited to summary statistics provided by biobanks; it can also be sourced from other databases. Please adjust the input method as needed to match the format of summary data from these alternative sources. Notably, the Mendelian randomization R package, TwoSampleMR, used in this study, supports direct access to data from the GWAS Catalog. For more information, please refer to: https://mrcieu.github.io/TwoSampleMR/ . 18. Perform Mendelian randomization analysis. a. Add the paths of these two files into 6.1.Mendelian_Randomization.R and run the code. #exposure .ss data exposure = read.table("/path_to/5.Mendelian_randomization_analysis/exposure_summaries.txt",header = T)[,1] #outcome summary statistics outcome = read.table("/path_to/5.Mendelian_randomization_analysis/outcome_summaries.txt",header = T) Note: The output will be a list of MR p-values generated using various Mendelian randomization methods. Note: In MR analysis, we need to perform 'clumping' on the GWAS summary statistics. In the R package TwoSampleMR, this function is called clump_data() . Note that using this function requires an internet connection, and access to the OpenGWAS server can be unstable in certain regions. If you encounter internet connectivity issues or prefer to use local linkage disequilibrium (LD) data for clumping, you can switch to the ld_clump() function from the ieugwasr package. For more information on local LD clumping, please refer to the following resource: https://mrcieu.github.io/ieugwasr/articles/local_ld.html . Prepare GWAS summary statistics files as exposure. Note: In this section, the summary statistics files we will use were generated during the ‘Heritability and Genetic Correlation Analysis’. If you have already completed the heritability analysis, you can consolidate all of the ‘.ss’ files into a single summary file, as demonstrated in ‘example_exposure_summaries.txt’ using: ls /path_to/4.Heritability_and_genetic_correlation_analysis/out/∗.ss > exposure_summaries.txt Alternatives: If you did not conduct the heritability or genetic correlation analysis, please refer to the first step of the ‘Heritability and Genetic Correlation Analysis’ to convert the necessary files. Download external summary statistic files as the outcome. Note: In our published work, 1 we used our GWAS results for pregnancy phenotypes as exposure variables and external GWAS results from BioBank Japan (BBJ) 3 as outcome variables to conduct Mendelian randomization analysis. It is important to note that the GWAS summary statistics from BBJ are not the only option available. Users can download and use any publicly available GWAS summary statistics or their own data, depending on the specific needs of the study. In this section, we use BBJ's GWAS summary data as an example of outcome data. To obtain this data, please download them from http://jenger.riken.jp/result . Instead of converting the downloaded results, we will provide the case and control numbers externally. a. Find the case and control numbers on the BBJ webpage. b. Summarize the paths of the downloaded files along with the corresponding case and control numbers into a single summary file, as shown in example_outcome_summaries.txt: ID Study                 Case  Ctrl 1 /path_to/outcome/phenocode-Atrial_Fibrillation.tsv.gz    8180  28612 2 /path_to/outcome/phenocode-Rheumatoid_Arthritis_Asian.tsv.gz  3636  15554 8 /path_to/outcome/phenocode-Open-angle_glaucoma_autosome.tsv.gz 3980  18815 11 /path_to/outcome/phenocode-Colorectal_Cancer.tsv.gz   6692  27178 14 /path_to/outcome/phenocode-Type_2_Diabetes.tsv.gz    36614 155150 Alternatives: The selection of outcome data is not limited to summary statistics provided by biobanks; it can also be sourced from other databases. Please adjust the input method as needed to match the format of summary data from these alternative sources. Notably, the Mendelian randomization R package, TwoSampleMR, used in this study, supports direct access to data from the GWAS Catalog. For more information, please refer to: https://mrcieu.github.io/TwoSampleMR/ . Note: In our published work, 1 we used our GWAS results for pregnancy phenotypes as exposure variables and external GWAS results from BioBank Japan (BBJ) 3 as outcome variables to conduct Mendelian randomization analysis. It is important to note that the GWAS summary statistics from BBJ are not the only option available. Users can download and use any publicly available GWAS summary statistics or their own data, depending on the specific needs of the study. In this section, we use BBJ's GWAS summary data as an example of outcome data. To obtain this data, please download them from http://jenger.riken.jp/result . Instead of converting the downloaded results, we will provide the case and control numbers externally. Find the case and control numbers on the BBJ webpage. Summarize the paths of the downloaded files along with the corresponding case and control numbers into a single summary file, as shown in example_outcome_summaries.txt: ID Study                 Case  Ctrl 1 /path_to/outcome/phenocode-Atrial_Fibrillation.tsv.gz    8180  28612 2 /path_to/outcome/phenocode-Rheumatoid_Arthritis_Asian.tsv.gz  3636  15554 8 /path_to/outcome/phenocode-Open-angle_glaucoma_autosome.tsv.gz 3980  18815 11 /path_to/outcome/phenocode-Colorectal_Cancer.tsv.gz   6692  27178 14 /path_to/outcome/phenocode-Type_2_Diabetes.tsv.gz    36614 155150 Alternatives: The selection of outcome data is not limited to summary statistics provided by biobanks; it can also be sourced from other databases. Please adjust the input method as needed to match the format of summary data from these alternative sources. Notably, the Mendelian randomization R package, TwoSampleMR, used in this study, supports direct access to data from the GWAS Catalog. For more information, please refer to: https://mrcieu.github.io/TwoSampleMR/ . ID Study                 Case  Ctrl 1 /path_to/outcome/phenocode-Atrial_Fibrillation.tsv.gz    8180  28612 2 /path_to/outcome/phenocode-Rheumatoid_Arthritis_Asian.tsv.gz  3636  15554 8 /path_to/outcome/phenocode-Open-angle_glaucoma_autosome.tsv.gz 3980  18815 11 /path_to/outcome/phenocode-Colorectal_Cancer.tsv.gz   6692  27178 14 /path_to/outcome/phenocode-Type_2_Diabetes.tsv.gz    36614 155150 Alternatives: The selection of outcome data is not limited to summary statistics provided by biobanks; it can also be sourced from other databases. Please adjust the input method as needed to match the format of summary data from these alternative sources. Notably, the Mendelian randomization R package, TwoSampleMR, used in this study, supports direct access to data from the GWAS Catalog. For more information, please refer to: https://mrcieu.github.io/TwoSampleMR/ . Perform Mendelian randomization analysis. a. Add the paths of these two files into 6.1.Mendelian_Randomization.R and run the code. Add the paths of these two files into 6.1.Mendelian_Randomization.R and run the code. #exposure .ss data exposure = read.table("/path_to/5.Mendelian_randomization_analysis/exposure_summaries.txt",header = T)[,1] #outcome summary statistics outcome = read.table("/path_to/5.Mendelian_randomization_analysis/outcome_summaries.txt",header = T) Note: The output will be a list of MR p-values generated using various Mendelian randomization methods. Note: In MR analysis, we need to perform 'clumping' on the GWAS summary statistics. In the R package TwoSampleMR, this function is called clump_data() . Note that using this function requires an internet connection, and access to the OpenGWAS server can be unstable in certain regions. If you encounter internet connectivity issues or prefer to use local linkage disequilibrium (LD) data for clumping, you can switch to the ld_clump() function from the ieugwasr package. For more information on local LD clumping, please refer to the following resource: https://mrcieu.github.io/ieugwasr/articles/local_ld.html . Timing: 2 h ours Timing: 2 h ours In our published article on pregnancy glycemic traits, 2 we conducted a Transcriptome-Wide Association Study (TWAS) using the gene ESR1. This analysis utilized expression data downloaded from GTEx, along with genetic variants and phenotypic data from our samples. In this section, we will use the findings from our published work as an example to illustrate the process of performing TWAS. The outcome will be the p -value of the analysis. 19. Prepare the files for analyzing. a. Prepare the expression data in different tissues from https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . We use data from subcutaneous adipose tissue as an example. b. Run the code 7.1.Expression_prepare.R after adjusting the input file paths: #load the required db file we downloaded filename <- paste0("/path_to/ctimp_Adipose_Subcutaneous.db") c. Prepare the PLINK2 formatted variant file. Note: In this section, we will utilize the ‘.vcf’ file generated in Step 5 (the third step of ‘Genotype Imputation Using STITCH’). You can obtain the variants associated with the target gene using the following command: bcftools view -r chr6:151640495-152153274 vcf_file > output_vcf plink2 --export Av --vcf output_vcf -out ESR1.traw d. The phenotypes of the samples in ‘.vcf’ file. Note: Here, we can utilize the phenotype file from the GWAS process. Please refer to phenotype_example.txt for details. 20. Execute the 7.2.TWAS.R script. a. Modify the input of the script: ##disease gdm <- as.data.frame(fread("pheno_quantitative.csv",header=T)) ##indiviual-level genotype matrix on ESR1 traw <- as.data.frame(fread("ESR1.traw",header = T)) #the formatted expression data eqtl <- as.data.frame(fread(paste0("ctimp_Adipose_Subcutaneous.db.txt"),header = T)) Note: This protocol illustrates the results of a GWAS analysis focused on a single gene within a specific tissue type, yielding a single outcome. Readers are encouraged to modify the code as needed to adapt to analytical scenarios involving multiple tissues and genes. Prepare the files for analyzing. a. Prepare the expression data in different tissues from https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . We use data from subcutaneous adipose tissue as an example. b. Run the code 7.1.Expression_prepare.R after adjusting the input file paths: #load the required db file we downloaded filename <- paste0("/path_to/ctimp_Adipose_Subcutaneous.db") c. Prepare the PLINK2 formatted variant file. Note: In this section, we will utilize the ‘.vcf’ file generated in Step 5 (the third step of ‘Genotype Imputation Using STITCH’). You can obtain the variants associated with the target gene using the following command: bcftools view -r chr6:151640495-152153274 vcf_file > output_vcf plink2 --export Av --vcf output_vcf -out ESR1.traw d. The phenotypes of the samples in ‘.vcf’ file. Note: Here, we can utilize the phenotype file from the GWAS process. Please refer to phenotype_example.txt for details. Prepare the expression data in different tissues from https://predictdb.org/post/2020/09/10/gtex-v8-cross-tissue-imputation-utmost-prediction-models/ . We use data from subcutaneous adipose tissue as an example. Run the code 7.1.Expression_prepare.R after adjusting the input file paths: #load the required db file we downloaded filename <- paste0("/path_to/ctimp_Adipose_Subcutaneous.db") #load the required db file we downloaded filename <- paste0("/path_to/ctimp_Adipose_Subcutaneous.db") Prepare the PLINK2 formatted variant file. Note: In this section, we will utilize the ‘.vcf’ file generated in Step 5 (the third step of ‘Genotype Imputation Using STITCH’). You can obtain the variants associated with the target gene using the following command: bcftools view -r chr6:151640495-152153274 vcf_file > output_vcf plink2 --export Av --vcf output_vcf -out ESR1.traw Note: In this section, we will utilize the ‘.vcf’ file generated in Step 5 (the third step of ‘Genotype Imputation Using STITCH’). You can obtain the variants associated with the target gene using the following command: bcftools view -r chr6:151640495-152153274 vcf_file > output_vcf plink2 --export Av --vcf output_vcf -out ESR1.traw The phenotypes of the samples in ‘.vcf’ file. Note: Here, we can utilize the phenotype file from the GWAS process. Please refer to phenotype_example.txt for details. Note: Here, we can utilize the phenotype file from the GWAS process. Please refer to phenotype_example.txt for details. Execute the 7.2.TWAS.R script. a. Modify the input of the script: Modify the input of the script: ##disease gdm <- as.data.frame(fread("pheno_quantitative.csv",header=T)) ##indiviual-level genotype matrix on ESR1 traw <- as.data.frame(fread("ESR1.traw",header = T)) #the formatted expression data eqtl <- as.data.frame(fread(paste0("ctimp_Adipose_Subcutaneous.db.txt"),header = T)) Note: This protocol illustrates the results of a GWAS analysis focused on a single gene within a specific tissue type, yielding a single outcome. Readers are encouraged to modify the code as needed to adapt to analytical scenarios involving multiple tissues and genes. Timing: 1–2 days Timing: 1–2 days In this step, we will conduct drug target enrichment analysis, as outlined in the article on pregnancy glycemic traits. 2 We primarily utilize the Anatomical Therapeutic Chemical (ATC) classification system for drug categorization, along with the Drug-Gene Interaction Database (DGIdb) version 5 interactions file as our inputs. The drug target analysis will be performed using MAGMA software, and we will visualize the results by plotting the area under the curve (AUC). 21. MAGMA annotation. Note: This step involves initializing the MAGMA software by downloading the required data files from its website ( https://cncr.nl/research/magma/ ) for preliminary annotation. Please confirm that the files have been successfully downloaded. a. Run the script 8.1.MAGMA_annotation.sh after adjusting the arguments: #magma_annotation magma=/path_to/magma #data required for initializing MAGMA g1000=/path_to/Required_reference_data/magma/g1000_eas.bim gene_loc=/path_to/Required_reference_data/magma/NCBI37.3.gene.loc Note: Please download the relevant files corresponding to the population used in this analysis. 22. Gene analysis with SNP p-values, modify the 8.2.gene_analysis.sh. Note: In this section, we will use summary statistics files generated in the first step of ‘Heritability and Genetic Correlation Analysis.’ If you have already conducted the heritability analysis, you can include the corresponding '.ss' file in the '--pval' argument. The input files required for this step include the 'bfile' downloaded in the last step, the GWAS summary statistics file for the phenotype of interest, and the output file from the last step. #step2_Gene analysis with SNP p-values magma=/path_to/magma g1000=/path_to/Required_reference_data/magma/g1000_eas #the output of 8.1.MAGMA_annotation.sh genes=/path_to/example_annotation.genes.annot #ss file Phenotype=/path_to/pheno.ss #use for SNP and pvalue, ncol for samplesize $magma \ --bfile $g1000 \ --pval $Phenotype use=4,14 ncol=10 \ --gene-annot $genes \ --gene-model multi \ --out GWAS_pheno_example Note: The 'use = ' argument of ‘--pval’ specifies the columns corresponding to the SNP and p-value, while the 'ncol = ' argument indicates the column representing the sample size. 23. Gene-set analysis. Note: In this step, we will input the previously downloaded DGIdb version 5 interactions file into MAGMA after making the necessary formatting adjustments. The input for this step includes the results from previous step and the DGIdb interactions file. a. Run the 8.3.Gene_set_analysis.sh script after filling in the required information: #step3_Gene_set_analysis magma=/path_to/magma interactions=/path_to/Required_reference_data/magma/interactions.tsv #convert the interactions.tsv cut -f 3,9 $interactions |awk '{if($2!=""&&$1!="") print $_}' - > interactions_sorted && \ sed -i '1d' interactions_sorted && \ sort -k2 interactions_sorted > interactions_sorted #gene set analysis $magma \ --gene-results /path_to/GWAS_pheno_example.genes.raw \ --set-annot interactions_sorted col=1,2 \ --out GWAS_pheno_example.drug 24. AUC calculation. a. Execute 8.4.AUC_calculation.R. Note: In this step, we will summarize all results and calculate the AUC. We need to match and classify the previously obtained files from the ATC classification system website with our prior results and the DGIdb version 5 interactions file. Finally, we will calculate the AUC for potential drugs and plot the results. This step primarily involves three input files. You can edit the relevant sections in the 8.4.AUC_calculation.R script and run it directly, or execute it in RStudio on your local device. #load the required files #the interaction file downloaded from DGIdb interact <- fread('/path_to/Required_reference_data/magma/interactions.tsv') #the ATC summary file db <- fread('/path_to/Required_reference_data/magma/atc_cata.txt') #the output gsa.out file produced by step 3 file_lines <- readLines('GWAS_pheno_example.drug.gsa.out') Note: The first input file for 8.4.AUC_calculation.R is the DGIdb interactions file. The second input file is the ATC classification file, which is typically acquired from a website as an Excel document. You can convert this file to a CSV format or any other format compatible with the fread() function from the R data.table package. The final input file is the output from the third step; please select the file with the '.gsa.out' extension as the input. MAGMA annotation. Note: This step involves initializing the MAGMA software by downloading the required data files from its website ( https://cncr.nl/research/magma/ ) for preliminary annotation. Please confirm that the files have been successfully downloaded. a. Run the script 8.1.MAGMA_annotation.sh after adjusting the arguments: #magma_annotation magma=/path_to/magma #data required for initializing MAGMA g1000=/path_to/Required_reference_data/magma/g1000_eas.bim gene_loc=/path_to/Required_reference_data/magma/NCBI37.3.gene.loc Note: Please download the relevant files corresponding to the population used in this analysis. Note: This step involves initializing the MAGMA software by downloading the required data files from its website ( https://cncr.nl/research/magma/ ) for preliminary annotation. Please confirm that the files have been successfully downloaded. Run the script 8.1.MAGMA_annotation.sh after adjusting the arguments: #magma_annotation magma=/path_to/magma #data required for initializing MAGMA g1000=/path_to/Required_reference_data/magma/g1000_eas.bim gene_loc=/path_to/Required_reference_data/magma/NCBI37.3.gene.loc Note: Please download the relevant files corresponding to the population used in this analysis. #magma_annotation magma=/path_to/magma #data required for initializing MAGMA g1000=/path_to/Required_reference_data/magma/g1000_eas.bim gene_loc=/path_to/Required_reference_data/magma/NCBI37.3.gene.loc Note: Please download the relevant files corresponding to the population used in this analysis. Gene analysis with SNP p-values, modify the 8.2.gene_analysis.sh. Note: In this section, we will use summary statistics files generated in the first step of ‘Heritability and Genetic Correlation Analysis.’ If you have already conducted the heritability analysis, you can include the corresponding '.ss' file in the '--pval' argument. The input files required for this step include the 'bfile' downloaded in the last step, the GWAS summary statistics file for the phenotype of interest, and the output file from the last step. #step2_Gene analysis with SNP p-values magma=/path_to/magma g1000=/path_to/Required_reference_data/magma/g1000_eas #the output of 8.1.MAGMA_annotation.sh genes=/path_to/example_annotation.genes.annot #ss file Phenotype=/path_to/pheno.ss #use for SNP and pvalue, ncol for samplesize $magma \ --bfile $g1000 \ --pval $Phenotype use=4,14 ncol=10 \ --gene-annot $genes \ --gene-model multi \ --out GWAS_pheno_example Note: The 'use = ' argument of ‘--pval’ specifies the columns corresponding to the SNP and p-value, while the 'ncol = ' argument indicates the column representing the sample size. Gene-set analysis. Note: In this step, we will input the previously downloaded DGIdb version 5 interactions file into MAGMA after making the necessary formatting adjustments. The input for this step includes the results from previous step and the DGIdb interactions file. a. Run the 8.3.Gene_set_analysis.sh script after filling in the required information: #step3_Gene_set_analysis magma=/path_to/magma interactions=/path_to/Required_reference_data/magma/interactions.tsv #convert the interactions.tsv cut -f 3,9 $interactions |awk '{if($2!=""&&$1!="") print $_}' - > interactions_sorted && \ sed -i '1d' interactions_sorted && \ sort -k2 interactions_sorted > interactions_sorted #gene set analysis $magma \ --gene-results /path_to/GWAS_pheno_example.genes.raw \ --set-annot interactions_sorted col=1,2 \ --out GWAS_pheno_example.drug Note: In this step, we will input the previously downloaded DGIdb version 5 interactions file into MAGMA after making the necessary formatting adjustments. The input for this step includes the results from previous step and the DGIdb interactions file. Run the 8.3.Gene_set_analysis.sh script after filling in the required information: #step3_Gene_set_analysis magma=/path_to/magma interactions=/path_to/Required_reference_data/magma/interactions.tsv #convert the interactions.tsv cut -f 3,9 $interactions |awk '{if($2!=""&&$1!="") print $_}' - > interactions_sorted && \ sed -i '1d' interactions_sorted && \ sort -k2 interactions_sorted > interactions_sorted #gene set analysis $magma \ --gene-results /path_to/GWAS_pheno_example.genes.raw \ --set-annot interactions_sorted col=1,2 \ --out GWAS_pheno_example.drug #step3_Gene_set_analysis magma=/path_to/magma interactions=/path_to/Required_reference_data/magma/interactions.tsv #convert the interactions.tsv cut -f 3,9 $interactions |awk '{if($2!=""&&$1!="") print $_}' - > interactions_sorted && \ sed -i '1d' interactions_sorted && \ sort -k2 interactions_sorted > interactions_sorted #gene set analysis $magma \ --gene-results /path_to/GWAS_pheno_example.genes.raw \ --set-annot interactions_sorted col=1,2 \ --out GWAS_pheno_example.drug AUC calculation. a. Execute 8.4.AUC_calculation.R. Execute 8.4.AUC_calculation.R. Note: In this step, we will summarize all results and calculate the AUC. We need to match and classify the previously obtained files from the ATC classification system website with our prior results and the DGIdb version 5 interactions file. Finally, we will calculate the AUC for potential drugs and plot the results. This step primarily involves three input files. You can edit the relevant sections in the 8.4.AUC_calculation.R script and run it directly, or execute it in RStudio on your local device. #load the required files #the interaction file downloaded from DGIdb interact <- fread('/path_to/Required_reference_data/magma/interactions.tsv') #the ATC summary file db <- fread('/path_to/Required_reference_data/magma/atc_cata.txt') #the output gsa.out file produced by step 3 file_lines <- readLines('GWAS_pheno_example.drug.gsa.out') Note: The first input file for 8.4.AUC_calculation.R is the DGIdb interactions file. The second input file is the ATC classification file, which is typically acquired from a website as an Excel document. You can convert this file to a CSV format or any other format compatible with the fread() function from the R data.table package. The final input file is the output from the third step; please select the file with the '.gsa.out' extension as the input.

The authors declare no competing interests.

This study was supported by the National Natural Science Foundation of China (No. 32171441 ).

When executing the script 2.2.run_STITCH.sh, a common issue that may arise is illustrated below: [2024-09-30 00:23:12] There are 0 variants in the left buffer region -249999 <= position < 1 [2024-09-30 00:23:12] There are 0 variants in the central region 1 <= position <= 5000000 [2024-09-30 00:23:12] There are 0 variants in the right buffer region 5000000 < position <= 5250000 Error in validate_region_to_impute_when_using_regionStart(L, regionStart, : There are fewer than 2 SNPs to impute, i.e., there is 1 SNP to impute. In this case, imputation is really just genotyping. STITCH could support genotyping but does not, and note that this kind of defeats the point of imputation. Please use your favourite genotyper e.g. GATK to genotype these SNPs. If you strongly disagree please file a bug report and this can be re-examined Calls: STITCH ... shrink_region -> validate_region_to_impute_when_using_regionStart Execution halted [2024-09-30 00:23:12] There are 0 variants in the left buffer region -249999 <= position < 1 [2024-09-30 00:23:12] There are 0 variants in the central region 1 <= position <= 5000000 [2024-09-30 00:23:12] There are 0 variants in the right buffer region 5000000 < position <= 5250000 Error in validate_region_to_impute_when_using_regionStart(L, regionStart, : There are fewer than 2 SNPs to impute, i.e., there is 1 SNP to impute. In this case, imputation is really just genotyping. STITCH could support genotyping but does not, and note that this kind of defeats the point of imputation. Please use your favourite genotyper e.g. GATK to genotype these SNPs. If you strongly disagree please file a bug report and this can be re-examined Calls: STITCH ... shrink_region -> validate_region_to_impute_when_using_regionStart Execution halted Technically, this is not an error but rather an expected response from STITCH when it detects no variant loci within the specified input range. Such intervals are typically found in the centromeric or telomeric regions of chromosomes, where the absence of read coverage is normal. If you encounter this situation, it is advisable to remove these intervals during the final file merging process. However, if these occurrences are unexpected within your dataset, it is recommended to validate your results using alternative software. For any anomalies, you can initiate discussions on the STITCH GitHub page. When executing the script 2.2.run_STITCH.sh, you may encounter another type of error, as illustrated below: ∗∗∗ caught segfault ∗∗∗ address 0x9b812004, cause 'memory not mapped' ∗∗∗ caught segfault ∗∗∗ address 0x9b812004, cause 'memory not mapped' As noted by the software’s developer in issue #86 ( https://github.com/rwdavies/STITCH/issues/86 ) and related discussions, this error can occur intermittently due to a known problem with the multithreading package used by STITCH. Despite the developer’s ongoing attempts to address this issue, it remains unresolved. A practical workaround we have identified involves adjusting the '--ncores' parameter in the 2.2.run_STITCH.sh script (as illustrated in the figure). If an error arises while using a single core, switching to two or more cores, or vice versa, can often resolve most of these issues upon re-execution.