Cell type-specific network analysis in Diversity Outbred mice identifies genes potentially responsible for human bone mineral density GWAS associations

  1. Luke J Dillard
  2. Gina Calabrese
  3. Larry Mesner
  4. Charles Farber  Is a corresponding author
  1. Department of Genome Sciences, University of Virginia, United States
  2. Department of Biochemistry and Molecular Genetics, School of Medicine, University of Virginia, United States
5 figures, 1 table and 4 additional files

Figures

Figure 1 with 1 supplement
Analysis of single-cell RNA-seq (scRNA-seq) data for bone marrow-derived stromal cells cultured under osteogenic conditions (BMSC-OBs) derived from 80 Diversity Outbred (DO).

(A) Uniform Manifold Approximation and Projection (UMAP) of 139,392 single cells (BMSC-OBs). Cell numbers and corresponding percentages for the fifteen (15) annotated cell clusters are listed in parenthesis to the right of the annotated cluster name. (B) Dot plot (Marsh et al., 2023) portraying representative and highly expressed genes for all annotated cell clusters. Dot color indicates the scaled gene expression while the size of the dot corresponds to the percentage of cells of a given cluster that express a given gene. (C) The raw proportional abundances of seven (7) mesenchymal cell clusters and one (1) cluster (Hem) representing the remaining cells (i.e. mainly hematopoietic immune cells) across all 80 DO mice. (D) UMAP plots for mesenchymal lineage cell clusters for DO mouse 50 and DO mouse 233. (E) CELLECT (CELL-type Expression-specific integration for Complex Traits) cell-type prioritization results displaying the Bonferroni adjusted p-values for the cell clusters. The OB1, Ocy, and marrow adipogenic lineage progenitor (MALP) cell clusters (red) were significantly enriched (padj<0.05, red dashed line) for BMD heritability (padj = 0.018, 0.010, 0.006, respectively).

Figure 1—figure supplement 1
Plots displaying the distribution of the total number of cells from each mouse (N=80).

(A) Density plot portraying the distribution of the total number of cells from each mouse after processing of the single-cell RNA-seq (scRNA-seq) data. (B) Boxplot displaying the distribution of the total number of cells for each mouse (Min: 723, 1st Qu: 1316, median: 1690, mean: 1742, 3rd Qu: 2118, Max: 3652). (C) Quantile-quantile plot (Q-Q plot) with 95% confidence interval. Shapiro-Wilk normality test: p-value = 0.1061; W=0.97425.

Figure 2 with 1 supplement
Overview of the network analysis pipeline.

Step 1: For all seven (7) of the mesenchymal lineage cell clusters (mesenchymal progenitor cell [MPC], late mesenchymal progenitors [LMP], osteoblast progenitor [OBP], OB1, OB2, osteocyte-like cell [Ocy], marrow adipogenic lineage progenitor [MALP]), cell type-specific co-expression modules were generated using iterative Weighted Gene Co-expression Network Analysis (iterativeWGCNA). Step 2: Bayesian networks were learned to generate directed networks and model causal interactions between co-expressed genes. Step 3: Differentiation driver genes (DDGs) were identified by extracting subnetworks (i.e. large three-step neighborhood) for each gene in each cell type-specific Bayesian network and highlighting those subnetworks that were enriched (padj<0.05) for trajectory-specific tradeSeq genes for the cell-type boundary. Step 4: DDGs (and associated networks) were prioritized if the DDG was identified previously as an expression/splicing quantitative trait loci (eQTLs/sQTLs) that colocalized with BMD genome-wide association studies (GWAS) associations. Created with Biorender.com.

Figure 2—figure supplement 1
Scale-free topology and mean connectivity graphs for the cell type-specific iterative Weighted Gene Co-expression Network Analysis (iterativeWGCNA).

A soft thresholding power of 14 was selected for the generation of all co-expression modules for all clusters, which was the point at which R2 exceeded a threshold of 0.85.

Pseudotime trajectory inference analysis and establishment of cell-type boundaries for tradeSeq analysis.

(A) Three (3) trajectories (two adipogenic, one adipogenic) were inferred from the mesenchymal cell clusters of the bone marrow-derived stromal cells cultured under osteogenic condition (BMSC-OB) single-cell RNA-seq (scRNA-seq) data using Slingshot. All trajectories originate from the mesenchymal progenitor cell (MPC) and end in either osteogenic (osteocyte-like cells [Ocy], OB2) or adipogenic (marrow adipogenic lineage progenitor [MALP]) cell fates. (B) For each of the trajectories, cell-type boundaries were generated using pseudotime values along the trajectories, which encompass the majority of cells of a cell-type mapping to their respective trajectory. (C) Normalized gene expression of select genes associated with each cluster is represented in feature plots (top) and each gene plotted as a function of pseudotime (bottom) for all pseudotime trajectories (color corresponds to cell-type annotation observed throughout). Vertical lines (red) represent the cell-type (pseudotime) boundaries established for each cell type (label). The cell-type boundary for OB1 and OB2 is represented as one red line/label for visualization purposes.

Figure 4 with 1 supplement
TradeSeq-identified genes associated with bone marrow-derived stromal cells cultured under osteogenic condition (BMSC-OB) differentiation exhibit expression quantitative trait locus (eQTL) effects.

(A) Pkm was identified as a significant global tradeSeq-identified gene (padj = 8.35 × 10–8) for late mesenchymal progenitor (LMP) cells along an osteogenic trajectory (LMP_to_OBP) (left). S100a1 was identified as a significant global tradeSeq-identified gene (padj=0.023) for OBP cells along osteogenic trajectory 1 (OBP_to_OB1) (right). (B) Plots indicating the cell type-specific eQTLs signal for both Pkm and S100a1. A negative eQTL effect on Pkm expression was observed in LMPs for Diversity Outbred (DO) mice with a PWK haplotype background at the Pkm locus (left). A positive eQTL effect on the expression of S100a1 was observed in OBPs for DO mice with a 129 haplotype background at the S100a1 locus, while a negative effect was observed for NZO mice (right). (C) The expression of Pkm and S100a1 based on DO mouse (expression values transformed via variance stabilizing transformation [VST], as described in Methods). Genotype abbreviations correspond to DO haplotype background (legend) at the respective gene locus. Mice with at least one PWK allele (genotype abbreviation G) tend to have decreased expression of Pkm (left). Mice with at least one 129 allele (genotype abbreviation C) tend to have increased expression of S100a1, while NZO mice (genotype abbreviation E) have decreased expression (right). (D) PWK mice had a significant reduction in mature osteoblasts (OB1) and osteocyte-like cells (Ocy) proportions relative to other mice (p=0.030 and p=0.026, respectively; t-test), while LMP proportions were unaffected. Asterisks represent any of the other haplotype backgrounds. 129 mice showed a significant decrease in LMP proportion and increase in OB1 proportion (p=0.008 and p=0.016, respectively; t-test), but OBP proportions were unaffected. No significant effects on cell-type proportions were observed in NZO mice (Figure 4—figure supplement 1).

Figure 4—figure supplement 1
Tests of significance for cell-type proportions for NZO mice.

Mice with at least one NZO allele at the S100a1 locus (N=14) had no significant difference in cell-type proportions (p>0.05; t-test) as compared to mice with other DO haplotype background at this locus. Asterisks represent any of the other haplotype backgrounds.

Fgfrl1 and Tpx2 are prioritized differentiation driver genes (DDGs) and putative drivers of mesenchymal differentiation.

(A) Fgfrl1 was identified as a DDG of a network for late mesenchymal progenitors (LMPs) differentiating along an adipogenic trajectory. The network is enriched (padj = 7.5 × 10–3) for trajectory-specific tradeSeq-identified genes for the LMP_to_MALP cell-type boundary (Hnmt, St5, Igfbp4, Cyp1b1, Pdzrn4, Mark1). Fgfrl1 was previously identified as a gene that has bone mineral density (BMD) genome-wide association studies (GWAS) associations that colocalize with an expression quantitative trait locus (eQTL) in the cultured fibroblast Genotype-Tissue Expression (GTEx) tissue. (B) An increase in the expression of tradeSeq-identified genes coincides with the LMP_to_MALP cell-type boundary in which they were identified as significant. (C) Tpx2 was identified as a DDG of a network for LMPs differentiating along an osteogenic trajectory. The network is enriched (padj = 5.7 × 10–7) for tradeSeq-identified genes for the LMP_to_OBP cell-type boundary (Tpx2, Top2a, Kif4, Iqgap3, Prc1, Kif11, Ect2, Sgo2a, Ube2c). Tpx2 is both a tradeSeq gene and previously identified as a gene that has BMD GWAS associations that colocalize with an eQTL in the Testis GTEx tissue. (D) An increase in the expression of tradeSeq-identified genes coincides with the LMP_to_OBP cell-type boundary in which they were identified as significant. (E) Box plot displaying whole-body BMD measurements (excluding skull) from the International Mouse Knockout Consortium (IMPC) for Tpx2 mutant mice, which exhibited a significant increase in BMD (genotype p-value = 1.03 × 10–3) in both male and female mice (N=8 (M) and 8 (F) mutants; N=2574 (M) and 2633 (F) controls).

Tables

Table 1
Prioritized differentiation driver genes (DDGs) that have bone mineral density (BMD) genome-wide association studies (GWAS) associations that colocalize with splicing/expression QTL (eQTL/sQTL) identified in a Genotype-Tissue Expression (GTEx) project tissue.

The tissue with the most significant colocalization (RCP and/or H4PP) is listed for each DGG (26 total, 21 distinct), as determined from Al-Barghouthi et al., 2022, and Abood et al., 2023, for eQTL and sQTL, respectively (Al-Barghouthi et al., 2022; Abood et al., 2023). RCP=Regional Colocalization Probability (GWAS and eQTL colocalization). H4P=H4 Posterior Probability (GWAS and sQTL colocalization).

TrajectoryCell-type boundaryDDGGTEx Tissue with strongest eQTL colocalization(RCP)GTEx Tissue with strongest sQTL colocalization(H4PP)eGene identified from scRNA-seq of the 80 DO mice
1LMP to OBPTet1Adipose (Visceral);
0.3191
1LMP to OBPTpx2Testis;
0.2031
1LMP to OBPCdk1Pituitary;
0.7795
1LMP to OBPTtyh3Liver;
0.9350
1LMP to OBPOlfml3Artery (aorta);
0.8048
1LMP to OBPIzumo4Brain (hypothalamus);
0.9182
1LMP to OBPSec24dNerve (tibial);
0.2677
1LMP to OBPTmem263Adipose (subcutaneous);
0.5704
Cultured cells (fibroblasts);
0.9716
1LMP to OBPLmf2Adrenal gland;
0.8181
1LMP to OBPTln2Esophagus (muscularis);
0.9697
1OBP to OB1Kremen1Heart (left ventricle);
0.8686
2OBP to OB2Kremen1Heart (left ventricle);
0.8686
2OBP to OB2Ebf1Testis;
0.8760
2OBP to OB2Lrp4Pancreas;
0.7943
3LMP to MALPTtyh3Liver;
0.9350
3LMP to MALPFgfrl1Cultured cells (fibroblasts);
0.1611
3LMP to MALPEbf1Testis;
0.8760
3LMP to MALPPpp1r12bNerve (tibial);
0.8807
3LMP to MALPRhojCultured cells (fibroblasts);
0.352
Breast;
0.7844
3LMP to MALPTln2Esophagus (muscularis);
0.9697
3MALP to endAdh1Esophagus (gastroesophageal junction);
0.9999
3MALP to endFgfrl1Cultured cells (fibroblasts);
0.1611
3MALP to endAdcy5Esophagus (gastroesophageal junction);
0.8456
3MALP to endCnn2Spleen;
0.7743
3MALP to endMxra8Pituitary;
0.7545
3MALP to endTimp2Testis;
0.9429

Additional files

Supplementary file 1

Supplementary file 1a-h.

(a) Differentially expressed genes (DEGs) for all clusters of the bone marrow-derived stromal cells cultured under osteogenic condition (BMSC-OB) scRNA-seq cell clusters. DEGs were calculated on all clusters of the BMSC-OB scRNA-seq data using the FindAllMarkers function from the Seurat R package. (b) DEGs between the OB1 and OB2 clusters of the scRNA-seq data. DEGs were calculated using the FindMarkers function from the Seurat R package. Positive values for average log2 fold change (avg_log2FC) indicate that a gene is more highly expressed in OB1. (c) BMSC-OB cell-type proportion analysis for the 80 Diversity Outbred (DO) mice. The raw proportions (top) and asin-transformed proportions (bottom) of each of the BMSC-OB cell types were calculated from the total number of cells contributed by each mouse using the Propeller R package. All non-mesenchymal lineage cell types (i.e. hematopoietic lineage cells) are aggregated as a group (Hem) for each mouse. (d) Correlation of cell proportions to bone trait metrics captured from the 80 DO mice. Raw (top) and transformed (bottom) cell-type proportions were correlated using Pearson and Spearman to bone trait metrics (55 total) captured on all mice from the 80 DO mice. (e) Bone trait abbreviations and units of measurement. (f) CELLECT cell-type prioritization table. Beta is a regression effect size estimate for given annotation. Beta SE is the standard error for the regression coefficient. The p-value is the one-sided test (beta>0) association between bone mineral density (BMD) genome-wide association study (GWAS) signal heritability and each annotated cell type. p-Values were adjusted using the Bonferroni correction method. MALP=marrow adipogenic lineage precursors; Ocy=osteocyte-like cell; OB1=osteoblast population 1; MPC=mesenchymal progenitor cell; LMP=late mesenchymal progenitor; OBP=osteoblast progenitor; OB2=osteoblast population 2; EC=endothelial cell; MF1=macrophage population 1; MO=monocyte; BC=B-cell; GC=granulocyte; OC=osteoclast-like cell; TC=T-cell; MF2=macrophage population 2. (g) Summary of results from the iterative Weighted Gene Co-expression Network Analysis (iterativeWGCNA). A total of 535 co-expression modules were generated using the mesenchymal lineage cell clusters (7 total) of the BMSC-OB scRNA-seq data, yielding an average of 76 modules per cell cluster. A total of 8810 Bayesian networks were generated from the co-expression modules. (h) Genes within each module generated from the iterativeWGCNA.

https://cdn.elifesciences.org/articles/100832/elife-100832-supp1-v1.xlsx
Supplementary file 2

Supplementary file 2a-k.

(a) Summary of tradeSeq-identified genes. For each cell-type (pseudotime) boundary associated with a specific trajectory (9 total), a global and trajectory-specific test was performed using the startVsEndTest function from the tradeSeq R Package. The number of genes identified for each test and for each boundary is displayed, as well as the number of tradeSeq-identified genes that were also identified as eGenes from the expression quantitative trait locus (eQTL) mapping of the 80 Diversity Outbred (DO) mice (73 total). (b) TradeSeq-identified genes from the trajectory-specific analysis for osteogenic trajectory 1. All significant trajectory-specific tradeSeq-identified genes (padj ≤0.05) across all cell-type boundaries (5 total; MPC, LMP, OBP, OB1, Ocy) associated with osteogenic trajectory 1. Associated eQTL information is also displayed for the gene if it was an eGene identified in the cell type from the cell type-specific eQTL analysis (if ‘NA’ is present, the gene was not identified as an eGene). (c) TradeSeq-identified genes from the trajectory-specific analysis for osteogenic trajectory 2. All significant trajectory-specific tradeSeq-identified genes (padj≤0.05) across all cell-type boundaries (2 total; OBP, OB2) associated with osteogenic trajectory 2. Associated eQTL information is also displayed for the gene if it was an eGene identified in the cell type from the cell type-specific eQTL analysis (if ‘NA’ is present, the gene was not identified as an eGene). (d) TradeSeq-identified genes from the trajectory-specific analysis for the adipogenic trajectory. All significant trajectory-specific tradeSeq-identified genes (padj≤0.05) across all cell-type boundaries (2 total; LMP, MALP) associated with the adipogenic trajectory. Associated eQTL information is also displayed for the gene if it was an eGene identified in the cell type from the cell type-specific eQTL analysis (if ‘NA’ is present, the gene was not identified as an eGene). (e) TradeSeq-identified genes from the global analysis. All significant global tradeSeq-identified genes (padj≤0.05) across all cell type (pseudotime) boundaries (9 total). Associated eQTL information is also displayed for the gene if it was an eGene identified in the cell type from the cell type-specific eQTL analysis (if ‘NA’ is present, the gene was not identified as an eGene). (f) Results from the cell type-specific eQTL analysis on the mesenchymal lineage cell types identified in the scRNA-seq data from the 80 DO mice. All significant eQTLs (LOD>9.68 for autosomal chromosomes; LOD>9.49 for X chromosome) and eGenes for the mesenchymal cell clusters (563 total). Chr=chromosome of eQTL, Pos=eQTL peak position, LOD=logarithm of the odds score, ci (low/hi)=LOD support intervals, Start=start position of gene (GRCm38), End=end position of gene (GRCm38), dist_start=distance of eQTL to start. (g) Results from tests of significance for cell-type proportions. Tests of significance on the transformed cell-type proportions were performed using the Propeller R package, and nominal p-values are reported. Sample batch (pool containing cells from mice in preparation for scRNA-seq) and sex were modeled as covariates. For the Pkm example, a t-test was performed for all mice with at least one PWK haplotype background (PWK/*; asterisk meaning any DO haplotype) at a Pkm locus (n=15) against all remaining mice (n=65). For the S100a1 example, a one-way ANOVA was performed on four groups: mice with at least one 129/* haplotype background (n=26) or NZO/* (n=10), heterozygous for both (129/NOZ, n=4), or any other DO haplotype background at the locus (n=40). Additionally, t-tests were performed on the 129/* and NZO/* haplotype background individually against all other mice aggregated as a group. (h) Summary of differentiation driver gene (DDG) network analysis. The number of DDGs and associated networks that were enriched (padj≤0.05) with more genes in the trajectory-specific tradeSeq genes for each cell-type boundary (408 total). The number of DDGs that had a corresponding human homolog with a human BMD GWAS association (that colocalizes with expression and/or splicing quantitative trait loci [eQTLs/sQTLs] from the Genotype-Tissue Expression (GTEx) Project) is also displayed (26 total, 21 distinct). Three of the DDGs were also tested by the International Mouse Knockout Consortium (IMPC) and had a significant BMD phenotype when knocked out. (i) All significant DDG network analysis for osteogenic trajectory 1 (178 total). The enrichment of each DDG Bayesian network for tradeSeq-identified genes (identified for each cell-type boundary along osteogenic lineage 1) is displayed as nominal and adjusted p-values, as well as the co-expression module in which the DDG was identified. The data can be filtered to highlight DDGs that are: a tradeSeq-identified gene for the cell boundary, a gene that was identified by Al-Barghouthi et al., 2022, as having eQTL that also colocalizes with BMD GWAS associations, a gene that was identified by Abood et al., 2023, as having sQTL that also colocalizes with BMD GWAS associations, a gene that was tested by the IMPC and had a significant effect on BMD when knocked out, or gene that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data. Other columns of the data contain information at the network level: Number_neighbors=number of nodes (genes) in Bayesian network, Number_tradeSeq_gene_neighbors=number of genes that were also tradeSeq-identified genes for the cell-type boundary, All network neighbors=all genes in network, Neighbors_eQTL_Al_Barghouthi=all genes in the network that were identified from Al-Barghouthi et al., 2022, Neighbors_sQTL_Abood=all genes in the network that were identified from Abood et al., 2023, Neighbors_IMPC_BMD_gene=all genes in the network that were tested by the IMPC and had a significant effect on BMD when knocked out, Neighbors_DO_eGene=all genes in the network that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data, Neighbors_tradeSeq_gene_for_boundary=all genes that were also tradeSeq-identified genes for the cell-type boundary. (j) All significant DDG network analysis for osteogenic trajectory 2 (55 total). The enrichment of each DDG Bayesian network for tradeSeq-identified genes (identified for cell-type boundary along osteogenic lineage 2) is displayed as nominal and adjusted p-values, as well as the co-expression module in which the DDG was identified. The data can be filtered to highlight DDGs that are: a tradeSeq-identified gene for the cell boundary, a gene that was identified by Al-Barghouthi et al., 2022, as having eQTL that also colocalizes with BMD GWAS associations, a gene that was identified by Abood et al., 2023, as having sQTL that also colocalizes with BMD GWAS associations, a gene that was tested by the IMPC and had a significant effect on BMD when knocked out, or gene that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data. Other columns of the data contain information at the network level: Number_neighbors=number of nodes (genes) in Bayesian network, Number_tradeSeq_gene_neighbors=number of genes that were also tradeSeq-identified genes for the cell-type boundary, All network neighbors=all genes in network, Neighbors_eQTL_Al_Barghouthi=all genes in the network that were identified from Al-Barghouthi et al., 2022, Neighbors_sQTL_Abood=all genes in the network that were identified from Abood et al., 2023, Neighbors_IMPC_BMD_gene=all genes in the network that were tested by the IMPC and had a significant effect on BMD when knocked out, Neighbors_DO_eGene=all genes in the network that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data, Neighbors_tradeSeq_gene_for_boundary=all genes that were also tradeSeq-identified genes for the cell-type boundary. (k) All significant DDG network analysis for the adipogenic trajectory (175 total). The enrichment of each DDG Bayesian network for tradeSeq-identified genes (identified for cell-type boundary along the adipogenic trajectory) is displayed as nominal and adjusted p-values, as well as the co-expression module in which the DDG was identified. The data can be filtered to highlight DDGs that are: a tradeSeq-identified gene for the cell boundary, a gene that was identified by Al-Barghouthi et al., 2022, as having eQTL that also colocalizes with BMD GWAS associations, a gene that was identified by Abood et al., 2023, as having sQTL that also colocalizes with BMD GWAS associations, a gene that was tested by the IMPC and had a significant effect on BMD when knocked out, or gene that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data. Other columns of the data contain information at the network level: Number_neighbors=number of nodes (genes) in Bayesian network, Number_tradeSeq_gene_neighbors=number of genes that were also tradeSeq-identified genes for the cell-type boundary, All network neighbors=all genes in network, Neighbors_eQTL_Al_Barghouthi=all genes in the network that were identified from Al-Barghouthi et al., 2022, Neighbors_sQTL_Abood=all genes in the network that were identified from Abood et al., 2023, Neighbors_IMPC_BMD_gene=all genes in the network that were tested by the IMPC and had a significant effect on BMD when knocked out, Neighbors_DO_eGene=all genes in the network that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data, Neighbors_tradeSeq_gene_for_boundary=all genes that were also tradeSeq-identified genes for the cell-type boundary.

https://cdn.elifesciences.org/articles/100832/elife-100832-supp2-v1.xlsx
Supplementary file 3

Supplementary file 3a-g.

(a) PANTHER Gene Ontology (GO) Enrichment analysis for differentiation driver genes (DDGs) identified for the late mesenchymal progenitor (LMP) to osteoblast progenitor (OBP) cell-type boundary (osteogenic trajectory 1). (b) PANTHER GO Enrichment analysis for DDGs identified for the OBP to OB1 cell-type boundary (osteogenic trajectory 1). (c) PANTHER GO Enrichment analysis for DDGs identified for the OBP to OB2 cell-type boundary (osteogenic trajectory 2). (d) PANTHER GO Enrichment analysis for DDGs identified for the LMP to MALP cell-type boundary (adipogenic trajectory). (e) PANTHER GO Enrichment analysis for DDGs identified for the MALP to the end (of the trajectory) cell-type boundary (adipogenic trajectory). (f) Prioritized DDG network analysis for the adipogenic trajectory (26 total, 21 distinct). The enrichment of each prioritized DDG Bayesian network for tradeSeq-identified genes (identified for the cell-type boundary along the associated trajectory) is displayed as nominal and adjusted p-values, as well as the co-expression module in which the DDG was identified. The data can be filtered to highlight DDGs that are: a tradeSeq-identified gene for the cell boundary, a gene that was identified by Al-Barghouthi et al., 2022, as having eQTL that also colocalizes with BMD GWAS associations, a gene that was identified by Abood et al., 2023, as having sQTL that also colocalizes with BMD GWAS associations, a gene that was tested by the IMPC and had a significant effect on BMD when knocked out, or gene that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data. Other columns of the data contain information at the network level: Number_neighbors=number of nodes (genes) in Bayesian network, Number_tradeSeq_gene_neighbors=number of genes that were also tradeSeq-identified genes for the cell-type boundary, All network neighbors=all genes in network, Neighbors_eQTL_Al_Barghouthi=all genes in the network that were identified from Al-Barghouthi et al., 2022, Neighbors_sQTL_Abood=all genes in the network that were identified from Abood et al., 2023, Neighbors_IMPC_BMD_gene=all genes in the network that were tested by the IMPC and had a significant effect on BMD when knocked out, Neighbors_DO_eGene=all genes in the network that was identified here as an eGene in the cell type-specific eQTL analysis using the 80 DO scRNA-seq data, Neighbors_tradeSeq_gene_for_boundary=all genes that were also tradeSeq-identified genes for the cell-type boundary. (g) Expression specificity scores (ESµ) for each gene for each cell cluster of the BMSC-OB scRNA-seq data for the 80 DO mice. ESµ scores are generated during the CELLEX portion of the CELLECT analysis pipeline. ESµ values range from 0 to 1 and are a gene’s marginal likelihood of being specifically expressed in a given cell type, where 1 is the most specific and 0 is not specific.

https://cdn.elifesciences.org/articles/100832/elife-100832-supp3-v1.xlsx
MDAR checklist
https://cdn.elifesciences.org/articles/100832/elife-100832-mdarchecklist1-v1.docx

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Luke J Dillard
  2. Gina Calabrese
  3. Larry Mesner
  4. Charles Farber
(2026)
Cell type-specific network analysis in Diversity Outbred mice identifies genes potentially responsible for human bone mineral density GWAS associations
eLife 13:RP100832.
https://doi.org/10.7554/eLife.100832.3