sc-DRS: 连接疾病GWAS与单细胞数据的分析方式

git clone https://github.com/martinjzhang/scDRS.git
cd scDRS
pip install -e .

python -m pytest tests/test_CLI.py -p no:warnings

wget https://figshare.com/ndownloader/files/34300925 -O data.zip
unzip data.zip && \
mkdir -p data/ && \
mv single_cell_data/test/* data/ && \
rm data.zip && rm -r single_cell_data/test/

import scdrs
import scanpy as sc
sc.set_figure_params(dpi=125)
from anndata import AnnData
from scipy import stats
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import os
import warnings

warnings.filterwarnings("ignore")

adata =  sc.read_h5ad("scDRS/data/expr.h5ad")

df_gs = pd.read_csv("single_cell_data/test/geneset.gs", sep="\t", index_col=0)`#改成自己下载的geneset.gs文件所在路径

df_gs = df_gs.loc[
    [
        "PASS_Schizophrenia_Pardinas2018",
        "spatial_dorsal",
        "UKB_460K.body_HEIGHTz",
    ],
    :,
].rename(
    {
        "PASS_Schizophrenia_Pardinas2018": "SCZ",
        "spatial_dorsal": "Dorsal",
        "UKB_460K.body_HEIGHTz": "Height",
    }
)
display(df_gs)

df_gs.to_csv("data/processed_geneset.gs", sep="\t")

scdrs compute-score \
    --h5ad-file /Users/han/scDRS/data/expr.h5ad \
    --h5ad-species mouse \
    --gs-file /Users/han/single_cell_data/test/processed_geneset_capitalized.gs \
    --gs-species mouse \
    --cov-file /Users/han/scDRS/data/cov.tsv \
    --flag-filter-data True \
    --flag-raw-count True \
    --flag-return-ctrl-raw-score False \
    --flag-return-ctrl-norm-score True \
    --out-folder /Users/han/single_cell_data/test/

    trait: pd.read_csv(f"single_cell_data/test/{trait}.full_score.gz", sep="\t", index_col=0)
    for trait in df_gs.index
}

for trait in dict_score:
    adata.obs[trait] = dict_score[trait]["norm_score"]

sc.set_figure_params(figsize=[3, 3], dpi=300)
sc.pl.umap(
    adata,
    color="level1class",
    ncols=1,
    color_map="RdBu_r",
    vmin=-5,
    vmax=5,
)

sc.pl.umap(
    adata,
    color=dict_score.keys(),
    color_map="RdBu_r",
    ncols=3,
    vmin=-5,
    vmax=5,
    s=20,
)

for trait in ["SCZ", "Height"]:
    !scdrs perform-downstream \
        --h5ad-file scDRS/data/expr.h5ad \
        --score-file single_cell_data/test/{trait}.full_score.gz \
        --out-folder single_cell_data/test/ \
        --group-analysis level1class \
        --flag-filter-data True \
        --flag-raw-count True

!cat single_cell_data/test/SCZ.scdrs_group.level1class | column -t -s $'\t'

dict_df_stats = {
    trait: pd.read_csv(f"/Users/han/single_cell_data/test/{trait}.scdrs_group.level1class", sep="\t", index_col=0)
    for trait in ["SCZ", "Height"]
}

dict_celltype_display_name = {
    "pyramidal_CA1": "Pyramidal CA1",
    "oligodendrocytes": "Oligodendrocyte",
    "pyramidal_SS": "Pyramidal SS",
    "interneurons": "Interneuron",
    "endothelial-mural": "Endothelial",
    "astrocytes_ependymal": "Astrocyte",
    "microglia": "Microglia",
}

fig, ax = scdrs.util.plot_group_stats(
    dict_df_stats={
        trait: df_stats.rename(index=dict_celltype_display_name)
        for trait, df_stats in dict_df_stats.items()
    },
    plot_kws={
        "vmax": 0.2,
        "cb_fraction":0.12
    }
)


ax.grid(False) 

plt.show()

#提取 CA1 pyramidal并进行重新聚类
adata_ca1 = adata[adata.obs["level2class"].isin(["CA1Pyr1", "CA1Pyr2"])].copy()
sc.pp.filter_cells(adata_ca1, min_genes=0)
sc.pp.filter_genes(adata_ca1, min_cells=1)
sc.pp.normalize_total(adata_ca1, target_sum=1e4)
sc.pp.log1p(adata_ca1)

sc.pp.highly_variable_genes(adata_ca1, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata_ca1 = adata_ca1[:, adata_ca1.var.highly_variable]
sc.pp.scale(adata_ca1, max_value=10)
sc.tl.pca(adata_ca1, svd_solver="arpack")

sc.pp.neighbors(adata_ca1, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata_ca1, n_components=2)

# 匹配 scDRS 分数
for trait in dict_score:
    adata_ca1.obs[trait] = dict_score[trait]["norm_score"]
sc.pl.umap(
    adata_ca1,
    color=dict_score.keys(),
    color_map="RdBu_r",
    vmin=-5,
    vmax=5,
    ncols=3,
    s=20,
    save="Subgroup-Disease.png"
)

df_plot = adata_ca1.obs[["Dorsal", "SCZ", "Height"]].copy()
df_plot["Dorsal quintile"] = pd.qcut(df_plot["Dorsal"], 5, labels=np.arange(5))

fig, ax = plt.subplots(figsize=(3.5, 3.5))
for trait in ["SCZ", "Height"]:
    sns.lineplot(
        data=df_plot,
        x="Dorsal quintile",
        y=trait,
        label=trait,
        err_style="bars",
        marker="o",
        ax=ax,
    )
ax.set_xticks(np.arange(5))
ax.set_xlabel("Dorsal quintile")
ax.set_ylabel("Mean scDRS disease score")
fig.show()

spatial_col = "Dorsal"
for trait in ["SCZ", "Height"]:
    df_score = dict_score[trait].reindex(adata_ca1.obs.index)
    ctrl_cols = [col for col in df_score.columns if col.startswith("ctrl_norm_score")]
    n_ctrl = len(ctrl_cols)

    # Pearson's r between trait score and spatial score
    data_r = stats.pearsonr(df_score["norm_score"], adata_ca1.obs[spatial_col])[0]

    # Regression: control score ~ spatial score
    ctrl_r = np.zeros(len(ctrl_cols))
    for ctrl_i, ctrl_col in enumerate(ctrl_cols):
        ctrl_r[ctrl_i] = stats.pearsonr(df_score[ctrl_col], adata_ca1.obs[spatial_col])[
            0
        ]
    pval = (np.sum(data_r <= ctrl_r) + 1) / (n_ctrl + 1)

    print(f"{trait} v.s. {spatial_col}, Pearson's r={data_r:.2g} (p={pval:.2g})")

•	基于GWAS导出的基因集，为每个单细胞计算疾病相关性评分（Disease Relevance Score, DRS）；
•	鉴定疾病富集的细胞亚群，揭示特定细胞在疾病发生发展中的作用；
•	支持不同疾病、多种组织、复杂细胞群体之间的系统性比较；
•	辅助发现潜在治疗靶点或生物标志物，加速疾病机制研究。