scDRS (single-cell disease-relevance score)

If your gene sets (.gs file) are in a different species than your single-cell data (.h5ad file), you must convert the gene symbols. scDRS provides a mapping utility (typically using mouse_human_homologs.txt) to translate between mmusculus and hsapiens.

Supported conversions:

• Human (hsapiens) $\rightarrow$ Mouse (mmusculus)
• Mouse (mmusculus) $\rightarrow$ Human (hsapiens)

# Example logic for mapping genes between species
# dic_map is created from mouse_human_homologs.txt
for trait in df_gs.index:
    gs_gene_list = df_gs.loc[trait, 'GENESET'].split(',')
    # Intersect GS genes with the mapping dictionary keys
    h5ad_gene_list = [dic_map[x] for x in set(gs_gene_list) & set(dic_map.keys())]
    df_gs.loc[trait, 'GENESET'] = ','.join(h5ad_gene_list)

scDRS scales linearly with both the number of cells and the number of control gene sets for both computation and memory usage.

Example Scale: A single-cell dataset with one million cells requires approximately 3 hours of computation and 60GB of memory.

To ensure statistically valid and powerful scDRS results, follow these recommendations:

Gene Set Size

• Ideal size: Moderate size, specifically >50 genes and <20% of all genes.
• Minimum threshold: Performance is generally reasonable as long as the gene set size is >=10 genes.
• Other uses: You can use non-GWAS gene sets to identify cells with excess expression of genes within that set.

GWAS and scRNA-seq Data

• GWAS: Use data with a heritability z-score > 5 or a sample size > 100K to ensure reasonable discovery rates.
• scRNA-seq: Use data with a diverse set of cells potentially relevant to the disease. Note that a smaller number of cells should not negatively affect scDRS power.

To install version 1.0.0, clone the repository, checkout the tag, and install in editable mode.

git clone https://github.com/martinjzhang/scDRS.git
cd scDRS
git checkout -b rv1 v1.0.0
pip install -e .

You can compute scDRS scores using provided bash scripts. For the TMS FACS dataset with 74 diseases/traits, ensure you set the DATA_PATH environment variable to your local folder containing the TMS FACS data and .gs files before running the script.

export DATA_PATH=/path/to/scDRS_data_release_092121
./reproduce_compute_score.tms_facs_with_cov.magma_10kb_1000.sh

If you require the versions used in the initial submission (Binary GWAS gene sets instead of Weighted), use these links:

• Binary GWAS Gene Sets: Download .gs files for 74 diseases and complex traits from Figshare.
• scDRS Results (Binary): Download .score.gz and .full_score.gz files for TMS FACS + 74 diseases/traits from Figshare.

After installing version 0.1, run the core scDRS tests using pytest to ensure the installation is correct.

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

The job.simulation directory contains tools to generate data and compare scoring methods.

Workflow

1. Generate Data: Use generate_null_simulation_data.ipynb for null data or generate_causal_simulation_data.ipynb for causal configurations.
2. Compute Scores: Scripts are available for various methods:
   • compute_simu_score.sh (scDRS)
   • compute_simu_score_scanpy.sh (Seurat)
   • compute_simu_score_vision.sh (Vision)
   • call_R_vam.sh (VAM)
   • compute_perturb_simu_score.sh (Causal simulations)
3. Visualize: Use make_figure.null_simulation.ipynb or make_figure.causal_simulation.ipynb to generate comparison figures.

The following Jupyter notebooks are available for specific analysis tasks related to the CA1 Pyramidal heterogeneity experiment:

• Replicate Figure 5AB: Use main_figure.ipynb.
• TMS FACS dataset analysis: Use tms.ipynb to analyze neurons in the TMS FACS dataset.
• Zeisel et al. 2015 dataset analysis: Use zeisel.ipynb to analyze the Zeisel et al. 2015 dataset.
• Spatial coordinate verification: Use spatial_verify.ipynb to verify inferred spatial coordinates.

Calculate LD scores for specific datasets (ts_facs, tms_facs, tms_droplet) using Slurm sbatch arrays.

To run the calculation for all cell types in a dataset:

for dset in ts_facs tms_facs tms_droplet ; do
    n_ct=$(cat out/${dset}/ct_list.txt | wc -l)
    sbatch --array=1-${n_ct} calc_ldscore.sh ${dset}
done

To check for unfinished jobs and only submit missing LD score files (out/${dset}/ldscore/${ct}.22.l2.ldscore.gz):

for dset in ts_facs tms_facs tms_droplet; do
    n_ct=$(cat out/${dset}/ct_list.txt | wc -l)
    for i in $(seq 1 $n_ct); do
        ct=$(sed "${i}q;d" out/${dset}/ct_list.txt)
        if [ ! -f out/${dset}/ldscore/${ct}.22.l2.ldscore.gz ]; then
            sbatch --array=$i calc_ldscore.sh ${dset}
        fi
    done
done