单细胞转录组实战06: pySCENIC转录因子分析（docker）

cd output/05.SCENIC
micromamba activate SC
nohup zsh ~/Project/SC10X/src/run_SCENIC.sh &> run_SCENIC.sh.log &

docker运行脚本

• docker脚本

#>>>scenic_docker.sh>>>
n_jobs=20
input_dir=/home/victor/DataHub/SCENIC
output_dir=/home/victor/Project/SC10X/output/05.SCENIC
# STEP 1/3: Gene regulatory network inference, and generation of co-expression modules
docker run --rm \
    -v ${input_dir}:/input_data \
    -v ${output_dir}:/output_data \
    aertslab/pyscenic:0.12.1 pyscenic grn \
        --num_workers ${n_jobs} \
        --method grnboost2 \
        --output /output_data/grn.csv \
        --sparse \
        /input_data/X.loom \
        /input_data/hs_hgnc_tfs.txt
# STEP 2/3: Regulon prediction (cisTarget)
docker run --rm \
    -v ${input_dir}:/input_data \
    -v ${output_dir}:/output_data \
    aertslab/pyscenic:0.12.1 pyscenic ctx \
        /output_data/grn.csv \
        /input_data/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
        --mask_dropouts \
        --sparse \
        --annotations_fname /input_data/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
        --expression_mtx_fname /output_data/X.loom \
        --mode "custom_multiprocessing" \
        --output /output_data/ctx.csv \
        --num_workers ${n_jobs}
#<<<scenic_docker.sh<<<

cd output/05.SCENIC
nohup zsh ~/Project/SC10X/src/scenic_docker.sh &> scenic_docker.sh.log &

from pathlib import Path
import operator
import cytoolz
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
from pyscenic.utils import load_motifs
from pyscenic.aucell import aucell
from pyscenic.binarization import binarize
from pyscenic.plotting import plot_binarization, plot_rss
from pyscenic.transform import df2regulons
import bioquest as bq #https://jihulab.com/BioQuest/bioquest

OUTPUT_DIR='output/05.SCENIC'
Path(OUTPUT_DIR).mkdir(parents=True,exist_ok=True)

adata = sc.read_h5ad('output/03.inferCNV/adata.h5')
adata_raw = adata.raw.to_adata()

def filter_regulons(motifs, db_names=("hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings",)):
    从ctx.csv中筛选重要的regulons，之后再运行AUCell
    motifs.columns = motifs.columns.droplevel(0)
    def contains(*elems):
        def f(context):
            return any(elem in context for elem in elems)
        return f
    # For the creation of regulons we only keep the 10-species databases and the activating modules. We also remove the
    # enriched motifs for the modules that were created using the method 'weight>50.0%' (because these modules are not part
    # of the default settings of modules_from_adjacencies anymore.
    lg = np.fromiter(map(cytoolz.compose(operator.not_, contains('weight>50.0%')), motifs.Context), dtype=np.bool) & \
        np.fromiter(map(contains(*db_names), motifs.Context), dtype=np.bool) & \
        np.fromiter(map(contains('activating'), motifs.Context), dtype=np.bool)
    motifs = motifs.loc[lg,:]
    # We build regulons only using enriched motifs with a NES of 3.0 or higher; we take only directly annotated TFs or TF annotated
    # for an orthologous gene into account; and we only keep regulons with at least 10 genes.
    regulons = list(filter(lambda r: len(r) >= 10, 
        df2regulons(motifs[(motifs.NES >= 3.0) 
        & ((motifs['Annotation'] == 'gene is directly annotated')
        | (motifs['Annotation'].str.startswith('gene is orthologous to')
            & motifs['Annotation'].str.endswith('which is directly annotated for motif')))
    # Rename regulons, i.e. remove suffix.