Single-cell multiomics of neuronal activation reveals context-dependent genetic control of brain disorders
Data files
Mar 02, 2025 version files 318.36 MB
-
All_other_supplementary_Tables_(1-3__5__7__8_10_11__15__16-18_22_26-29_31_32_33).xlsx
28.74 MB
-
README.md
27.42 KB
-
Table_S12_ABC_EnhancerPredictions_threshold0.021_self_promoter.xlsx
62.06 MB
-
Table_S13._Differential_TF_motif_activity_analysis_using_a_one-tailed_Wilcoxon_signed-rank_test.xlsx
93.69 KB
-
Table_S14._Gene_regulatory_networks_identified_in_three_cell_types.xlsx
3.50 MB
-
Table_S19_Signficant_eQTLs_in_each_cell_type_and_time_point.xlsx
64.12 MB
-
Table_S20._Dynamic_cGenes_in_each_cell_type_and_time_point.xlsx
877.12 KB
-
Table_S21._eQTL-based_cTWAS_results.xlsx
70.51 KB
-
Table_S23._cPeak_and_caQTL_identified_by_tensorQTL_in_each_cell_type_and_time_point.xlsx
8.14 MB
-
Table_S24._ASoC_SNPs_in_each_cell_type_and_time_point.xlsx
113.96 MB
-
Table_S25._Dynamic_caQTLs_for_each_cell_type_with_permutation_(n_1000)_testing_result.xlsx
1.57 MB
-
Table_S30._caQTL-based_cTWAS_results.xlsx
82.04 KB
-
Table_S6._Differentially_accessible_peaks_upon_KCI_stimulation_table.xlsx
31.93 MB
-
Table_S9_OCR-gene_pairs_identified_through_single-cell_level_peak-gene_correlation_(co-activation)_analysis.xlsx
3.18 MB
Abstract
Despite hundreds of genetic risk loci identified for neuropsychiatric disorders (NPD), most causal variants/genes remain unknown. A major hurdle is that disease risk variants may act in specific biological contexts, e.g., during neuronal activation, which is difficult to study in vivo at the population level. Here, we modeled neuronal activation in human iPSC-induced excitatory and inhibitory neurons from 100 donors. Single-cell multiomics analyses of over a million neurons uncovered complex activity-dependent transcriptomic and epigenomic regulation and significantly expanded the repertoire of stimulation-specific causal variants/genes for NPD. We identified thousands of genetic variants associated with activity-dependent gene expression (i.e., eQTL) and chromatin accessibility (i.e., caQTL). These caQTL explained considerably larger proportions of NPD heritability than the eQTL. Integrating the multiomic data with GWAS revealed NPD risk variants/genes whose effects were only detected upon stimulation. Interestingly, multiple lines of evidence support a role of activity-dependent cholesterol metabolism in NPD. Our work highlights the power of cell stimulation to reveal context-dependent “hidden” genetic effects.
https://doi.org/10.5061/dryad.0zpc8677w
Description of the data and file structure
Single-cell multiomics of neuronal activation reveals context-dependent genetic control of brain disorders.
Lifan Liang1,†, Siwei Zhang2,5,†, Zicheng Wang1,†, Hanwen Zhang2,†, Chuxuan Li2,3,†, Alexandra C. Duhe2, Xiaotong Sun1, Xiaoyuan Zhong1, Alena Kozlova2, Brendan Jamison1,2, Whitney Wood2, Zhiping P. Pang4, Alan R. Sanders2,5, Xin He1,‡, Jubao Duan2,5,‡
Affiliations:
1Department of Human Genetics, The University of Chicago, Chicago, IL 60637, USA.
2Center for Psychiatric Genetics, Endeavor Health Research Institute, Evanston, IL 60201, USA.
3Graduate Group in Genomics and Computational Biology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA
4Department of Neuroscience and Cell Biology, Child Health Institute of New Jersey, Rutgers Robert Wood Johnson Medical School, New Brunswick, NJ 08901, USA.
5Department of Psychiatry and Behavioral Neuroscience, The University of Chicago, Chicago, IL 60637, USA.
† These authors contributed equally
‡ Corresponding authors: Xin He (xinhe@uchicago.edu), Jubao Duan (jduan@uchicago.edu)
Files and variables
File: Table_S9_OCR-gene_pairs_identified_through_single-cell_level_peak-gene_correlation_(co-activation)_analysis.xlsx
Description: Table S9. OCR-gene pairs identified through single-cell level peak-gene correlation (co-activation) analysis.
Variables
- Gene: Gene symbol;
- Peak: correlated peak position, chromosome:start-end;
- Beta: calculated coefficient from the linear model;
- SE: standard error calculated from the simple linear model;
- Z: z-score;
- Pvalue: P value;
- FDR: False Discovery Rate. The adjusted p-values corrected for multiple testing (Benjamini-Hochberg).
File: Table_S13._Differential_TF_motif_activity_analysis_using_a_one-tailed_Wilcoxon_signed-rank_test.xlsx
Description: Table S13. Differential TF motif activity analysis using a one-tailed Wilcoxon signed-rank test.
Variables
- TF: transcription factor motif name;
- Pval-test_1hrVS0hr: p-value of the Wilconxon signed-rang test, 1 hr stimulation vs unstimulated;
- Pval-test_6hrVS0hr: p-value of the Wilconxon signed-rang test, 6 hr stimulation vs unstimulated;
- Pval-test_6hrVS1hr: p-value of the Wilconxon signed-rang test, 6 hr stimulation vs 1 hr stimulation;
- FDR-test_1hrVS0hr: BH-adjusted p-value (FDR) of the Wilconxon signed-rang test, 1 hr stimulation vs unstimulated;
- FDR-test_6hrVS0hr: BH-adjusted p-value (FDR) of the Wilconxon signed-rang test, 6 hr stimulation vs unstimulated;
- FDR-test_6hrVS1hr: BH-adjusted p-value (FDR) of the Wilconxon signed-rang test, 6 hr stimulation vs 1 hr stimulation.
File: Table_S12_ABC_EnhancerPredictions_threshold0.021_self_promoter.xlsx
Description: Table S12. Enhancer-gene pairs predicted by the ABC (Activity-by-Contact) model, with an ABC score ≥ 0.021. Predictions from all nine contexts are included.
Variables
- Chr: chromosome of the predicted enhancer;
- Start: start position of the predicted enhancer;
- End: end position of the predicted enhancer;
- Name: Name of the predicted enhancer;
- TargetGene: predicted target gene of the corresponding enhancer region;
- TargetGeneTSS : genomic distance (bp) between the corresponding enhancer region and the TSS site of its target gene;
- Context: time x cell type combination of the corresponding enhancer-gene pairs found;
- ABC.Score: the final score of Activity-by-Contact model.
File: Table_S25.Dynamic_caQTLs_for_each_cell_type_with_permutation(n_1000)_testing_result.xlsx
Description: Table S25. Dynamic caQTLs for each cell type with permutation (n=1000) testing result.
Variables
- phenotype_id: phenotype id (as in cis-mod);
- variant_id: SNP id;
- empirical p_val: raw empirical p-value, as calculated by TensorQTL;
- empirical fdr: Adjusted empirical p-values (q-values), as calculated by the qvalue package.
File: Table_S23._cPeak_and_caQTL_identified_by_tensorQTL_in_each_cell_type_and_time_point.xlsx
Description: Table S23. cPeak and caQTL identified by tensorQTL in each cell type and time point.
Variables
- phenotype_id: phenotype id (as in cis-mod);
- num_var: Number of variants in cis-window;
- beta_shape1: Parameter of the fitted Beta distribution;
- beta_shape2: Parameter of the fitted Beta distribution;
- true_df: Degrees of freedom used to compute p-values;
- pval_true_df: Nominal p-value based on true_df;
- variant_id: SNP id;
- start_distance: Distance between the variant and phenotype start position (e.g., TSS);
- end_distance: Distance between the variant and phenotype end position (only present if different from start position);
- ma_samples: Number of samples carrying at least on minor allele;
- ma_count: Number of minor alleles;
- af: In-sample ALT allele frequency of the variant;
- pval_nominal: Nominal p-value of the association between the phenotype and variant;
- slope: Regression slope;
- slope_se: Standard error of the regression slope;
- pval_perm: Empirical p-value from permutations;
- pval_beta: Beta-approximated empirical p-value;
- qval: Storey q-value corresponding to pval_beta;
- pval_nominal_threshold: Nominal p-value threshold for significant associations with the phenotype.
File: Table_S19._Signficant_eQTLs_in_each_cell_type_and_time_point.xlsx
Description: Table S19. Significant eQTLs in each cell type and time point.
Variables
- Snps: the SNP used for eQTL analysis (Matrix-eQTL);
- Gene: the gene used for eQTL analysis;
- Statistic: sample correlation (|r|) as calculated by the simple linear model;
- Pvalue: raw p-value as calculated by the simple linear model;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg);
- Beta: coefficient calculated from the simple linear model;
- Se: standard error calculated from the simple linear model.
File: Table_S24._ASoC_SNPs_in_each_cell_type_and_time_point.xlsx
Description: Table S24. ASoC SNPs in each cell type and time point.
Variables
- CHROM: SNP chromosome;
- POS: SNP position;
- ID: SNP rsid;
- REF: reference allele;
- ALT: alternative allele;
- QUAL: quality score;
- NCALLED: number of total samples collated in analysis;
- REF_N: reference allele count;
- ALT_N: alternative allele count;
- DP: depth;
- pVal: nominal p-value of the binominal test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
File: Table_S6._Differentially_accessible_peaks_upon_KCI_stimulation_table.xlsx
Description: Table S6. Differentially accessible peaks upon KCI stimulation in each cell type.
Variables
- Chr: chromosome;
- Start: peak start position;
- End: peak end position;
- Strand: strandness information.
- Ave-log2FC: average log2 fold change of the corresponding peak between groups;
- P-value : P-value associated with the t-test used during analysis;
- Adj_p: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
File: Table_S14._Gene_regulatory_networks_identified_in_three_cell_types.xlsx
Description: Table S14. Gene regulatory networks identified in three cell types.
Variables
- TF: transcript factor;
- Target Gene: TF-targeting genes;
- Correlation: Spearman’s correlation coefficient.
File: Table_S20._Dynamic_cGenes_in_each_cell_type_and_time_point.xlsx
Description: Table S20. Dynamic eGenes in each cell type and time point.
Variables
- Snps: the SNP used for dynamic eQTL analysis (CellRegMap);
- Gene: the gene used for dynamic eQTL analysis (CellRegMap);
- p_inter: raw p-value from the interaction (GxC) test;
- Bonf: p-value adjusted at the gene level using the Bonferroni procedure to control the family‐wise error rate;
- FDR: the Bonf values further adjusted across genes, Storey’s method.
File: Table_S21._eQTL-based_cTWAS_results.xlsx
Description: Table S21. eQTL-based cTWAS results.
Variables
- PIP_sum: sum of posterior inclusion probability (PIP);
- PIP_0hr_GABA: PIP from 0hr in GABA, the other 8 columns are similar;
- GWAS_Pval_loci: the minimum P value in the loci from GWAS;
- eQTL_Pval: P value of the eQTL from the group with the highest PIP;
- Zgene_Pval: Imputed Z score of the gene with the highest PIP;
- Cell type: the cell type of the gene if over 60% of PIP comes from this cell type;
- PIP_diff : difference of PIP between the resting state (0hr) and the stimulated states (1/6 hour) across all cell types;
- Stimulated: a gene is dynamic if PIP_diff > 50%, otherwise static”;
- Rsid: rsid for the SNPs supporting the risk gene;
- Position: genomic position of the supporting SNP in hg38;
- ref_allele: reference allele;
- eff_allele: alternative allele;
- context: the context from which the SNP supports the eGene.
File: Table_S30._caQTL-based_cTWAS_results.xlsx
Description: Table S30. caQTL-based cTWAS results.
Variables
- condition: the condition with the highest PIP;
- snp: the SNP from the condition with the highest PIP;
- eqtl.gene: potential target genes from colocalized eQTL;
- range: location of the risk peak;
- p2g.gene: potential target gene derived from the correlation between peaks and genes from sc-multi omics profile.;
- overlap.gene: genes that overlap with the peak;
- tss.gene: potential target genes with transcription start site (TSS) closest to the peak;
- tss.dist: distance between the risk peak and the TSS;
- ABC.gene: potential target genes of the risk peaks according to our ABC scores;
- ABC.max: the target gene with the maximum ABC score;
- rsid: rsid for the SNPs supporting the risk gene;
- position: genomic position of the supporting SNP in hg38;
- ref_allele: reference allele;
- eff_allele: alternative allele;
- context: the condition from which the SNP supports the eGene.
File: All_other_supplementary_Tables_(1-3__5__7__8_10_11__15__16-18_22_26-29_31_32_33).xlsx
Description: All_other_supplementary_Tables_(1-3__5__7__8_10_11__15__16-18_22_26-29_31_32_33)
Table S1. Donor information of all the iPSC lines and their summary statistics of sn-multiomics data for all three time points of KCI stimulation.
Variables
- Cell line: Cell line ID representing anonymized individual donors;
- Aff: schizophrenia patients/healthy individuals. Case: SCZ patients; control: healthy individuals;
- Sex: individual sex at birth. M: male; F: female;
- Age: age at blood draw, numerical;
- Ancestry: ancestry information of individual donors;
- Co-culture batch: culture group ID used during neuron differentiation. 3-5 individuals per co-culture batch;
- Seq_batch: group ID describing the batch of Illumina NextSeq sequencing. Factorial;
- Time_point: property describing the stimulation status of each library. 0hr: unstimulated; 1hr: 1 hour post-stimulation; 6hr: 6 hours post-stimulation;
- GABA_count: the number of GABAergic neurons (barcodes) in this sample;
- Nmglut_count: the number of NEFM- glutamatergic neurons in this sample;
- Npglut_count: the number of NEFM+ glutamatergic neurons in this sample;
- intermediateGlut_count: the number of intermediate neurons in this sample;
- unidentified: the number of cells (barcodes) with uncertain identity;
- GABA_proportion: the fraction of GABAergic neurons in this sample;
- Nmglut_proportion: the fraction of NEFM- glutamatergic neurons in this sample;
- Npglut_proportion: the fraction of NEFM+ glutamatergic neurons in this sample;
- intermediateGlut_proportion: the fraction of intermediate neurons in this sample;
- unidentified_proportion: the fraction of cells (barcodes) with uncertain identity;
- total_umber: the gross number of valid cells (barcodes) identified in this sample;
- Used in SCZ analysis: whether this sample has been used during the differential expression analysis between SCZ/control groups.
**Table S2. Summary statistics of each sequencing library before and after QC. Each sequencing library is for cells of a co-cultured 2-4 iPSC lines at a particular time point of stimulation. **
Variables
- Sample ID: unique library IDs assigned to the corresponding sample (as separate 10x Genomics scGEX+ATAC-seq libraries), constructed as [Sequencing batch]_[sequencing group id]-[time];
- Estimated number of cells: QC-passed cell counts (as barcodes) of the corresponding library from the output of 10x Genomics Cell Ranger ARC;
- ATAC Sequenced read pairs: QC-passed scATAC-seq counts of the corresponding library;
- ATAC mean raw read pairs per cell: Mean of QC-passed scATAC-seq counts per cell of the corresponding library;
- ATAC Fraction of high-quality fragments overlapping TSS: the fraction of QC-passed reads overlapping transcription start sites (TSS) of the corresponding library;
- GEX Sequenced read pairs: QC-passed scRNA-seq (Gene expression, GEX) counts of the corresponding library;
- GEX mean raw read pairs per cell: The Mean values of QC-passed scRNA-seq (GEX) counts per cell of the corresponding library;
- GEX Median UMI counts per cell: The median values of QC-passed scRNA-seq (GEX) counts per cell of the corresponding library;
- GEX Reads mapped confidently to genome: gross number of GEX read pairs mapped confidently to human GRCh38 genome of the corresponding library;
- GEX Reads mapped confidently to transcriptome: gross number of GEX read pairs mapped confidently to human GRCh38 transcriptomeof the corresponding library;
- After_demux: cell counts after genotype-based de-multiplex. Representing cells (barcodes) confidently assigned to individuals;
- After_final_QC: Counts of cells (barcode) that are both confidently assigned to individuals and passed QC control. Representing the barcodes used in all subsequent analysis;
- Note: Note.
Table S3. SnRNA-seq post-QC summary data for all 100 lines.
Variables
- Line: cell line, same as in Table S1;
- Cell_counts: gross number of cells that have been assigned to this cell line;
- n_UMI: gross number of Unique Molecular Identifier (UMIs) that have been assigned to this cell line;
- UMI.per.cell: Mean values of UMIs per cell (barcode) of the corresponding cell line;
Table S5. NPD gene sets used for enrichment analyses.
Variables
- Prioritized PGC3 GWAS genes (n=65) (Nature 2022)
- SZ SCHEMA rare risk genes_Q meta <0.05 (Nature 2022)
- ASD risk genes (N=102) (Cell 2020);
- ASD genes (N=185) (Nature Genetics 2022);
- NDD genes (n=664) (Nature Genetics 2022);
- BP GWAS genes (Nature Genetics 2021);
- MDD GWAS genes (Nature Neuroscience 2019);
- PTSD GWAS genes (Nature Genetics 2021).
Table S7. Genes included in each gene module (cluster).
Variables
- Gene: gene symbol (GENCODE v35);
- Cluster: cluster number assigned to the corresponding gene.
Table S8. GO terms (biological process) enriched for each gene module.
Variables
- Term: GO term and ID;
- Overlap: Gene hits, sample/background;
- P.value: computed p-value of the Fisher’s exact test;
- Adjusted.P.value: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg);
- Odds.Ratio: The ratio of the proportion of genes in the query gene set annotated to the corresponding GO term divided by the proportion of genes in the background gene set of the corresponding GO term;
- Combined.Score: calculated as ln(p)*z. P represents the p-value computed using the Fisher exact test. Z represents the z-score computed by assessing the deviation from the expected rank.
- Genes: the list of all gene hits of the corresponding GO term according to the query gene set;
- Cluster: functional cluster assigned to the corresponding GO term.
Table S10. Micro-C summary statistics.
Variables
- Lib: Micro-c library used;
- Total: gross read count and (as the percentage of) total reads of each library;
- total_unmapped: gross unmapped read count and (as the percentage of) total reads;
- total_single_sided_mapped: gross count of read pairs that have only one read mapped to the genome and (as the percentage of) total reads;
- total_mapped: gross count of read pairs that mapped to the genome and (as the percentage of) total reads;
- total_dups: gross count of duplicated read pairs in the library and (as the percentage of) total reads;
- total_nodups: gross count of read pairs after de-duplication and (as the percentage of) total reads;
- cis: gross count of intra-chromosome interactions found and (as the percentage of) total reads;
- trans: gross count of inter-chromosome interactions found and (as the percentage of) total reads;
- cis_1kb+: gross count of intra-chromosome interactions found with separation of more than 1 kb distance and (as the percentage of) total reads;
- cis_2kb+: gross count of intra-chromosome interactions found with separation of more than 2 kb distance and (as the percentage of) total reads;
- cis_4kb+: gross count of intra-chromosome interactions found with separation of more than 4 kb distance and (as the percentage of) total reads;
- cis_10kb+: gross count of intra-chromosome interactions found with separation of more than 10 kb distance and (as the percentage of) total reads;
- cis_20kb+: gross count of intra-chromosome interactions found with separation of more than 20 kb distance and (as the percentage of) total reads;
- cis_40kb+: gross count of intra-chromosome interactions found with separation of more than 40 kb distance and (as the percentage of) total reads;
Table S11. Micro-C chromatin contacts in each sample. Listed are 5 kb bins.
Variables
- bait_chr: chromosome of the bait sequence;
- bait_start: start position of the bait sequence;
- bait_end: end position of the bait sequence;
- bait_name: name of the bait sequence;
- otherEnd_chr: chromosome of the other end of the fragment;
- otherEnd_start: start position of the other end of the fragment;
- otherEnd_end: end position of the other end of the fragment;
- otherEnd_name: name of the fragment, if any;
- N_reads: gross count of the fragments found;
- Score: the final CHiCAGO score of the corresponding fragment.
Table S15. Target genes regulated by the ASD risk TFs in each cell type.
Variables
- TF: transcript factor;
- Cell type: cell type in which the target genes were identified;
- Number of target genes: the number of genes targeted by the corresponding TF;
- Target genes: the full list of target genes.
Table S16. GO-term (Biological Process) enrichment for four ASD risk-associated TFs. The FDR was recalculated to account for multiple testing across the four TFs.
Variables
TF: transcription factor;
- Term: GO term and ID;
- Overlap: Gene hits, sample/background;
- P.value: computed p-value of the Fisher’s exact test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg);
- Odds.Ratio: The ratio of the proportion of genes in the query gene set annotated to the corresponding GO term divided by the proportion of genes in the background gene set of the corresponding GO term;
- Combined.Score: calculated as ln(p)*z. P represents the p-value computed using the Fisher exact test. Z represents the z-score computed by assessing the deviation from the expected rank.
- Genes: the list of all gene hits of the corresponding GO term according to the query gene set.
Table S17. GO-term (Molecular Function) enrichment of target genes shared by at least three of the ASD risk TFs. Only "molecular function" terms are enriched.
Variables
TF: transcription factor;
- Term: GO term and ID;
- Overlap: Gene hits, sample/background;
- P.value: computed p-value of the Fisher’s exact test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg);
- Odds.Ratio: The ratio of the proportion of genes in the query gene set annotated to the corresponding GO term divided by the proportion of genes in the background gene set of the corresponding GO term;
- Combined.Score: calculated as ln(p)*z. P represents the p-value computed using the Fisher exact test. Z represents the z-score computed by assessing the deviation from the expected rank.
- Genes: the list of all gene hits of the corresponding GO term according to the query gene set.
Table S18. TFs with target genes enriched for ASD genes in each cell type.
Variables
- TF: transcription factor;
- Overlap: Gene hits, sample/background;
- Odds.Ratio: The ratio of the proportion of genes in the query gene set annotated to the corresponding GO term divided by the proportion of genes in the background gene set of the corresponding GO term;
- P-value: computed p-value of the Fisher’s exact test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
Table S22. GO term enrichment of cTWAS NPD genes.
Variables
- Term: GO term and ID;
- Overlap: Gene hits, sample/background;
- P.value: computed p-value of the Fisher’s exact test;
- Adjusted.P.value: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg);
- Odds.Ratio: The ratio of the proportion of genes in the query gene set annotated to the corresponding GO term divided by the proportion of genes in the background gene set of the corresponding GO term;
- Genes: the list of all gene hits of the corresponding GO term according to the query gene set.
Table S26. Percentage of ASoC SNPs that can be assigned to a Micro-C based cis-target gene in a matching context.
Variables
- Type: cell type;
- ASoC: ASoC set found in the stimulation state of the corresponding cell types;
- Micro-C: micro-c interaction found in stimulation time of the corresponding cell types (categorised as 0 hr and 1/6 hr and represented unstimulated and stimulated states);
- Interacting bins/with targets: the count of interaction pairs of the corresponding cell type and stimulation stage combinations;
- % of ASoC: interacting ASoCs as the percentage of total ASoCs at the specific cell type and stimulation stage combination.
Table S27. ASoC status of Schizophrenia (SCZ) GWAS risk SNPs and their LD proxies (R2>0.8) in each cell type and time point. Brain eQTL and Micro-C target annotations are from column BN to BT.
Variables
- GWAS_index: rsid of the proxy SNP;
- ASOC_snp: SNP rsid;
- NCALLED: number of total samples collated in analysis;
- REF_N: reference allele count;
- ALT_N: alternative allele count;
- DP: depth;
- pVal: nominal p-value of the binominal test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
- Brain_eQTL: rsid of identified brain eQTL SNPs;
- CommonMind/GTEx eQTL: rsid of identified CommonMind/GTEx eQTL SNPs;
- GTEx eQTL: rsid of identified GTEx eQTL SNPs;
- iPSC neuron Eqtl: rsid of identified iPSC-derived neuron eQTL SNPs;
- Micro_C: target gene from Micro-C analysis;
- Time: stimulation stages where the corresponding SNP was found significant;
- Target: potential targeted genes, if any.
Table S28. ASoC status of Bipolar disorder (BP) GWAS risk SNPs and their LD proxies (R2>0.8) in each cell type and time point. Brain eQTL and Micro-C target annotations are from column BN to BT.
Variables
- GWAS_index: rsid of the proxy SNP;
- ASOC_snp: SNP rsid;
- NCALLED: number of total samples collated in analysis;
- REF_N: reference allele count;
- ALT_N: alternative allele count;
- DP: depth;
- pVal: nominal p-value of the binominal test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
- Brain_eQTL: rsid of identified brain eQTL SNPs;
- CommonMind/GTEx eQTL: rsid of identified CommonMind/GTEx eQTL SNPs;
- GTEx eQTL: rsid of identified GTEx eQTL SNPs;
- iPSC neuron Eqtl: rsid of identified iPSC-derived neuron eQTL SNPs;
- Micro_C: target gene from Micro-C analysis;
- Time: stimulation stages where the corresponding SNP was found significant;
- Target: potential targeted genes, if any.
Table S29. ASoC status of Major depression disorder (MDD) GWAS risk SNPs and their LD proxies (R2>0.8) in each cell type and time point. Brain eQTL and Micro-C target annotations are from column BN to BT.
Variables
- GWAS_index: rsid of the proxy SNP;
- ASOC_snp: SNP rsid;
- NCALLED: number of total samples collated in analysis;
- REF_N: reference allele count;
- ALT_N: alternative allele count;
- DP: depth;
- pVal: nominal p-value of the binominal test;
- FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg).
- Brain_eQTL: rsid of identified brain eQTL SNPs;
- CommonMind/GTEx eQTL: rsid of identified CommonMind/GTEx eQTL SNPs;
- GTEx eQTL: rsid of identified GTEx eQTL SNPs;
- iPSC neuron Eqtl: rsid of identified iPSC-derived neuron eQTL SNPs;
- Micro_C: target gene from Micro-C analysis;
- Time: stimulation stages where the corresponding SNP was found significant;
- Target: potential targeted genes, if any.
Table S31. Result of single cell DE analysis in neurons of each context (Cell type x time point after stimulation) between 28 SCZ cases and matched controls. BH_FDR values were derived from the MAST p-value for each context.
Variables
- geneIDs: gene symbols (GENCODE v35);
- log2FC: log2 fold change of the corresponding gene between groups;
- pvalue: p-values associated with the t-statistic of the corresponding gene;
- padj: Adjusted p-values, corrected for multiple testing (Bonferroni method);
- pct.case: percentage of cells expressing the corresponding gene in SCZ case group;
- pct.ctrl: percentage of cells expressing the corresponding gene in SCZ control group;
- BH_FDR: Adjusted p-values, corrected for multiple testing (Benjamini-Hochberg method);
- Category: the cell types and stimulation stages where the SCZ DEG analysis is performed;
- C6_Cholesterol genes: is the corresponding gene found in the C6 cholesterol gene list;
- PGC3_SCZ genes: is the corresponding gene found in the PGC3 SCZ gene list;
- caQTL_cTWAS gene: is the corresponding gene found as a member of the caQTL-cTWAS list;
- SZ SCHEMA rare risk genes: is the corresponding gene found in the SCZ SCHEMA rare risk genes list;
- ASD_102 genes: is the corresponding gene found in the autism spectrum disorder (ASD) risk gene list.
Table S32. sgRNA sequences and PCR primer sequences as well as qPCR assays.
Variables
- sgRNA sequences, PCR primer sequences and Assay IDs of ThermoFisher qPCR assays.
Table S33. GWAS datasets used in enrichment tests (MAGMA, TORUS, sLDSC and cTWAS).
Variables
- Trait: GWAS disease trait;
- Id: database name in abstract form;
- Link: original publication;
- download link: download link of the corresponding GWAS datasets.
Access information
Access information
Data was derived from the following sources:
- NCBI SRA PRJNA1194194
- NCBI GEO GSE286488
- Zenodo repository https://zenodo.org/records/14577991
Human iPSC lines and cell culture
We initially started with 107 iPSC lines (European ancestry) of which 100 lines were successfully differentiated into both excitatory and inhibitory neurons (Table S1). Of the 100 iPSC lines used for data production, 58 with their IDs starting with “CD” were reprogrammed at Rutgers University Cell and DNA Repository (RUCDR)-NIMH Stem Cell Center using the cryopreserved lymphocytes (CPLs) of donors of Molecular Genetics of Schizophrenia (MGS) cohort, and have been used in previous studies(62, 63, 79-81). The rest were purchased from the California Institute of Regenerative Medicine (CIRM). 28 iPSC lines are from SCZ cases and 72 are from healthy controls. The donors’ average age is 52 (+/-17) years old, and 56 of them are males. All iPSC lines were verified for their identify based on the matching of snRNA/ATAC-seq-inferred genotypes with their known genotypes. Cells were maintained in feeder-free mTeSR Plus (Stemcell# 100-0276) and passaged using ReLeSR (Stemcell# 100-0483) every 4-6 days following the vendor’s instructions.
Lentivirus preparation
Lentiviral particles were prepared using low-passage 293T cells maintained in DMEM media supplemented with 10% FBS. 24 hrs prior to transfection, 95% confluent 293T cells were dissociated using Accutase and replated at 1:3 ratio to achieve 70-80% confluence the next day. On the day of transfection, lentiviral plasmids including FUW-rtTA (Addgene# 20342), TetO-Ascl1-puro (Addgene# 97329), Dlx2-hygro (Addgene# 97330), pTet-O-Ngn2-puro (Addgene# 50247) were co-transfected with pMDLg/pRRE (Addgene #12251), pMD2.G (Addgene #12259) and pRSV-Rev (Addgene #12253) at 1:1:1:1 molar ratio using FuGENE HD (Promega). 18 hrs post transfection, DMEM with 10% FBS media was replaced with fresh mTeSR Plus media. 48 hrs post transfection, supernatant containing viral particles was collected and cell debris were removed by centrifugation at 500 g for 5 min. Supernatant was aliquoted into low-binding tubes and stored at -80C.
Neuronal differentiation and co-culture
We differentiated iPSC lines into both excitatory and inhibitory neurons and co-cultured with rat glial cells. For excitatory neuronal differentiation, we used the method for deriving Ngn2-induced glutamatergic neurons(34) with minor modifications. On DIV (days in vitro) 0, 60-80% confluent iPSCs were dissociated into single cells using Accutase and replated into Matrigel-coated 6-well plate at 5 105 cells per well in mTeSR Plus media with 5M ROCK inhibitor, together with appropriate amount of rtTA virus and Ngn2-puro virus. On DIV1, media was refreshed with mTeSR Plus containing 5M ROCK inhibitor and 2g/ml doxycycline. On DIV2 to DIV4, cells were treated with neural culture media (Neuralbasal media with 1 B27Plus and 1 Glutamax) supplemented with 2g/ml puromycin, 2g/ml doxycycline to remove non-transduced cells. On DIV5 to DIV6, puromycin was withdrawn and cells were treated with neural culture media with 2g/ml doxycycline. On DIV7, cells were ready for replating.
For generating excitatory neurons, we used the method for deriving Ascl1/Dlx2-GABAergic neurons(35) with minor modifications. On DIV0, 60-80% confluent iPSCs were dissociated into single cells using Accutase and replated into Matrigel-coated 6-well plate at 7 105 cells per well in mTeSR Plus media with 5M ROCK inhibitor together with rtTA virus, Ascl1-puro virus and Dlx2-hygro virus. On DIV1, media was refreshed with mTeSR Plus containing 5M ROCK inhibitor and 2g/ml doxycycline. From DIV2 to DIV4, cells were treated with neural culture media supplemented with 2g/ml puromycin, 150g/ml hygromycin and 2g/ml doxycycline to remove non-transduced cells. On DIV5 to DIV6, cells were treated with neural culture media with 2g/ml doxycycline and 2M AraC. Cells were ready for replating on DIV7.
To co-culture the excitatory and inhibitory neurons, on DIV7, separately cultured Ngn2-glutamatergic neurons, Ascl1/Dlx2-GABAergic neurons, and primary rat cortical astrocytes (Thermofisher; N7745100) were dissociated using Accutase at 37C for 15-20 min. Ngn2-glutamatergic neurons, Ascl1/Dlx2-GABAergic neurons and astrocytes were replated at 5:5:1 ratio onto Matrigel pre-coated 12-well plate in neural culture media supplemented with 2 g/ml doxycycline, 10 ng/ml BDNF, 10 ng/ml GDNF, 10 ng/ml NT-3, and 1% FBS. For each co-culture, we pooled cells from 3-4 donor lines with equal proportion. On DIV8, we refilled cells with more media. From DIV9 to DIV33, media was refreshed every three days with half volume changes. Doxycycline was withdrawn at DIV14 and 1 µM AraC was added in the media from DIV8 to DIV17 to ensure neuron purity. DIV33 neurons were used for KCI stimulation.
KCI stimulation of neuronal cultures
To model neuronal activation, we followed a previously used protocol to treat the cells with KCI.(32) Briefly, one day before stimulation, old media was aspirated and DIV33 neuron co-culture were silenced overnight in neural culture media supplemented with 1 uM TTX and 100uM DL-AP5. The next day, we added 0.45 volume (e.g., 0.45 ml for 1 ml original culture media) of warmed depolarization solution (10 mM HEPES, 170 mM KCl, 1 mM MgCl2, 2 mM CaCl2) to initiate KCI stimulation. We prepared three conditions with different treatment durations (0 hr, 1 hr and 6 hrs). After stimulation, cells were dissociated for processing for 10 Genomics sequencing library preparation.
Single nuclei multiomics sequencing library preparation
After stimulation, cells were briefly washed with 1 PBS and dissociated in Accutase at 37C for 40 min with gentle shaking. After Accutase incubation, we pipetted cells multiple times and filtered suspension twice with 40m tip filter (Sigma# BAH136800040-50EA) to obtain single cells. The washed cells were subjected to fresh nuclei isolation following 10 Genomics’ protocols (CG000365 and CG000338). For each library, around 15,000 freshly isolated nuclei were loaded for GEM capture targeting 10,000 nuclei recovery. The pre-amplified DNAs (for snATAC-seq) and cDNAs (for snRNA-seq) were shipped to the University of Minnesota Genomics Center (UMGC) for sn-multiomics sequencing library construction following 10x Genomics’ standard protocol.
Immunocytochemistry
For Immunocytochemistry, iPSCs/neurons were fixed with 4% PFA in PBS at room temperature for 15 min. After three brief washes in PBS, cells were permeabilized with 0.5% Triton X-100 in PBS for 15 min at room temperature and further blocked with 3% BSA and 0.1% Triton X-100 in PBS at room temperature for 1 hr or 4ºC overnight. After blocking, samples were incubated with primary antibodies diluted with blocking buffer at room temperature for 1 hr or 4ºC overnight. After 3 PBS washes, samples were incubated with secondary antibodies diluted in blocking buffer at room temperature for 1 hr followed by 3 more PBS washes. Then samples were incubated in PBS containing 1 g/ml DAPI (4', 6-diamidino-2-phenylindole) at room temperature for 10 min. After DAPI staining, samples were washed once with PBS and mounted on glass slides. The images were taken by a Nikon ECLIPSE C2 confocal microscope.
BDNF OCR peak deletion by CRISPR/Cas9 editing
For BDNF OCR peak deletion, two gRNAs flanking the targeted region were cloned into pSpCas9(BB)-2A-Puro (PX459) V2.0 (Addgene# 62988) respectively. Two iPSC lines were used for editing. For editing, 24 hrs prior to transfection, 60-80% confluent iPSCs were dissociated into single cells using Accutase and replated into Matrigel-coated 60 mm dish at 4.5 105 cells per well in mTeSR Plus media with 5M ROCK inhibitor. On the day of transfection, 3 g of plasmid DNA carrying gRNA1 and 3 g of plasmid DNA carrying gRNA2 were introduced into iPSCs using LipofectamineSTEM (Thermofisher) following vendor’s instruction at 1:2 DNA:reagent ratio. 24-48 hrs post transfection, cells were selected with 0.5 g/ml puromycin; 48-72 hrs post transfection, cells were selected with 0.25 g/ml puromycin. Afterwards, antibiotics were withdrawn, and cells were maintained in regular mTeSR Plus until colonies reached appropriate size for picking. After colony picking, genomic DNAs (gDNAs) from collected cell pellets was isolated using QuickExtract DNA Extraction Solution for PCR amplification. Amplified DNAs were loaded on 1% agarose gel for electrophoresis to examine fragment size. gDNAs of the colonies with confirmed peak deletion were used for Sanger sequencing confirmation. The confirmed colonies (2-3) with the expected peak deletion were sub-cloned and expanded for cryopreservation. Please refer to Table S32 for gRNA and primer sequences.
RNA isolation and qPCR
Total RNA was isolated using RNeasy Plus Kits (QIAGEN) following vendor instructions. For qPCR, 300 ng-1 g RNAs were reverse transcribed using High-Capacity cDNA Reverse Transcription Kit (Thermofisher) following vendor instruction. cDNAs were diluted with nuclease-free water at 1:10 ratio and qPCR reactions were prepared using Taqman Universal PCR Master Mix (Thermofisher). Reactions were loaded on Roche LightCycler 480 system in a 384-well white plate. The 2–∆∆Ct method was used for RNA expression quantification, with GAPDH as endogenous control. Please refer to Table S32 for qPCR assay information.
Multiomics sequencing and data quality control (QC)
10 Genomics Chromium single cell Multiome sequencing libraries were sequenced at UMGC on Novaseq S4 platform targeting pair-end (2 150bp) 50 K reads per nuclei for ATAC library and 25 K reads per nuclei for gene expression library. Briefly, after raw data collection, Illumina’s BCL2FASTQ software was used to demultiplex and assemble the fastq files corresponding to reads (R1/R2) and indices (i5/i7) for Cell Ranger ARC. The fastq files were subsequently processed by 10 Genomics Cell Ranger ARC (v2.0.2) and aligned twice to both the human GRCh38.p14 genome and a contingent of human GRCh38.p14 and mouse GRCm38 provided by 10 Genomics for efficient identification and removal of rodent astrocytes. For snATAC-seq, the per-library Transcription Start Sites (TSS) enrichment score is usually 7-9, and the Fraction of high-quality Reads overlapping Peaks (FRiP) is usually 40-80%. We excluded cells with >15% of the reads mapped to mitochondrial genes, cells with > 8,000 or < 400 number of features, and cells with > 40,000 or <500 number of UMI counts (filtered by 400 < nCount < 8,000 using Seurat).
Multiomic data analyses
Barcode level identification of cell line identity
For each library that contained cells derived from 3-4 iPSC lines, barcode-level identification was performed separately for snRNA-seq (GEX) and snATAC-seq data with genotyping information of all donors. Briefly, the BAM files generated from GEX and ATAC assays (gex_possorted_bam.bam and atac_possorted_bam.bam) were processed by demuxlet(82) with known genotype information provided as .vcf files. Barcodes (cells) identified as singlets with P1 likelihood (p1LLK) < 1 10-8 were collected, and only barcodes positively identified in both GEX and ATAC assays were retained for downstream analysis.
snRNA-seq data analyses
snRNA-seq (GEX) data were processed by extracting gene expression data from the .h5 files generated by Cell Ranger ARC. Only barcodes (cells) confidently assigned to individual iPSC lines were retained. With Seurat 5.1.0(83), To assign cell type identity, Leiden clustering was performed at the library level based on the first 30 PCs (or appropriate as determined by the elbow plot) and subsequently comparing cell clusters and their marker gene expressions (GAD1, GAD2, SLC17A6, SLC17A7, NEFM). The library-level gene expression matrices were subsequently collated in Seurat 5.1.0(83) as one large object with cell line-specific metadata assigned. Harmony(84) was used to integrate the library-level gene matrices and removed sequencing batch-derived effects.
snATAC-seq data analyses
For snATAC-seq data, the aggregation function of Cell Ranger ARC was used to merge and generate new fragment files for MACS2-based peak calling with the CallPeaks() function in Signac(85) with default settings. We made two merged fragment files, one containing all libraries from batch 024 (18 cell lines) and the other from batch 018, 024, 025, and 029 (76 cell lines). From each fragment file, we iterated through the combinations of cell types (NEFM+ glut, NEFM- glut, GABA) and stimulation stages (0 hr, 1 hr, 6 hrs), which generated nine peak sets (cell type time). We further generated a union peak set from batch 024 by setting CallPeaks(combine.peaks = TRUE) for peak analysis across different cell types. Finally, peaks that fell within the ENCODE blacklisted regions were removed. We also removed peaks that fell within chromosomes X and Y and the mitochondrial genome.
Approximation of developmental stages for iPSC-derived neurons
We used the reference single-cell dataset of human embryonic and postmortem brains(38) to evaluate our iPSC-derived neurons and their comparable developmental stages. Briefly, the excitatory neuron (ExNeu) and interneuron (IN) populations were extracted using the identity assigned in the original publications. The excised data were then processed as documented in the original publication with Harmony-assisted data normalization, scaling, and dimension reduction using the first 30 PCs to generate the anchor gene set, UMAP dimensions, and unimodal UMAP coordinates for projection. The three main cell types identified (npglut, nmglut, GABA) from snRNA-seq were used as queries. Each cell type was firstly randomly subset to 10,000 cells to reduce computational demands. Subsequently, each cell subset was separately projected by the MapQuery() function to its corresponding cell map (npglut/nmglut to ExNeu and GABA to IN) to show their approximate developmental stages in the human brain.
PCA analysis
For PCA analysis, we generated pseudo-bulk count matrices from merged snRNA-seq data using the AggregatedExpression() function from Seurat, and each sample represented one cell line in one of three main cell types and its stimulation stage. To further eliminate the interference from sequencing batch-specific effects, we ran regression using the ComBat-seq() function from the R package sva(86) using sequencing batch information as the factor. The sva-adjusted pseudo-bulk count matrices were log-transformed and estimated for their observation-level weights using the voom() function from the R package limma(87). Finally, we calculated the principal components (PCs) with the prcomp() function and plotted the samples using their first two PCs using the fviz_pca_ind() function from the R package factoextra (https://cran.r-project.org/package=factoextra).
Pseudo-bulk RNA-seq differential gene expression analysis
Gene differential expression (DE) analysis was performed using the R package limma. Briefly, the ComBat-seq-adjusted pseudo-bulk count matrices generated from PCA analysis were used as the initial input. Low-expression genes whose Count Per Million (CPM) value was less than < 1 in half of lines in any of the 3 time points were removed. The matrices were log-transformed and observation-level weights were calculated by voom(). Subsequently, we made the design matrix of known covariates (time, batch, age, sex, SCZ-affection status, and the percentage of corresponding cell types). We fitted the linear model using the lmFit() function of limma. After checking the mean-variance trend (SA plot) for gene distributions, the DE values (log2FC, p-value, and FDR) were calculated using the topTable() function with corresponding coefficients (1 hr vs 0 hr, 6 hrs vs 0 hr, 6 hrs vs. 1 hr) for all three cell types.
MAGMA gene set enrichment analysis
We performed MAGMA analysis using MAGMA version 1.08b(88) to evaluate the enrichment for the GWAS risk of several psychiatric disorders (SCZ, Neuroticism, ASD, BP, and MDD)(7, 10, 89-92), as well as the GWAS set of Crohn’s disease that served as a control set(93) (Table S33). The method was adapted from our previous publications with modifications(63). Specifically, we compiled the MAGMA-required gene annotation data files based on the more recent GRCh38/hg38 genome to get the gene-SNP annotation file. With the gene-SNP annotation file, we then performed gene-level analysis on SNP p-values using the reference SNP data of 1,000 Genomes European panel (g1000_eur_hg38, --bfile) and the pre-computed SNP p-values from each disorder’s GWAS data set. The sample size (ncol=) was derived from either GWAS data frames or specified according to the affiliated README data. Subsequently, the result files (--gene-results) from the gene-level analysis were read in for competitive gene-set analysis (--set-annot), where we used default setting (correct=all) to control for gene sizes in the number of SNPs and the gene density (a measure of within-gene L.). The gene-set analysis produced the output files (.gsa.out) with competitive gene-set analysis results that contained the effect size (BETA) and the statistical significance of the enrichment of each gene set (upregulated/downregulated) for each disorder’s GWAS data set.
Gene ontology enrichment analysis
GO term analysis was performed for the DEG sets specific to each cell type and expression pattern. Briefly, the list of DE genes was used as the SynGo(94) input and all the expressed genes in the corresponding cell type were used as the background gene list. A similar approach was used to calculate the enrichment of different gene expression patterns against a multitude of known NPD gene sets.
Differential peak accessibility analysis
DA peak accessibility was performed using aggregated pseudo-bulk counts based on the three cell-type-specific peak intervals (nmglut, npglut, GABA). Briefly, the AggregatedExpression() function was used to generate cell line peaks matrices using the cell-type-specific peak intervals (approximately 300 K peaks per cell type). The ComBat_seq() function from sva was subsequently applied to correct group bias for each count matrix. We removed all low-count peaks (CPM < 1 in at least half of the cell lines), and approximately 200 K of the original 300 K peaks survived filtration. The matrices were log-transformed and observation-level weights were calculated by voom(). Subsequently, we applied design matrices of known covariates (time, batch, age, sex, SCZ-affection status, and the percentage of corresponding cell types) as we did in DE gene analysis. We fitted the linear model using the lmFit() function of limma. The DA values (log2FC, p-value, and FDR) of peaks for each cell type and contrast were calculated using the topTable() function with corresponding coefficients (1 hr vs 0 hr, 6 hrs vs 0 hr).
Stratified linkage disequilibrium score regression (sLDSC) for GWAS enrichment analysis
sLDSC(95) analysis was performed by using the hg38 version of European genotype data (SNPs) from 1000 Genomes Phase 3 and v2.2 baseline linkage disequilibrium/weights. Briefly, linkage disequilibrium score estimations were pre-calculated from the hg38 version of the 1000 Genomes EUR file set (w_hm3_no_hla.snplist), window size 1 cM (ld-wind-cm 1). We used the summary statistics of major psychiatric disorders and non-psychiatric diseases for partition heritability, with several data sets lifted over from hg19 to hg38 when necessary. Disease-specific regressions were performed independently using hm3 SNP weights against each disease for cell-type-specific analysis. Integration of snATAC-Seq Profiles and Peak Calling for caQTL and gene regulatory network analyses To integrate the snATAC-seq data, we utilized ArchR. Quality control measures were applied by filtering cells with a Transcription Start Site (TSS) enrichment score of at least 4 and fragment counts of at least 1,000. After quality control and integration of expression data, we retained a total of 552,653 cells. To call peaks from the snATAC-Seq data, pseudo-bulk replicates were generated for each cell type and for each time point. To ensure equal representation of all cell types, we followed the default setting of ArchR, and sampled a maximum of 500 cells and 25 million fragments from each pseudo-bulk replicate. With the pseudo-bulk replicates, 501-bp fixed-width peaks were called using MACS2(96) and integrated in ArchR. Peaks were required to be reproducible in at least two samples, and sex chromosomes were excluded in the analysis. To ensure high-quality peak detection, we limited the number of peaks to a maximum of 250,000 for each cell type and time point context. Consensus peaks were obtained through an iterative overlap peak merging approach in ArchR. This method allowed for robust and reproducible peak identification across different samples and conditions.
Single-cell Embeddings
To reduce computational cost, we downsampled 18 cell lines (115,855 cells in total) from Batch 024 for all analyses related to gene regulation, including dimensionality reduction, chromVAR analysis, pseudo-time inference, peak-gene mapping, and gene regulatory network construction. Using ArchR, we performed dimensionality reduction on snRNA-Seq and snATAC-Seq data separately, utilizing the iterative latent semantic indexing approach. To create a more robust two-dimensional representation using UMAP, we combined the embeddings of expression and peak accessibility, subsequently correcting for batch effects using Harmony.
Inference of Pseudo-time
Pseudo-time trajectories were estimated for the GABA, nmglut, and npglut cell types separately to model the transitions between neuron states from unstimulated neurons to neurons stimulated for 1 hr and 6 hrs. Based on the UMAP constructed in the previous section, cell-type-specific trajectories were inferred using the addTrajectory() function in ArchR. We fine-tuned the sparsity of the spline fit and quality filters for the trajectory of each cell type to avoid overfitting of smoothing spline. Cells were divided into 100 bins based on pseudo-time, and the average activities were calculated by pooling cells with similar pseudo-time values.
Co-activation analysis to link OCRs to genes
To link OCRs to potential target genes, we assessed the correlation of gene expression and normalized chromatin accessibility for all ORCs within 500 kb of genes. All differentially expressed genes were included in the analysis. For each cell type, peaks present in at least 5% of the cells were considered. Gene expression was regressed on peak accessibility using a negative binomial mixed regression implemented using lme4 R package (DOI: https://doi.org/10.18637/jss.v067.i01) with the following formula:
Expression_g ~ Accessibility_p+CellLine_ID + log(Gex_LibrarySize)+log(Gex_MitoRatio)
where:
- Expression_g represents the raw counts of gene expression for gene g,
- Accessibility_p denotes the normalized chromatin accessibility of peak p,
- CellLine_ID is a categorical variable representing cell line identity, included as a random intercept to account for donor-specific effects,
- Gex_LibrarySize refers to the single-cell level library size from the snRNA-Seq data, captured by the number of reads per cell,
- Gex_MitoRatio captures the single-cell level percentage of mitochondrial reads from the snRNA-Seq experiment, serving as a measure control for cell quality.
To control for multiple testing, the Benjamini–Hochberg procedure was applied to calculate the false discovery rate (FDR). OCR-gene pairs were considered significant and biologically relevant if they had an FDR ≤ 0.05 and a positive regression coefficient.
Activity-by-Contact (ABC) analysis
We employed the Activity-by-Contact (ABC) model(46, 97) to predict enhancer-gene regulatory relationships. For our analysis, we utilized snATAC-Seq data from all 18 samples in batch 024. Additionally, we incorporated time-specific bulk Micro-C data derived from one cell line (CD-11), averaged across all cell types, which allowed us to capture the dynamic nature of chromatin interactions during different neuronal stages. Our analysis aimed to predict enhancer-gene pairs across nine distinct contexts, defined by different combinations of time and cell type. To achieve this, we called peaks for each specific context and subsequently predicted ABC scores under these conditions. An enhancer-gene pair was classified as significant if its ABC score was greater than or equal to 0.021, a threshold recommended by the software to achieve 70% recall.
Gene Module Analysis
To cluster genes based on their expression trajectory in pseudo-time, we first selected 5,221 differentially expressed and highly variable genes, using criteria of FDR ≤ 0.05 and abs (LogFC) > 1 in at least one differential expression test. We then normalized expression for each gene separately, so that genes with similar trajectory patterns would be clustered together, irrespective of their absolute expression levels. Using the average expression values from 100 pseudo-time points across three cell types, we performed K-means clustering with 15 cluster centers to group the genes into 15 distinct modules. Gene Ontology enrichment analysis was conducted using the Enrichr package(98-100), and GWAS enrichment analysis was performed in MAGMA. We inferred chromatin accessibility trajectories for the 15 gene modules by mapping all peaks to the genes within each module, as described in the previous section.
Defining candidate TFs regulating early and late responses
We used chromVAR to estimate single-cell level enrichment of TF activity. With motif annotations curated from Cis-BP, we estimated chromVAR motif deviations implemented in the ArchR toolkit. We also had single-cell expression from scRNA-seq for all TFs. For each TF, we assessed its role in early and late responses by testing differential motif activity and expression between conditions, one cell type at a time. Specifically, a TF was considered a candidate TF for early response if it showed both higher motif activity and increased expression at 1 hr compared to 0 hr. Similarly, a TF was considered a candidate TF for late response if it showed higher motif activity and elevated expression at 6 hrs compared to both 1 hr and 0 hr. Differential expression analysis was performed using the limma-voom analysis, and differential motif analysis was performed using one-tail Wilcoxon signed-rank test using single-cell level motif deviations estimated by chromVAR.
Gene Regulatory Network Analysis
To construct the gene regulatory network, we followed these steps: (1) Identify Potential TF Regulators: We identified a set of potential TF regulators by selecting TFs that were differentially expressed in at least one condition. We also ensured that a TF’s motif activity (estimated by chromVAR) and its expression activity were positively correlated. Specifically, we filtered TFs based on the correlation between their expression trajectories and motif activity trajectories, retaining only those with a Spearman’s correlation > 0.3. (2) Check Motif Presence in Open Chromatin Regions (OCRs): We examined the presence of the TF motifs in OCRs linked to each gene, using an FDR threshold of ≤ 0.1. (3) Correlate TF Motif Activity with Target Gene Expression: For motifs present in at least one peak associated with a gene, we correlated the motif activity of TFs with the expression of the target genes. We retained the TF-gene pairs that exhibited a Spearman’s correlation > 0.5 between the motif activity of the TF and the expression of the target gene. This methodology allows us to identify key regulatory relationships between TFs and their target genes, revealing insights into the underlying gene regulatory networks the cellular transitions observed in the pseudo-time trajectories.
eQTL mapping
To associate the genotypic variation and gene expression variation within each context (combination of cell types and time points), we used the scRNA read count matrix of 36,601 features and 548,800 cells from 100 cell lines and three cell types (i.e., GABA, nmglut, npglut). We first generated pseudo-bulk matrices by summing up the read counts by each feature for cells within a certain time point, a certain cell type, and a certain cell line. We removed samples with less than 1 million total reads. We also removed genes that were not protein-coding, in mitochondria, with less than 2,000 reads across all cells, or without HGNC symbol. After filtering, we retained the pseudo-bulk read count matrix of 14,818 genes by 824 samples. Then, we performed trimmed mean of M values (TMM) normalization and inverse Normal Transformation. The resulting gene expression matrix served as input for MatrixeQTL, along with covariates of cell type composition, gender, age, SCZ affection status, 5 genotype PCs, and 15 gene expression PCs. For genotypes, we removed SNP with minor allele frequency (MAF) below 5%. 5,966,820 SNPs remained for eQTL testing. We used cis-variants 250 kb up- and downstream of the gene body. This cis window size and the number of gene expression PCs were chosen to optimize the number of eGenes (i.e., genes with at least one significant eQTL) discovered. Our scheme for defining significant eQTL is detailed by the standard procedure of MatrixeQTL(101).
In addition to the eQTL mapping in each context, we also performed eQTL mapping by aggregating all contexts. We removed the 5 cell lines from batch 22 and summed up the read counts per cell line per contexts. TPM normalization was applied to the summed reads and then we further summed up all contexts within the same cell line. We proceeded to eQTL mapping with the resulting 95 samples (cell lines) by applying TMM, INT, and MatrixeQTL. This is the same procedure as described above, except that the number of expression PCs was 21 instead of 15. eQTL results from this approach are called “pseudo-bulk eQTL” in subsequent sections.
Dynamic eQTL testing
We tested the genetic effect heterogeneity across 3 times points, one cell type at a time. We used the following procedure to choose the candidate set of gene-SNP pairs to be tested for each cell type. We considered only eGenes found in at least one time point. Genes that were only significant in pseudo-bulk eQTL were not considered. For each of these eGenes, we started with its eQTL(s) with the smallest nominal P values per condition. Then we took a union of these top eQTLs along with the pseudo-bulk eQTL. There are at most four eQTLs for each eGene in a cell type. We then performed LD pruning with r2<0.1 on these eQTLs within the same eGene.
For the candidate gene-SNP pairs, we used CellRegMap (102) to test the heterogeneity of genetic effects across different contexts (time points in our study) from eQTLs. We have n samples (95 cell lines in our case), and each sample belongs to one of K contexts (3 in our case). We use Z (n×K matrix) as the design matrix, encoding whether a sample belongs to a certain context. Under the CellRegMap model, the genetic effect of a variant on expression of a gene has two components: an average effect denoted as β, and a context-specific effect, denoted as u_k for the k-th context. The model assumes that the context-specific effects follow a normal distribution, u_k~N(0,σ^2). To test if the genetic effects vary across context means testing if σ=0.
We write the regression model in a matrix form, equivalent to the assumption described above:
y = Gβ + Gγ + Wψ + μ + ϵ
where y is gene expression jointly normalized across three time points for each cell type (dimension is n = 95x3), G is genotype, W are covariates, μ is the individual random effect, ϵ is the iid error term, and β and ψ are fixed effects corresponding to genotypes and covariates, respectively. The focus of the model is γ, the random slope for genotypes. γ is a vector with n dimensions, following multivariate Normal distribution:γ ∼ MVN(0,σ^2 ZZ^T)
CellRegMap follows StructLMM's idea of using Rao's score test to evaluate the null hypothesis that σ is 0. According to CellRegMap, this is the highest power test when many contexts were jointly tested. We used Bonferroni to adjust nominal P values in each gene, retain eQTLs with the minimum adjusted P-value per gene, and then applied Benhamini-Hochberg FDR control to all the remaining eQTL.
TORUS enrichment analyses of eQTLs and ASoC variants
We applied a Bayesian hierarchical model (TORUS) to perform SNP-based enrichment analysis, testing the enrichment of certain features with eQTL and ASoC variants (62, 103). Briefly, TORUS associates the prior probability that a variant is causal to gene expression with the annotations of the variant with a logistic regression model. TORUS then uses this prior in a simple fine-mapping model (assuming a single causal variant per genetic locus) of eQTL data. TORUS uses Maximum Likelihood to estimate the parameters of the logistic regression. In dynamic eQTL enrichment analysis, we used the upregulated peaks in the matched cell type. To assess the enrichment of ASoC in eQTL, we applied the ASoC in matched condition as the neuronal activity eQTL. ASoC SNPs were treated as binary annotations: 1 being within the significant ASoC and 0 otherwise. We used SNPs that did not pass the ASoC testing (FDR>5%) for comparison. For testing disease GWAS enrichment, TORUS assumes that every variant is a risk variant or not, represented by a binary indicator variable (1 or 0). The prior probability of the indicator of a SNP being 1 depends on its annotations. Here, we tested ASoC SNPs categorised by their cell type time against a selection of diseases/traits GWAS databases. All the annotations are encoded as binary (1 if an SNP is found in the corresponding GWAS database, 0 otherwise). We performed univariate analysis in TORUS to assess the enrichment of each ASoC cell type time combinations. The GWAS datasets used for enrichment/TORUS analysis were from multiple sources, including NPDs, neuro-related traits, and control disorders/traits, as listed in Table S33. Analysis of genetic effect sharing between neuronal activity eQTL and GTEx and ASoC. We used Storey’s π1 analysis to analyze the sharing between our eQTL and other summary statistics. This is a common approach to investigate sharing of eQTL under different cell types. For each context in our eQTL and each tissue in GTEx brain eQTL, we obtained all eQTL (FDR<0.05) in a context and estimated the proportion of non-null tests (π1, Pi1) based on the binomial p values of these gene-SNP pairs in a GTEx brain tissue. A similar approach was taken when we computed π1 for our dynamic eQTL in GTEx brain tissues. We also used Storey’s π1 analysis to investigate the sharing between our eQTL and ASoCs. For each ASoC SNP with FDR<0.05, we located its nearest gene and obtained the nominal P values of the corresponding eQTL in the matched context if possible. Then, we estimated non-null proportions of these P values.
caQTL Mapping
We calculated pseudo-bulk counts from ATAC-seq for 94 GABAergic neuron cell lines and 95 Glutamatergic cell lines under each time point. The peak insertion counts were initially normalized for library size using the trimmed mean of M values (TMM) method. Subsequently, we applied a rank-based inverse normal transformation (INT) to standardize the accessibility of each peak to a standard normal distribution, utilizing the RNOmni package(104). These counts served as the input for caQTL analysis. We employed TensorQTL (DOI: https://doi.org/10.1186/s13059-019-1836-7) for caQTL mapping. To adjust for batch effects and confounding factors, we included three genotype principal components and between four to nine chromatin accessibility PCs as covariates in our analysis:
Accessibility_p ~ Genotype_i + P〖C^G〗_1+ ...+P〖C^A〗_n
where Accessibility_p denotes the normalized chromatin accessibility of peak p, Genotype_i genotype of a cis-SNP i, P〖C^G〗_n are genotype principle components, and P〖C^A〗_n are the chromatin accessibility principle components,
For each peak, we conducted 1,000 permutations and computed beta-approximated P-values of the top QTL. To account for multiple testing, we calculated Q-values using the qvalue package (DOI: 10.18129/B9.bioc.qvalue). Peaks with a Q-value ≤ 0.05 were considered significant and designated as cPeaks (caQTL peaks). This threshold controlled the FDR and ensured the reliability of our findings. To investigate time-dependent QTL effects for each cell type, we included the union of cPeaks from all time points. This inclusion ensured that all cPeaks from different time points were considered. For each cPeak, we identified and tested the SNP from the strongest QTL across all time points. This SNP served as the representative variant for assessing QTL effects. Similar to caQTL mapping, we first normalized the library size using counts per million (CPM), then applied rank-based INT to each peak. To test time-dependent effects, we modified the caQTL test above, including Time, and Genotype-Time interaction term. We included 3 genotype PCs, 5 chromatin accessibility PCs, and their interaction effects with time as covariates and ran linear regression as:
Accessibility_p ~ Genotype_i+Time_t +Genotype_i*Time_t + P〖C^G〗_1+P〖C^G〗_1*Time_t+ ...+P〖C^A〗_5+P〖C^A〗_5*Time_t
where Time_t is a categorical variable of the duration of KCl stimulation.
We employed an ANOVA test to compare the full model (including the time-genotype interaction effects) with a reduced model (excluding the time-dependent effect). This allowed us to assess the significance of the time-dependent effect on QTLs. FDR correction was applied to adjust for multiple comparisons, ensuring the robustness of our results. To validate the test statistics and avoid inflation, we performed 1,000 permutations. This permutation approach verified that the observed test statistics were not artificially inflated, thereby enhancing the reliability of our time-dependent QTL effect analysis.
ASoC mapping
We applied an improved, two-step ASoC mapping with WASP to calibrate for alternate-allele bias effect and to gain more accurate results. BAM files generated from Cell Ranger ARC (atac_possorted_bam.bam) were firstly split and re-merged by cell line using barcodes from initial demultiplexing. Each line-specific BAM file included aligned reads of all cell types and three stimulation stages. GATK (4.2.6.0) HaplotypeCaller was applied to each BAM file (cell line), which identified and extracted all heterozygotic SNP sites and homozygotic SNP sites carrying alternative alleles. The extracted SNP coordination was subsequently used to generate the positional reference of WASP. Next, the original BAM files from Cell Ranger ARC were re-split at the cell line cell type stimulation stage (time) levels to generate raw BAM files for WASP calibration. For a set of BAM files from the same line, WASP calibration was performed at the individual BAM files using the line-specific SNP heterozygosity information to re-align reads carrying alternative alleles with bwa and remove any bad read pairs. A second round of GATK HaplotypeCaller was used to identify heterozygous SNP sites and their counts in each calibrated BAM file to generate corresponding VCF files as recommended by the GATK Best Practice (software.broadinstitute.org/Gatk/best-practices/) with VariantRecalibrator (-an DP -an Q.D. -an F.S. -an SOR -an M.Q. -an ReadPosRankSum -mode SNP -tranche 100.0 -tranche 99.5 -tranche 95.0 -tranche 90.0) and reference databases including HapMap v3.3 (priority = 15), 1000G_omni v2.5 (priority = 12), Broad Institute 1000G high confidence SNP list phase 1 (priority = 10), Mills 1000G golden standard INDEL list (priority = 12), and dbSNP v154 (priority = 2). Heterozygous sites with tranche level >95.0% were extracted. And only SNPs with corresponding rs# records found in dbSNP v154 were retained. The individual VCF files were subsequently merged at the cell type stimulation stage level to increase the statistical power (each merged output included 96-100 lines as cell type varied). After merging, Biallelic SNP sites (GT: 0/1) with minimum read depth count (DP) ≥ 30 and minimum reference or alternative allele count ≥ 2 were retained. The binomial p-values (non-hyperbolic) were calculated using the binom.test(x, n, P = 0.5, alternative = “two.sided”, conf.level = 0.95) from the R package, and Benjamini & Hochberg correction was applied to all qualified SNPs as the correcting factor of R function p.adjust(x, method = "fdr"). We set the threshold of ASoC SNP at FDR value = 0.05.
Homer enrichment analysis of TF-binding motifs for ASoC SNPs
We used HOMER(105) to assess the enrichment of TF binding motifs in sequences flanking ASoC SNPs (25 bp) adapted from our previous investigations(62). We used the JASPAR database (2022 release) with all 522 human TF motifs. For each cell type time, we used the ASoC SNPs (FDR < 0.05) as input intervals. The parameters used were findMotifsGenome.pl <input.bed> hg38 -cpg -mknown jasper2018.known <output>. HOMER provided the genetic background for the enrichment test, and the enrichment significance was derived from HOMER.
Brain eQTL enrichment analysis for ASoC SNPs
We performed eQTL enrichment analysis to investigate the differences in the association between ASoC and non-ASoC SNPs in GTEx cortex-eQTL(91) and PsychENCODE eQTL(106) databases with adaptation(62). For ASoC/non-ASoC SNPs in each set of cell type time combination, the SNPs were assigned to genes (eQTL targets) if the distance of the SNP and its corresponding gene TSS was less than 500 kb, and the SNP was associated with the expression of the corresponding gene at FDR < 0.05.
Micro-C analysis of ASoC SNPs and their targets The Micro-C experiment was conducted in co-cultured Ngn2-Glut and GABA neurons of two donor lines (CD11, CD12) at 0 hr, 1 hr, and 6 hrs post-stimulation. Dovetail™ Micro-C (Pan-promoter Capture) kit was used process the samples. To minimize background interference from cell-line discrepancies, we used the samples derived from a single line (CD11) from comparing different time points. The Micro-C data was performed by an external vendor and analysed using ChiCAGO(107) at 500 bp resolution. To assign the micro-C targets of each SNP, we firstly associated each SNP (0 hr, 1 hr, 6 hrs) with its nearest interacting micro-C fragment (from 0 hr, 1 hr, 6 hrs samples) if their distance was less than 2.5 kb (upstream/downstream). If such an association coould be established, we extracted the other end of the interacting micro-C fragment and checked if it overlapped with the promoter region (-2 kb / +1 kb of TSS). If the other fragment fell into the promoter region, we added the original SNP to the pool of ASoC SNPs with their corresponding identity.
Promoter and enhancer enrichment analysis for ASoC SNPs
We defined the interactions of ASoC SNPs and human enhancers/promoters according to our previously developed algorithms with adaptation(62). The analysis was similar to that implemented in GREAT and was described previously(66, 108, 109) for testing the folds of peak enrichment within the annotated genomic regions and epigenetically annotated regions. To adapt the original algorithm for SNP-based testing, we developed an updated version of the algorithm, and the enrichments were calculated using the formula below:
"ratio" □(=(("total number of SNPs in the cell type time combination overlap the features" /"genome size" ))/(("total lengths of the features" /"genome size" )∙("total SNP number of the cell type time combination" /"genome size" ) ))
For ASoC SNPs that were cell-type-specific or shared by three cell types, we analyzed enrichment considering SNP intervals as the unit similar to OCR peaks. The enrichment p-values were estimated based on the binomial test as implemented in GREAT and as previously described(108, 109).
p_"enrichment" =∑_(i=k_π)^n▒〖(■(n@i)) p_π^i (1-p_π )^(n-i) 〗
from which we calculated the p-value using the
binom.test(x, n, p, alternative=“greater”) function in R, where
x = numbers of SNPs that fell within the designated epigenetically annotated features;
n = total length of the SNP intervals (in bp);
p = total length of the designated epigenetically annotated features (in bp) / genome size (3.2 109 bp)
The Gene-based annotation of the genome was derived from GENCODE v35 as part of the built-in database of the updated HOMER package(105). The list of human forebrain/non-forebrain enhancers was from the VISTA enhancer browser(110). Human-gained enhancers were acquired from GSE63648 using data from either the frontal or occipital cortex at 12 PCW and marked by H3K27ac or H3K4Me2(111). The definitions of chromatin state were assembled using an imputed 25-state model derived from individual #E081 of fetal brain tissue by the Roadmap Epigenomics Project(109, 112). For estimating the enrichment of ASoC SNPs (Fig. 2A), the categories of epigenomic features used for promoter and enhancer annotations were promoter = TssA, PromU, PromD1 and PromD2; enhancers = TxReg, TxEnh5, TxEnh3, TxEnhW, EnhA1, EnhA2, EnhAF, EnhW1, EnhW2, EnhAc, and DNase. The abbreviation of different chromatin states was promulgated below, as used in(109):
TssA Active TSS
PromU Promoter Upstream TSS
PromD1 Promoter Downstream TSS with DNase
PromD2 Promoter Downstream TSS
Tx5’ Transcription 5’
Tx Transcription
Tx3’ Transcription 3’
TxWk Weak transcription
TxReg Transcription Regulatory
TxEnh5’ Transcription 5’ Enhancer
TxEnh3’ Transcription 3’ Enhancer
TxEnhW Transcription Weak Enhancer
EnhA1 Active Enhancer 1
EnhA2 Active Enhancer 2
EnhAF Active Enhancer Flank
EnhW1 Weak Enhancer 1
EnhW2 Weak Enhancer 2
EnhAc Enhancer Acetylation Only
Dnase DNase only
ZNF/Rpts ZNF genes & repeats
Het Heterochromatin
PromP Poised Promoter
PromBiv Bivalent Promoter
ReprPC Repressed PolyComb
Quies Quiescent
Forebrain_enh Forebrain enhancers
Non_forebrain_enh Non-forebrain enhancers
H3K27acF H3K27 acetylated regions found in the frontal lobe
H3K27acO H3K27 acetylated regions found in the occipital lobe
H3K4me2F H3K4me2 regions found in the frontal lobe
H3K4me2O H3K4me2 regions found in the occipital lobe
Integrative analysis of NPD GWAS with neuronal activity eQTL or caQTL
Causal TWAS (cTWAS) (56) was used to analyze the heritability of NPD mediated by the molecular traits from our QTL and ASoCs. Compared to other TWAS approaches, cTWAS employs a regression model that effectively reduces false positives by jointly analyzing all nearby variants and genes. Given the genetically predicted expression of all genes (as in standard TWAS), cTWAS assumes a regression model:
y=∑_j▒〖β_j X ̃_j 〗+∑_m▒〖θ_m G_m 〗+ϵ,ϵ~N(0,σ^2)
where βj and θm are the effect sizes of gene expression j and the variant m, respectively. cTWAS assumes a spike-and-slab prior for the effects of gene expression and variants, respectively. The parameters of these distributions are estimated through an expectation-maximization (EM) algorithm, using all the data from the genome. Once these parameters are estimated, cTWAS solves the regression model above with SuSiE, a fine-mapping method. The results would be Posterior Inclusion Probabilities (PIPs) for all the expression traits and variants. While the model is formulated using individual level data, the method can be used with summary statistics. In that case, a reference LD matrix matching the GWAS samples is used.
The cTWAS model above has only gene expression in a single context. But it is easy to extend the model to allow multiple groups of QTL in the same model, with exactly the same algorithm. For example, we could have eQTL from multiple tissues or cellular contexts; or one eQTL group and one caQTL group, in the model. We used this version of cTWAS in our study. This integration can combine evidence across multiple groups to improve power. For example, if we have eQTLs across several contexts in the model, then we can obtain the PIP of a gene in each of the contexts. These PIPs can be combined to obtain the probability that the gene is causal in at least one context (denoted as gene PIP). Meanwhile, we can also quantify the importance of each context, defined as the ratio of PIP of that context over the total PIP across all contexts. This allows us to assign causal context for a putative causal gene.
Another advantage of cTWAS is that it can estimate (1) the proportion of heritability mediated by each group of molecular traits and (2) fold enrichment of molecular traits compared to background SNPs. These parameters are derived from the prior parameter estimation step:
enrichment=π_G/π_v
mediated h^2 g=PVE_G/(PVE_G+PVE_v)
where π_G is the estimated prior of a group of molecular trait, π_v is the estimated prior of background SNPs, PVE_G is the proportion of variance explained by a group of molecular traits, and PVE_v is the proportion of variance explained by the background SNPs.
When we used cTWAS to integrate eQTL with NPD traits, we first extract top eQTL per gene per context and use these eQTL as the genetic prediction models (denoted as “weights”) of the gene expression traits. cTWAS then includes eQTL data across all 9 contexts in analysis. We note that even though we analyzed all 9 contexts together, a particular gene may occur only in a subset of 9 contexts, as it may not have eQTL in all contexts. However, one problem is that our eQTL study has incomplete power, thus for a particular gene, we may miss some relevant contexts. To address this issue, we used a relaxed threshold. Suppose a gene has eQTL in one context (say C1), but no eQTL in the second context (C2). If the top eQTL in C1 has p value < 0.1 in C2, then this eQTL would be included as the weight in C2. This allows us to include as many contexts as possible for a single gene.
To run cTWAS, we used the default setting, allowing at most 5 causal signals per locus. We used the LD reference of UK Biobank White British samples. For ADHD GWAS, we initially found a large number of putative causal genes, suggesting that LD mismatch between the reference and GWAS samples may be a problem. We thus ran ADHD analysis with the option of L = 1, i.e., allowing at most a single causal signal per locus. Under this setting, the results are independent of the reference LD matrix.
We used the same procedure when applying cTWAS to caQTL. The difference is that we complemented caQTL with ASoC SNPs. In the case where caQTL and ASoC have overlapped peaks, we chose the one with smaller nominal P values. In addition, we used a more stringent threshold of nominal P values at 0.05 instead of 0.10 to increase confidence in the results and reduce the number of molecular traits included.
Single-cell differential gene expression analysis in SCZ neurons
We performed DEG analysis between neurons from SCZ donors and those from controls to identify genes that exhibited SCZ case-control differential neuronal activation. Briefly, we first subset the snRNA-seq data of the 28 SCZ lines and the sex/age-matched 28 controls (17 males and 11 females in SCZ cases with an average age of 48; 16 males and 12 females in controls with an average age of 50). The subset Seurat object was first split by individual x cell type and processed by standard Seurat SCTransform for data transformation. Subsequently, all transformed libraries were re-integrated by Seurat IntegrateLayers using “HarmonyIntegration” method to construct the master Seurat object for MAST analysis using the zlm hurdle model. We used covariates including cell co-culture batch, sex, age, and cell type fraction in MAST test. The Seurat wrapper FindMarkers was used to find DE genes between SCZ case and control groups. Only genes expressed in at least 5% of cells in either comparison group were used. We recalculated FDR for the result of each context (cell type/ timepoint). FDR < 0.05 was used as cut-off for DEGs used for the gene set enrichment analyses.
Statistical analyses
For genetic analyses (eQTL and caQTL mapping and GWAS enrichment analyses), different statistical methods were used and specified in the method section above and in different figure legends. For other assays, unless otherwise specified, Student’s t-test (between two groups) or the Kruskal Wallis test with Dunn’s multiple comparisons and P-value adjustment (more than two groups) was used to determining significance between groups. Samples were assumed to be unpaired and have non-parametric distribution unless otherwise specified. Data were analyzed using R 4.3.2 and GraphPad Prism 10. Results were considered as significant if P < 0.05 (*: P < 0.05; **: P < 0.01; ***: P < 0.001; ****: P < 0.0001). All data are reported as mean ± SEM.