从Rec 96/085的选举终止获得了发展的人类肢体和颅骨组织样本，并获得了所有样本捐赠者（英格兰东部的书面和知情同意） ，并得到了剑桥中央研究伦理委员会的全面批准。简而言之
，在解剖过程中将样品悬浮在PBS中，在冰上-4°C 。从四肢解剖肩膀 ，臀部和膝关节
。对于肩关节
，在锁骨的远端进行了近端切口，并在肱骨的脖子上产生了远端切口。对于未形成独特骨特征的胚胎肩部样品 ，进行了近似值以捕获全藻藻和杂技关节 。对于颅样品（少于8个PCW）
，在两种情况下在额骨的后部边界分离的每个瓦尔瓦里亚和颅底都剖析了两个区域
。对于较旧的颅样品（超过8个PCW），剖析组织以分离可行的额叶 ，顶叶，脾
，筛，胸膜和颞骨。最初 ，将样品嵌入最佳切割温度培养基中
，并在异义烷干燥的冰浆中冷冻在-80°C下 。然后，使用Leica CM1950低温恒温器以10μm的厚度切除冷冻切片 ，并将其放置在ISS或森林中的Superfrost Plus载玻片（VWR）上
，或直接用于单核处理
。对于整个安装免疫染色的样品，通过从所有样本供体获得的书面和知情同意书中获得了妊娠终止的样品。样品由Inserm的Hudeca生物库提供 ，并符合法国法规 。法国生物医学研究机构（Agence de laBioMédecine
，Saint-Denis la Plaine，France; PFS19-012）和Inserm Ethics委员会（IRB00003888）授予了完全授权使用这些组织。 　　使用10倍基因组学Cartana HS库制备套件（1110-02 ，遵循协议D025）和原位测序试剂盒（3110-02
，遵循协议D100）进行ISS，该版本构成了Hybriss63的商业化版本
。探针面板的设计基于四肢细胞状态中的倍数变化阈值（补充表3）。简而言之 ，发育肢体的冷冻切除术在PBS中的3.7％甲醛（252549，默克）中固定30分钟
，然后在PBS中两次洗涤1分钟 。在0.1 m HCl（10325710，Fisher）37°C中以0.5 mg ml -1胃蛋白（P7012 ，默克）在37°C中短暂地消化切片
，持续15 s（5 pCW）或30 s（6个PCW及以上），然后在PBS中再次在PBS中洗涤两次
。脱水在70％和100％乙醇中 ，每次脱水2分钟后
，将直径为9毫米（50μl）固定杂交室（GBL621505-20EA，GRACE BIO-LABS）粘附在每个载玻片上 ，并用于保存后续反应混合物。在缓冲液WB3中补液后
，在37°C下进行缓冲液RM1中的探针杂交 。158-PLEX探针面板包括每个基因的5个挂锁探针，其序列是专有的（10x基因组Cartana）
。将截面用PBS-T（具有0.05％TWEEN-20）的PBS洗涤两次 ，然后在37°C下用缓冲液WB4进行30分钟
，再用PBS-T再次洗涤3次。在37°C下进行RM2的探针连接2小时，并用PBS-T将截面洗涤3次 ，然后在30°C下进行RM3中的滚动圆扩增18小时 。在PBS-T洗涤后，将所有滚动圆产物（RCP）与LM（与DAPI的Cy5-Cy5贴合混合物）杂交30分钟
，在室温下用PBS-T洗涤，并用70％和100％乙醇脱水
。除去杂交室 ，并用慢速金抗固定器（S36937
，Thermo）安装滑梯。在此阶段的CY5标记RCP的成像是确认RCP（“锚”）的质量控制步骤 在解码过程中产生并有助于识别斑点 。使用Perkin Elmer Opera Phenix加上1-μmZ-STEP尺寸的高素质筛选系统进行成像，使用63×（Na 1.15 ，0.097μm Pixel-1）水浸水物目标进行成像
。For channels: DAPI (excitation of 375 nm and emission of 435–480 nm), Atto 425 (excitation 425 nm and emission 463–501 nm), Alexa Fluor 488 (excitation 488 nm and emission 500–550 nm), Cy3 (excitation 561 nm and emission 570–630 nm) and Cy5 (excitation 640 nm和排放650–760 nm）。成像后
，每张载玻片在PBS中垂直垂直覆盖（轻轻地，搅动最小 ，使盖玻片“掉下来 ”以防止对组织损坏） 。该切片用70％和100％乙醇脱水
，并将新的杂交室固定到幻灯片上。使用100％甲酰胺（AM9342，Thermo）剥离了上一个周期 ，该甲酰胺每分钟新鲜施加5分钟
，然后用PBS-T洗涤
。使用两轮杂交进行条形码标记，首先是适配器探针池（随后的周期中 ，AP混合AP1-AP6），然后进行测序池（SP Mix
，使用ATTO 425自定义），每个在37°C下在中间和之后的PBS-T洗涤37°C持续1 h。该部分如前所述脱水 ，拆除腔室
，并将幻灯片安装和成像 。再重复五次，以生成七个周期（锚和六个条形码位）的完整数据集
。 　　在室温下在EDTA 0.5 m的EDTA中 ，通过在1周内孵育标本在室温下通过孵育。在孵育期间的一半将溶液更新 。在1天内将样品在PBS 1X中两次洗涤
。将样品在室温下以H2O中升高的甲醇（20％
，40％，60％和80％）脱水1小时。然后 ，将样品放在白光（11 W和3,000 K°）和滚动搅拌下（004011000
，IKA），在100％甲醇中用6％过氧化氢溶液（004011000 ，IKA）放置过夜 。将样品在室温下以下降浓度的甲醇（80％
，60％，40％和20％）的浓度将样品复合1小时 ，在1周内用0.5％Triton X-100（PBSGT）溶液中洗涤两次，并在0.2％PBS-明智素中封闭
。将样品转移到含有原代抗体（Osterix
，1/500; AB209484，ABCAM和胶原蛋白2 ，1/500; AB185430
，ABCAM）的溶液中，并在PBSGT中稀释 ，并在37°C下以搅拌20 rpm孵育14天。接下来是在室温下在PBSGT中洗1小时的六个洗涤 。接下来
，将二抗在PBSGT中稀释并通过0.22-μm滤波器
。将样品在二级抗体溶液中在37°C下孵育7天，并在室温下在PBSGT中1小时洗涤六次。 　　IDISCO+协议用于清除样品64 。将样品放入TPP（Techno塑料产品）管中 ，在旋转搅拌下（SB3
，Stuart）下甲醇（20％，40％ ，60％
，80％和100％（2x））脱水1小时。甲醇量等于样品体积的五倍
。接下来将样品孵育在67％DCM和33％MEOH的溶液中，然后在旋转器上的室温下100％DCM孵育30分钟 ，然后放入100％DBE。用IMSpector Pro 7.5.3录取软件（Miltenyi Biotec）控制的配备SCMOS摄像头5.5MP（2,560×2,160像素）配备了SCMOS摄像头5.5MP（2,560×2,160像素）的燃烧光显微镜（Miltenyi Biotec）成像 。厚度为4 µm的光板是由四个不同波长（488 nm，561 nm
，639 nm和785 nm）产生的
。使用了1×或4×目标，具有不同放大镜头×0.6 ，×1和×1.66的目标。样品由Miltenyi的样品支架支撑
，并将其放置在装满DBE的水箱中，并根据样品的大小从一个或两侧的激光灯片照亮 。在获取期间使用Lightspeed模式在合理的时间和合适的分辨率中获取这些图像
。3D图像瓷砖的镶嵌物在瓷砖之间的重叠组装为10％。图像以16位TIFF格式获取 。最初使用Macs IQ View软件处理图像 ，该软件执行了瓷砖的自动对齐
。使用ImarisfileConverter将堆栈图像转换为Imaris文件（.IM）。为了隔离Imaris中的特定结构
，我们使用了带有手动选择的表面工具，然后使用表面掩盖了图像 。通过使用Imaris中的功能快照和动画拍摄图像和视频
。Adobe Photoshop（v25.2）用于为缝合区域着色。 　　根据制造商的指示 ，使用Leica Bond RX使用Leica Bond RX处理与RNASCOPE多重荧光试剂试剂盒V2测定法（高级细胞诊断和生物技术）自动化染色 。可以在补充表7中找到探针。在染色之前
，将新鲜的冷冻切片在4％的PBS中在45分钟的4°C中固定45分钟，然后通过一系列的50％ ，70％
，100％，100％和100％乙醇脱水 ，持续5分钟
。在手动预处理后，自动处理包括在探针杂交之前对蛋白酶III的消化15分钟。使用OPAL 520
，OPAL 570和OPAL 650（Akoya Biosciences），TSA-Biotin（TSA Plus Biotin kit ，Perkin Elmer）和链霉亲和链霉素偶联的Atto 425（Sigma-Aldrich）
，使用OPAL 520，OPAL 570和OPAL 650（Akoya Biosciences） ，TSA-BIOTIN（SIGMA-ALDRICH）
，用于开发Rnascope探针通道 。与上述ISS一样，对染色部分进行了成像
。 　　如前所述 ，我们应用了新鲜供体组织的全细胞解离5。在提取细胞之前
，解剖样品组织（大约9个PCW肩关节）以获得骨样品，并将软组织显微解剖 。将所得的细胞悬浮液用DAPI（Invitrogen）染色 ，用于实时可持续性
，FGFR3抗体（1:50; MA5-38521，Thermo Fisher Scientific）和TACR3抗体（1:50; BS-0166R ，Thermo Fisher Scientific）和二级抗体
。使用Bigfoot Spectral Sorter（Thermo Fisher Scientific）及其专有软件，将FACS门控DAPI阳性单元细胞进行DAPI染色。然后进行FGFR3和TACR3的顺序门控以鉴定双阳性细胞 。使用人外周血单核细胞进行FGFR3和TACR3的阳性对照
，并将未染色的细胞用作阴性对照
。 　　我们使用了Li等人65中概述的ISS解码管道。该管道由五个不同的步骤组成 。首先，我们使用Perkin Elmer提供的Acapella脚本进行了图像缝制 ，该脚本为每个周期生成了所有通道的二维最大强度投影
。接下来
，我们使用MicroAligner66（V1.0.0）使用默认参数根据DAPI信号进行注册。对于细胞分割，我们利用了一种利用CellPose67（v3.0）作为分割方法的可扩展算法 。预期的细胞大小设置为直径为70像素 ，并进一步扩展了10像素以模仿细胞质。To decode the RNA molecules, we used the PoSTcode algorithm68 (v1.0) with the following parameters: rna_spot_size = 5, prob_threshold = 0.6, trackpy_percentile = 90 and trackpy_separation = 2. Furthermore, we assigned the decoded RNA molecules to segmented cells using STRtree (v2.0.6) and subsequently generatedAnndata对象遵循Virshup等人描述的方法
。69。最后
，只保留了超过四个RNA分子的细胞进行下游分析 。 　　根据制造商的方案进行了新鲜冷冻（10倍基因组学）的覆覆质细胞助手空间基因表达
。基于与液滴数据相关的骨形成的微环境（例如冠状缝合线）选择了感兴趣的区域，并相应地与细胞兼机器垫圈对齐。在进行visium cytassist方案以进行后续比对之前 ，使用Hamamatsu S60滑动扫描仪以×40放大倍率捕获图像 。库与间隔剂（10倍基因组学）映射
。 　　如先前的工作70所述
，通过冷冻样品从新鲜的冷冻样品中分离出单核。In brief, 10-μm sections were homogenized in homogenization buffer (250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM Tris-HCl, 1 mM dithiothreitol, 1× protease inhibitor, 0.4 U μl−1 RNaseIn, 0.2 U μl−1 SUPERaseIn and 0.1% Triton X-100 in nuclease-free water) using a glass Dounce组织研磨机套件（默克） 。将样品与10–20杆的悬一面的杵“ A
”分离，然后将10个紧身的杵“ B”击中
。将所得的混合物通过50μm的细胞过滤器 ，然后离心（500G 5分钟）
，然后将沉淀重悬于300μl的储存缓冲液中（1×PBS，4％BSA和0.2UμL-1 Protecor RNasein） ，并通过20μm细胞菌株。在用锥虫蓝溶液染色后，可视化细胞核并评估在显微镜下的生存力
，并根据制造商的方案进一步处理10倍基因组铬单细胞Multiome ATAC +基因表达 。每个反应的靶向核恢复16,000滴
。对于某些核样品，将来自不同样品供体的样品的混合物集合在一个反应​​中 ，然后通过基因型消除。使用Bioanalyzer高灵敏度DNA分析（Agilent）进行cDNA和最终文库的质量控制 。使用Novaseq 6000（Illumina）对库进行测序
，其测序深度为每滴20,000个读取对。 　　使用CellRanger-ARC软件（v2.0.0），将测序数据与人参考基因组（GRCH38-2020-A-2.0.0）对齐
。来自10倍多层泳道的称为条形码 ，使用SOOLCELL（v2.0）71
，通过BAM输出将来自多个样品供体的合并基因型的基因型取消。随后，通过基因型聚类 ，将汤鼠输出聚集
，以分配每个条形码 。使用默认的输入设置将visium数据映射到间距（v1.1.0），并根据捕获标记区域使用10X Genomics LoupeBrowser（v7.0）将低分辨率的细胞助手图像与处理幻灯片的HI分辨率显微镜图像对齐
。对于基因表达数据 ，SOODX（V1.6.0）72被应用以去除背景环境RNA。对于称为矩阵的Cellranger-Arc
，其中包含超过16,000滴（超过目标液滴回收的预期数），我们将估计的全球RHO值增加了0.1 ，以说明其他环境RNA的潜力 。将液滴用于200多个基因，小于5％的线粒体和核糖体读数
。Doublet删除如下所述。对于单细胞ATAC-SEQ
，我们应用ARCHR73（V1.0.2）来处理CellRanger-ARC的输出 。考虑到独特的核碎片，信噪比和片段尺寸分布的数量 ，进行了初始的人均质量控制
。此外
，带有转录开始站点富集的液滴< 7 and number of fragments < 1,000 were removed. Doublets were discarded using the default settings. Initial clustering was performed at a resolution of 0.2 with the top 40 dimensions from iterative latent semantic indexing. Then, pseudo-bulk replicates were made for each broad cluster per region from the initial clustering results. Peak calling (501-bp fixed-width peaks) was performed based on pseudo-bulk coverages by MACS2 (v2.2.7.1). Then, a cell-by-peak count matrix was obtained and exported. We applied muon74 (v0.1.2) for normalization, latent semantic indexing dimension reduction and clustering analysis using BBKNN75 (v1.5.1) to correct for batch effects from anatomical regions and sample donors to obtain an ATAC embedding. Gene scores based on chromatin accessibility around gene bodies were calculated. We then applied MultiVI76 (via scVI v0.6.8) to construct a joint embedding for snRNA-seq and single-cell ATAC-seq. We also applied EmptyDropMultiome77 (v1.0.0) to repeat droplet calling to identify nucleus-containing droplets in our Multiome data to reduce the ambient RNA noise (‘soup’). By generalizing EmptyDrops to the multi-omic setting, we used the smallest droplets to create an RNA and an ATAC ambient RNA ‘soup’ profile, and then tested each droplet for statistical deviations from each of these two profiles, retaining only droplets that were statistically significantly different from the soup profile. 　　All potential doublets detected in both RNA and ATAC modalities were removed from our data. For RNA data, Scrublet (v0.2.3)78 was applied to estimate doublet probability, and a score of more than 0.3 was used as a cut-off value. To apply a stringent doublet threshold, we conducted an adapted Scrublet workflow as previously described79. In brief, per-droplet Scrublet scores were first determined for CellRanger-ARC count matrices from each 10X Multiome (gene expression) lane independently. The droplets were then overclustered through the standard scanpy workflow using default parameters up to Leiden clustering. Each individual cluster was further clustered. A per-cluster median of Scrublet scores was computed. A normal distribution of doublet score, centred at the score median with a standard deviation estimated from the median absolute deviation, was used to compute P values for each of the clusters. After false discovery rate adjustment using Benjamini–Hochberg correction, a P >0.65被认为是优质细胞的截止值，因为双重距离是显着的异常值。对于ATAC数据 ，我们首先从ARCHR应用了Doublet检测方法来删除假定的ATAC Doublets 。此外
，同型和异型双重运动表的特征是在单个SNATAC-SEQ库上运行护身符（v1.1.0），以及带有Q的液滴< 0.01 were removed. 　　We adopted a hierarchical clustering approach by first conducting Leiden clustering on the global integrated scVI (v0.9.1; hidden layers = 256, latent variables = 52, dispersion = ‘gene-batch’) RNA embeddings to obtain broad clusters. To validate these, we used Celltypist to train a model on cell states in the embryonic limb bud5,80,81, and transferred labels onto our embedding for inspection. We utilized this information in addition to canonical marker genes to annotate broad clusters and subset sublineages. For sublineages (chondrocytes, fibroblasts, osteogenesis-related clusters, Schwann cells, immune cells and endothelial cells), we further embedded each subset using scVI (hidden layers = 256, latent variables = 52 and dispersion = ‘gene-batch’) and conducted Leiden clustering (resolution = 0.6), followed by differentially expressed gene (DEG) analyses (method = ‘wilcoxon’) to obtain cluster markers. We additionally utilized the inferred spatial location of cell states (described below) to inform annotations. 　　We applied the Python implementation of the MILO package (v0.1.1) for differential abundance testing (http://github.com/emdann/milopy)82. We used the scVI latent representation to create a k-nearest neighbour graph of droplets in the relevant compartment and subsequently applied milopy to allocate droplets to overlapping neighbourhoods, with these droplets originating from multiple samples (brc_code). Each neighbourhood was then annotated as a cluster based on majority voting. We binarized values for anterior–posterior positions and calvarium-appendicular covariates to allow testing across these variables. We then determined log fold-change values for differential abundance and false discovery rate values based on the Bejanmini–Hochberg correction. 　　We performed Cell2location (v0.1.4) for deconvolution of Visium CytAssist voxels using our annotated Multiome data as inputs. Sample donor was used as the batch variable, and each library was considered a covariate in the regression model. For spatial mapping, we estimated 30 cells per voxel based on histological data, and set a hyperparameter detection alpha value of 20 for per-voxel normalization. 　　ISS-Patcher is a package for approximating features not experimentally captured in low-dimensional data based on related high-dimensional data. It was developed as an approach to approximate expression signatures for genes missing in ISS data using matched snRNA-seq data as a reference in this study. First, a shared feature space between both datasets was identified by subsetting the 155–158 genes present in the ISS pool, followed by separate normalization to median total cell counts, log-transformation and z-scoring for both modalities. Then, the 15 nearest neighbours in scRNA-seq space were identified for each ISS cell with the Annoy Python package, and the genes absent from ISS were imputed as the average raw counts of the scRNA-seq neighbours. 　　Our Visium cranium sample was annotated with TissueTag8 using a semi-automatic mode to generate a one-dimensional maturation axis. Regions of the developing bone were first manually annotated based on haematoxylin and eosin features. Tissue regions that did not include bone-forming niches were excluded from annotation. The annotation categories that were stored included multiple regions of the coronal suture (level 0 to level 2 annotation), stemming from the central-most portion, an osteogenic front (level 3 annotation) with histological features of osteoprogenitors and osteogenic zones (level 4 to level 7 annotation) from the emergence of histological osteoblasts. All annotations were saved as TissueTag output format, which logs the annotation resolution, the pixels per micrometre and the pixel value interpretation of annotation names (for example, 0 = ‘suture’) and colours (for example, ‘osteogenic front’: ‘red’). To robustly and efficiently migrate TissueTag annotations to the Visium objects, we first transferred TissueTag annotations from pixel space to a high-resolution hexagonal grid space (15-µm spot diameter and 15-µm point-to-point centre distance with no gap between spots) based on the median pixel value of the centre of the spot (radius/4) in the annotated image. Next, to generate continuous annotations for Visium data, for each spot in the hexagonal high-resolution grid, we measured the mean Euclidean distance to the ten nearest points from each annotated structure in the level 0 annotation and the distance from the closest point for structures in level 1 annotation. All annotations were mapped to the Visium spots by proximity of the spot annotation grid to the nearest corresponding spot in the Visium array. 　　The SCENIC+83 (v1.0.0) pipeline was used to predict transcription factors and putative target genes as well as regulatory genomic regions with binding sites. The fragment matrix of peaks called with MACS2 and processed within ArchR73 together with the corresponding RNA count matrix were used as inputs. Meta-cells were created by clustering droplets into groups of around 10–15 droplets based on their RNA profiles and subsequent aggregation of counts and fragments. The pipeline was applied to subsets of the dataset corresponding to individual lineages: first, CisTopic (pycistopic v1.0.2) was applied to identify region topics and differentially accessible regions from the fragment counts as candidate regions for transcription factor binding. CisTarget (pycistarget v1.0.2) was then run to scan the regions for transcription factor-binding sites, and GRNBoost2 (arboreto v0.1.6)84 was used to link transcription factors and regions to target genes based on co-expression or accessibility. Enriched transcription factor motifs in the regions linked to target genes were used to construct transcription factor–region and transcription factor–gene regulons. Finally, regulon activity scores were computed with AUCell based on target gene expression and target region accessibility, and regulon specificity scores derived from them. Networks of transcription factors, regions and target genes (enhancer-driven GRNs) were constructed by linking individual regulons. Transcription factor–enhancer–gene links for all subsets (osteogenesis, chondrogenesis, fibrogenesis, early joint progenitors, immune and Schwann) can be found in Supplementary Table 9. 　　For pseudotime trajectory construction in the osteogenic subcompartment, non-cycling droplets were subsetted, and X_scVI was used as projections for palantir to obtain multiscale diffusion space. A neighbourhood graph was generated on the diffusion space using sc.pp.neighbors, and the first two principal components were used as initial positions to create ForceAtlas2 embeddings using sc.tl.draw_graph. scFates85 (v1.0.3) was used to predict a principle graph that captures the differentiation path. The force-directed embeddings and principle graph were exported into R, and monocle3 (v1.0.0)86,87 was used to compute differentially expressed genes along pseudotime using a graph-based test (morans’ I)87,88, which allows identification of genes upregulated at any point in pseudotime. The results were visualized with heatmaps using the complexHeatmap (v2.6.2)89 and seriation (v1.3.0)90 packages, after smoothing gene expression with smoothing splines in R (smooth.spline; d.f. = 12). Velocity analysis91 was performed using scvelo92 (v0.2.3). Spliced and unspliced read counts were computed with velocyto (v0.17.17) from the unprocessed data, before using scvelo.pp.moments, scvelo.tl.velocity and scvelo.tl.velocity_graph to compute velocities for the preprocessed droplets. cytoTRACE93 was used (through the CellRank94 (v2.0.2) implementation) to obtain another prediction of directionality, independent of RNA velocity (based on the assumption that the number of expressed genes decreases throughout differentiation). 　　To approximate the timing of cavitation onset, we computed a cavitation enrichment score using sc.tl.score_genes() in scanpy on a specific gene set within the mesenchymal and muscle compartments of the hip, shoulder and knee joints comprising CD44, HAS2, ABCC5, HMMR, MSN and UDPGD, derived from literature and Gene Ontology terms, which encompass hyaluronan biosynthetic processes and hyaluronan synthase activity. We excluded genes with low expression levels in our data, such as HAS3. For pathway analysis, we retrieved gene sets corresponding to all 18,640 Gene Ontology terms, and computed the correlation between their enrichment scores and cavitation enrichment scores. 　　CellOracle95 (v0.12.0) was used with the osteogenesis trajectory created with scFates85, and the regulons predicted with SCENIC+83 for the same cells were imported into CellOracle as a base GRN. Cells were aggregated into meta-cells of 10–15 cells, and linear models explaining transcription factor from target gene expression were fitted with CellOracle per cell cluster. Regulon-based transcription factor perturbation vectors were inferred using the cell cluster-specific models to predict effects of transcription factor overexpression and knockout. Diffusion pseudotime96 was then computed for intramembranous and endochondral ossification lineages separately by selecting corresponding starting points. The pseudotime gradients were used to derive pseudotime-based differentiation vectors, and the pseudotime-perturbation vector cross-product was computed to obtain perturbation scores. These perturbation scores indicate whether the in silico perturbation of a transcription factor is consistent with or opposes differentiation along a lineage (osteogenesis). The simulations were carried out systematically, overexpressing and knocking out all transcription factors in the GRN. For each transcription factor and condition, the perturbation scores were then averaged per cell cluster and summarized in a table to screen for transcription factors promoting or inhibiting osteogenesis. 　　fGWAS analysis97 was applied to identify disease-relevant cell clusters as described in detail55 (https://github.com/natsuhiko/PHM). The model makes use of full summary statistics from GWAS, linking single-nucleotide polymorphisms (SNPs) to genes, and captures a general trend between gene expression and disease association of linked loci for each cell cluster. At the same time, the model also corrects for linkage disequilibrium and other relevant factors. We used full GWAS summary statistics obtained from the EBI GWAS Catalog, open targets, and knee and hip osteoarthritis as well as total knee and hip replacement from ref. 62 (https://msk.hugeamp.org/downloads.html; Supplementary Table 8). 　　We used a network propagation98 approach to integrate GWAS summary statistics and cell cluster marker gene-based scores for prioritizing disease-relevant and cell cluster-specific subunits of our transcription factor network. First, scores per SNP were computed from downloaded summary statistics and weighted by linkage disequilibrium. Then, the scores were mapped to a GRN, here an enhancer-driven GRN computed with SCENIC+ for the corresponding lineage. As the used networks contain transcription factors and target genes, and also regions with transcription factor-binding sites as nodes, SNP scores were mapped to both genes and regions, representing distal regulatory elements. The scores were then propagated across the network using a random walk with restart (or personalized page rank) process. This integrates the contribution of individual SNPs, with signals converging around relevant network nodes. The procedure was repeated with 1,000 randomly permuted scores to compute permutation-test results and z-scores. Next, differential expression-based marker gene scores for each cell cluster were propagated in the same way, resulting in cell cluster specificity scores for each network node. The SNP and expression-based scores were then combined per cell cluster (as in ref. 99) by using the minimum for each node. The final scores were thresholded, and the resulting connected components were obtained as enriched subnetworks. The method has been compiled into a tool that we called SNP2Cell, which is available at https://github.com/Teichlab/snp2cell. 　　Ligand–receptor interactions were inferred using ‘cpdb_analysis_method.call’ in CellPhoneDB (v4.0.0). We included genes expressed in more than 10% of cells within each cluster. Inferred interactions with a P >删除0.001
。我们使用Nichenet（v1.1.1）来识别内侧软骨和膜内壁ni之间的不同相互作用。我们首先使用在Seurat中实现的Wilcoxon测试来计算两个壁ni的成骨簇和尖端细胞的DEG ，并使用每个群集的最小对数折叠变化来总结差异表达的配体和受体 。前1,000摄氏度用于计算配体活动。我们使用默认设置优先考虑配体 - 对象链接
。使用热图可视化前十个配体及其顶级得分受体。 　　从Chembl数据库（https://www.ebi.ac.uk/chembl/）收集了人类的药物和靶基因信息（HOMO SAPIENS） 。对于靶向致病性药物
，我们搜索了针对编码其报告靶标的基因的临床认可的分子，并从Chembl数据库中策划了65种临床批准的药物列表 ，该数据库提出了致畸性的警告（补充表6）
。如前所述58计算药物评分。随后
，我们根据广泛的临床实用性引入了每种药物的药物类别 。Drug2cell Python软件包可在GitHub（https://github.com/teichlab/drug2cell）上找到
。 　　首先，根据对多组数据执行的工作流程 ，使用Starsolo将Zhang等人的FASTQ文件使用STARSOLO转换为常见的基因组参考（GRCH38-2020-A-2.0.0）。应用细胞板来删除表示为模拟环境RNA的背景计数 。We intersected this matrix with barcodes from the post-quality control counts matrix from Zhang et al., and scVI was then used to integrate this with our snRNA-seq data, accounting for categorical covariates of sample donor and droplet technology (cell or nuclei), as well as continuous variables of total counts, the percentage of ribosomal and mitochondrial counts, and cell cycle scores使用Scanpy软件包计算的“ S_SCORE ”和“ G2M_SCORE
”
。然后使用从此获得的潜在变量来确定邻域，然后确定UMAP的维度降低。Zhang等人的集群标签 。然后用作CellHint.harmonize（）比对函数中的cellhint协调的标签
。CellHint.Treeplot（）用于检查和半自动对两个数据集的标签对齐。标记基因的基因表达谱用于验证两个数据集中簇的比对 。 　　有关研究设计的更多信息可在与本文有关的自然投资组合报告摘要中获得。