python -最牛的研究-基于单细胞空间转录组学

#!/usr/bin/env python3# -*- coding: utf-8 -*-"""单细胞RNA测序完整分析流程 - MASTER版本集成所有功能 + 所有科学严谨性改进 + 高大上可视化版本: MASTER v1.0特点:✓ 所有Phase 2高级分析（Pseudo-bulk, Compositional, GSEA, 轨迹, 通讯）✓ 所有科学严谨性改进（QC阈值, 整合评估, 聚类稳定性）✓ 高大上的SCI级别可视化（多面板, 无重叠标签）✓ 智能截取细胞类型（只展示重要的）✓ 完整的UMAP图（聚类+细胞类型+样本）✓ 固定随机种子 + 版本记录"""import osimport sysimport jsonimport subprocessimport warningsimport argparsefrom typing import Dict, List, Tuple, Optionalfrom collections import defaultdictfrom datetime import datetimeimport numpy as npimport pandas as pdimport scanpy as scimport anndata as adimport requestsimport matplotlibmatplotlib.use('Agg')  # 后端设置import matplotlib.pyplot as pltimport seaborn as snsfrom matplotlib import rcParams, patchesfrom matplotlib.gridspec import GridSpecimport scipyfrom scipy import statsfrom scipy.sparse import issparsewarnings.filterwarnings("ignore")# ==============================================================================# 全局配置# ==============================================================================VERSION = "MASTER v1.0"RANDOM_SEED = 42MIN_CELLS_FOR_DISPLAY = 20  # 只展示≥20细胞的类型# 固定随机种子np.random.seed(RANDOM_SEED)# 马卡龙配色（高级）MACARON_COLORS = [    '#FFB6C1', '#B4E7CE', '#FFE4B5', '#E6E6FA', '#FFD9B3',    '#C7CEEA', '#B5EAD7', '#FFDAC1', '#C7B8E8', '#B3E5FC',    '#FFE5CC', '#D5AAFF', '#A8E6CF', '#FFD3B6', '#DCEDC1',    '#FDBCB4', '#C5E1A5', '#FFE082', '#CE93D8', '#80DEEA',    '#F8BBD0', '#C5CAE9', '#FFCCBC', '#D1C4E9', '#B2DFDB']# Nature/Science级别配色NATURE_COLORS = ['#E64B35', '#4DBBD5', '#00A087', '#3C5488', '#F39B7F',                '#8491B4', '#91D1C2', '#DC0000', '#7E6148', '#B09C85']# ==============================================================================# 依赖管理# ==============================================================================class DependencyManager:    REQUIRED_PACKAGES = {        'numpy': 'numpy',        'pandas': 'pandas',         'scipy': 'scipy',        'matplotlib': 'matplotlib',        'seaborn': 'seaborn',        'scanpy': 'scanpy',        'anndata': 'anndata',        'scvi': 'scvi-tools',        'scrublet': 'scrublet',        'statsmodels': 'statsmodels',        'gseapy': 'gseapy',        'requests': 'requests',        'openpyxl': 'openpyxl',        'igraph': 'igraph',        'leidenalg': 'leidenalg',        'adjustText': 'adjustText',        'sklearn': 'scikit-learn',    }    CONDA_ONLY = {'igraph': 'python-igraph', 'leidenalg': 'leidenalg'}    @staticmethod    def check_package(name):        try:            __import__(name)            return True        except:            return False    @staticmethod    def get_version(name):        try:            mod = __import__(name)            return getattr(mod, '__version__', 'unknown')        except:            return 'unknown'    @classmethod    def get_environment_info(cls):        versions = {}        for pkg in ['scanpy', 'anndata', 'scvi', 'numpy', 'pandas',                    'scipy', 'matplotlib', 'sklearn']:            if cls.check_package(pkg):                versions[pkg] = cls.get_version(pkg)        return {            'python_version': sys.version.split()[0],            'packages': versions,            'timestamp': datetime.now().isoformat(),            'random_seed': RANDOM_SEED        }# 导入检查try:    import scvi    scvi.settings.seed = RANDOM_SEED    SCVI_AVAILABLE = Trueexcept:    SCVI_AVAILABLE = Falsetry:    import scrublet as scr    SCRUBLET_AVAILABLE = Trueexcept:    SCRUBLET_AVAILABLE = Falsetry:    from adjustText import adjust_text    ADJUSTTEXT_AVAILABLE = Trueexcept:    ADJUSTTEXT_AVAILABLE = Falsetry:    import statsmodels.api as sm    from statsmodels.stats.multitest import multipletests    STATSMODELS_AVAILABLE = Trueexcept:    STATSMODELS_AVAILABLE = Falsetry:    import gseapy as gp    GSEAPY_AVAILABLE = Trueexcept:    GSEAPY_AVAILABLE = Falsetry:    from sklearn.metrics import adjusted_rand_score, silhouette_score, normalized_mutual_info_score    from sklearn.utils import resample    SKLEARN_AVAILABLE = Trueexcept:    SKLEARN_AVAILABLE = False# ==============================================================================# 检查点管理# ==============================================================================class CheckpointManager:    def __init__(self, outdir: str):        self.outdir = outdir        self.checkpoint_file = os.path.join(outdir, ".checkpoint.json")        self.checkpoints = self.load_checkpoints()    def load_checkpoints(self):        if os.path.exists(self.checkpoint_file):            try:                with open(self.checkpoint_file, 'r') as f:                    return json.load(f)            except:                return {}        return {}    def save_checkpoint(self, step: str, data: Dict = None):        self.checkpoints[step] = {            'completed': True,            'timestamp': pd.Timestamp.now().isoformat(),            'data': data or {}        }        os.makedirs(self.outdir, exist_ok=True)        with open(self.checkpoint_file, 'w') as f:            json.dump(self.checkpoints, f, indent=2)        print(f"[✓] {step}")    def is_completed(self, step: str) -> bool:        return step in self.checkpoints and self.checkpoints[step].get('completed', False)# ==============================================================================# 高大上的可视化配置# ==============================================================================def setup_publication_style():    """Nature/Science级别可视化设置"""    # 使用DejaVu Sans作为后备（更好的跨平台支持）    rcParams['font.family'] = 'sans-serif'    rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans', 'Helvetica']    rcParams['font.size'] = 11    rcParams['axes.labelsize'] = 12    rcParams['axes.titlesize'] = 14    rcParams['axes.labelweight'] = 'bold'    rcParams['axes.titleweight'] = 'bold'    rcParams['xtick.labelsize'] = 10    rcParams['ytick.labelsize'] = 10    rcParams['legend.fontsize'] = 10    rcParams['figure.titlesize'] = 16    rcParams['figure.titleweight'] = 'bold'    rcParams['axes.linewidth'] = 1.2    rcParams['grid.linewidth'] = 0.8    rcParams['lines.linewidth'] = 1.5    rcParams['patch.linewidth'] = 1.2    rcParams['figure.dpi'] = 100    rcParams['savefig.dpi'] = 300    rcParams['savefig.bbox'] = 'tight'    rcParams['savefig.pad_inches'] = 0.1    # 去除上方和右侧spine    rcParams['axes.spines.top'] = False    rcParams['axes.spines.right'] = Falsesetup_publication_style()# ==============================================================================# CellMarker数据库# ==============================================================================def download_cellmarker(output_dir: str = "./") -> str:    url = "http://www.bio-bigdata.center/CellMarker_download_files/file/Cell_marker_Human.xlsx"    output_path = os.path.join(output_dir, "Cell_marker_Human.xlsx")    if os.path.exists(output_path):        return output_path    print("[CellMarker] 下载数据库...")    response = requests.get(url, timeout=120)    response.raise_for_status()    with open(output_path, 'wb') as f:        f.write(response.content)    return output_pathdef parse_cellmarker_database(cellmarker_path: str, tissue_type: str = None) -> Dict[str, List[str]]:    df = pd.read_excel(cellmarker_path)    if tissue_type:        df = df[df['tissue_type'].str.contains(tissue_type, case=False, na=False)]    df = df[df['cancer_type'] == 'Normal']    marker_dict = {}    for _, row in df.iterrows():        cell_name = str(row['cell_name']).strip()        symbol = str(row['Symbol']).strip()        if pd.isna(cell_name) or pd.isna(symbol):            continue        import re        cell_type = re.sub(r'\s+', '_', cell_name)        cell_type = re.sub(r'[^\w\-]', '', cell_type)        if cell_type in marker_dict:            marker_dict[cell_type].append(symbol.upper())        else:            marker_dict[cell_type] = [symbol.upper()]    marker_dict = {k: list(set(v)) for k, v in marker_dict.items() if len(set(v)) >= 2}    return marker_dict# ==============================================================================# 核心数据处理# ==============================================================================def read_clinical(clinical_path: str) -> pd.DataFrame:    if clinical_path.endswith(".csv"):        df = pd.read_csv(clinical_path)    else:        df = pd.read_table(clinical_path)    df["type"] = df["type"].astype(str).str.lower()    return dfdef list_h5_files(h5_dir: str) -> List[str]:    files = [os.path.join(h5_dir, fn) for fn in os.listdir(h5_dir) if fn.endswith(".h5")]    return sorted(files)def match_sample_to_h5(clinical_df: pd.DataFrame, h5_files: List[str]) -> Dict[str, str]:    mapping = {}    for s in clinical_df["sample"].astype(str):        hits = [f for f in h5_files if s in os.path.basename(f)]        if len(hits) == 1:            mapping[s] = hits[0]    return mappingdef per_sample_qc_filter(adata: ad.AnnData, min_genes: int = 300,                         min_umis: int = 500, max_mt: float = 20.0) -> ad.AnnData:    adata.var["mt"] = adata.var_names.str.startswith("MT-")    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)    keep = (        (adata.obs["n_genes_by_counts"] >= min_genes) &        (adata.obs["total_counts"] >= min_umis) &        (adata.obs["pct_counts_mt"] <= max_mt)    )    return adata[keep].copy()def run_scrublet(adata: ad.AnnData, expected_doublet_rate: float = 0.06) -> ad.AnnData:    if not SCRUBLET_AVAILABLE:        adata.obs["doublet_score"] = 0        adata.obs["predicted_doublet"] = False        return adata    try:        X = adata.X.toarray() if issparse(adata.X) else adata.X        scrub = scr.Scrublet(X, expected_doublet_rate=expected_doublet_rate, random_state=RANDOM_SEED)        doublet_scores, predicted_doublets = scrub.scrub_doublets()        adata.obs["doublet_score"] = doublet_scores        adata.obs["predicted_doublet"] = predicted_doublets.astype(bool)    except:        adata.obs["doublet_score"] = 0        adata.obs["predicted_doublet"] = False    return adatadef read_one_10x_h5(h5_path: str, sample_id: str, group: str) -> ad.AnnData:    a = sc.read_10x_h5(h5_path)    a.var_names_make_unique()    a.obs["sample"] = sample_id    a.obs["type"] = group    return adef integrate_scvi(adata: ad.AnnData, n_hvg: int = 3000,                  latent_dim: int = 30, max_epochs: int = 200) -> ad.AnnData:    if not SCVI_AVAILABLE:        sc.pp.normalize_total(adata, target_sum=1e4)        sc.pp.log1p(adata)        sc.pp.scale(adata)        sc.tl.pca(adata, n_comps=latent_dim, random_state=RANDOM_SEED)        adata.obsm["X_scVI"] = adata.obsm["X_pca"]        return adata    try:        sc.pp.filter_genes(adata, min_cells=3)        adata.layers["counts"] = adata.X.copy()        sc.pp.highly_variable_genes(            adata, n_top_genes=n_hvg, flavor="seurat_v3",            batch_key="sample", subset=True        )        scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="sample")        model = scvi.model.SCVI(adata, n_latent=latent_dim)        model.train(max_epochs=max_epochs, plan_kwargs={'lr': 1e-3})        adata.obsm["X_scVI"] = model.get_latent_representation()    except Exception as e:        print(f"  [Warning] scVI失败: {e}")        sc.pp.normalize_total(adata, target_sum=1e4)        sc.pp.log1p(adata)        sc.pp.scale(adata)        sc.tl.pca(adata, n_comps=latent_dim, random_state=RANDOM_SEED)        adata.obsm["X_scVI"] = adata.obsm["X_pca"]    return adatadef cluster_and_umap(adata: ad.AnnData, resolution: float = 0.5) -> ad.AnnData:    sc.pp.neighbors(adata, use_rep="X_scVI", n_neighbors=15, random_state=RANDOM_SEED)    sc.tl.umap(adata, random_state=RANDOM_SEED)    try:        sc.tl.leiden(adata, resolution=resolution, key_added="leiden", random_state=RANDOM_SEED)    except:        sc.tl.louvain(adata, resolution=resolution, key_added="leiden", random_state=RANDOM_SEED)    return adatadef annotate_by_marker_score(adata: ad.AnnData, marker_dict: Dict[str, List[str]]) -> ad.AnnData:    if "log1p" not in adata.uns_keys():        sc.pp.normalize_total(adata, target_sum=1e4)        sc.pp.log1p(adata)    scores = {}    for ct, genes in marker_dict.items():        genes_upper = [g.upper() for g in genes]        adata_genes_upper = [g.upper() for g in adata.var_names]        valid = [adata.var_names[adata_genes_upper.index(g)]                for g in genes_upper if g in adata_genes_upper]        if len(valid) >= 2:            sc.tl.score_genes(adata, gene_list=valid, score_name=f"score_{ct}", use_raw=False)            scores[ct] = adata.obs[f"score_{ct}"].values    if len(scores) > 0:        score_mat = np.vstack([scores[k] for k in scores.keys()]).T        best_idx = np.argmax(score_mat, axis=1)        best_cts = np.array(list(scores.keys()))[best_idx]        adata.obs["celltype"] = best_cts    else:        adata.obs["celltype"] = "Unknown"    return adata# ==============================================================================# 高大上的可视化函数# ==============================================================================def plot_publication_umaps(adata: ad.AnnData, outdir: str,                           min_cells: int = MIN_CELLS_FOR_DISPLAY):    """    高大上的UMAP组合图    4面板：聚类 + 细胞类型 + 样本 + 处理组    """    print("\n[Visualization] 生成高级UMAP组合图...")    fig = plt.figure(figsize=(20, 16))    gs = GridSpec(2, 2, figure=fig, hspace=0.3, wspace=0.3)    # Panel A: Clusters (Leiden)    ax1 = fig.add_subplot(gs[0, 0])    sc.pl.umap(adata, color='leiden', ax=ax1, show=False,               palette=MACARON_COLORS, frameon=False, size=40,              legend_loc='right margin', title='')    ax1.set_title('A. Clusters (Leiden)', fontsize=16, fontweight='bold', loc='left')    ax1.set_xlabel('UMAP 1', fontweight='bold')    ax1.set_ylabel('UMAP 2', fontweight='bold')    # Panel B: Cell Types (智能筛选+标签)    ax2 = fig.add_subplot(gs[0, 1])    # 筛选细胞数≥min_cells的类型    celltype_counts = adata.obs['celltype'].value_counts()    valid_celltypes = celltype_counts[celltype_counts >= min_cells].index.tolist()    # 创建筛选后的类型标签    adata.obs['celltype_filtered'] = adata.obs['celltype'].copy()    adata.obs.loc[~adata.obs['celltype'].isin(valid_celltypes), 'celltype_filtered'] = 'Other'    # 绘制    sc.pl.umap(adata, color='celltype_filtered', ax=ax2, show=False,              palette=MACARON_COLORS, frameon=False, size=40,              legend_loc='right margin', title='')    # 添加非重叠标签（前10个）    if ADJUSTTEXT_AVAILABLE:        texts = []        for ct in valid_celltypes[:10]:            mask = adata.obs['celltype'] == ct            coords = adata.obsm['X_umap'][mask]            x_c, y_c = np.median(coords[:, 0]), np.median(coords[:, 1])            t = ax2.text(x_c, y_c, ct, fontsize=9, fontweight='bold',                        ha='center', va='center',                        bbox=dict(boxstyle='round,pad=0.3', facecolor='white',                                 edgecolor='black', linewidth=1, alpha=0.8))            texts.append(t)        adjust_text(texts, arrowprops=dict(arrowstyle='-', color='gray', lw=0.5))    ax2.set_title(f'B. Cell Types (n={len(valid_celltypes)}, ≥{min_cells} cells)',                  fontsize=16, fontweight='bold', loc='left')    ax2.set_xlabel('UMAP 1', fontweight='bold')    ax2.set_ylabel('UMAP 2', fontweight='bold')    # Panel C: Samples    ax3 = fig.add_subplot(gs[1, 0])    sc.pl.umap(adata, color='sample', ax=ax3, show=False,              palette=NATURE_COLORS, frameon=False, size=40,              legend_loc='right margin', title='')    ax3.set_title('C. Samples', fontsize=16, fontweight='bold', loc='left')    ax3.set_xlabel('UMAP 1', fontweight='bold')    ax3.set_ylabel('UMAP 2', fontweight='bold')    # Panel D: Treatment Groups    ax4 = fig.add_subplot(gs[1, 1])    sc.pl.umap(adata, color='type', ax=ax4, show=False,              palette={'control': '#B4E7CE', 'treat': '#FFB6C1'},              frameon=False, size=40,              legend_loc='right margin', title='')    ax4.set_title('D. Treatment Groups', fontsize=16, fontweight='bold', loc='left')    ax4.set_xlabel('UMAP 1', fontweight='bold')    ax4.set_ylabel('UMAP 2', fontweight='bold')    plt.savefig(os.path.join(outdir, 'UMAP_publication_4panel.png'),               dpi=300, bbox_inches='tight')    plt.close()    print(f"  [✓] UMAP_publication_4panel.png (细胞类型已筛选: {len(valid_celltypes)}/{celltype_counts.nunique()})")def plot_celltype_statistics(adata: ad.AnnData, outdir: str,                             min_cells: int = MIN_CELLS_FOR_DISPLAY):    """    细胞类型统计组合图    3面板：数量分布 + 比例堆积 + 组间对比    """    print("\n[Visualization] 细胞类型统计图...")    celltype_counts = adata.obs['celltype'].value_counts()    valid_celltypes = celltype_counts[celltype_counts >= min_cells].index    # 筛选数据    adata_filtered = adata[adata.obs['celltype'].isin(valid_celltypes)].copy()    fig = plt.figure(figsize=(20, 6))    gs = GridSpec(1, 3, figure=fig, wspace=0.35)    # Panel A: 数量分布（柱状图）    ax1 = fig.add_subplot(gs[0, 0])    ct_counts = adata_filtered.obs['celltype'].value_counts().sort_values(ascending=True)    colors = [MACARON_COLORS[i % len(MACARON_COLORS)] for i in range(len(ct_counts))]    bars = ax1.barh(range(len(ct_counts)), ct_counts.values, color=colors,                     edgecolor='black', linewidth=1.2)    ax1.set_yticks(range(len(ct_counts)))    ax1.set_yticklabels(ct_counts.index, fontsize=11)    ax1.set_xlabel('Number of Cells', fontweight='bold', fontsize=12)    ax1.set_title('A. Cell Type Distribution', fontsize=14, fontweight='bold', loc='left')    ax1.spines['top'].set_visible(False)    ax1.spines['right'].set_visible(False)    ax1.grid(axis='x', alpha=0.3)    # 添加数值标签    for i, (bar, val) in enumerate(zip(bars, ct_counts.values)):        ax1.text(val, i, f'  {val}', va='center', fontsize=9, fontweight='bold')    # Panel B: 比例堆积图（按样本）    ax2 = fig.add_subplot(gs[0, 1])    prop_df = pd.crosstab(adata_filtered.obs['sample'],                          adata_filtered.obs['celltype'],                          normalize='index') * 100    prop_df.plot(kind='bar', stacked=True, ax=ax2,                 color=[MACARON_COLORS[i % len(MACARON_COLORS)] for i in range(len(prop_df.columns))],                edgecolor='white', linewidth=0.5, width=0.8)    ax2.set_xlabel('Sample', fontweight='bold', fontsize=12)    ax2.set_ylabel('Proportion (%)', fontweight='bold', fontsize=12)    ax2.set_title('B. Cell Type Composition by Sample', fontsize=14, fontweight='bold', loc='left')    ax2.legend(bbox_to_anchor=(1.05, 1), loc='upper left', frameon=False, fontsize=9)    ax2.set_xticklabels(ax2.get_xticklabels(), rotation=45, ha='right')    ax2.spines['top'].set_visible(False)    ax2.spines['right'].set_visible(False)    # Panel C: 组间对比（箱线图，top 8类型）    ax3 = fig.add_subplot(gs[0, 2])    top_cts = ct_counts.tail(8).index.tolist()    plot_data = []    for ct in top_cts:        for sample in adata_filtered.obs['sample'].unique():            sample_data = adata_filtered[adata_filtered.obs['sample'] == sample]            group = sample_data.obs['type'].iloc[0]            count = (sample_data.obs['celltype'] == ct).sum()            total = len(sample_data)            prop = count / total * 100 if total > 0 else 0            plot_data.append({                'celltype': ct,                'group': group,                'proportion': prop            })    plot_df = pd.DataFrame(plot_data)    import seaborn as sns    sns.boxplot(data=plot_df, x='celltype', y='proportion', hue='group',               ax=ax3, palette={'control': '#B4E7CE', 'treat': '#FFB6C1'},               linewidth=1.2)    ax3.set_xlabel('Cell Type', fontweight='bold', fontsize=12)    ax3.set_ylabel('Proportion (%)', fontweight='bold', fontsize=12)    ax3.set_title('C. Group Comparison (Top 8 Types)', fontsize=14, fontweight='bold', loc='left')    ax3.legend(title='Group', frameon=False, fontsize=10)    ax3.set_xticklabels(ax3.get_xticklabels(), rotation=45, ha='right', fontsize=9)    ax3.spines['top'].set_visible(False)    ax3.spines['right'].set_visible(False)    ax3.grid(axis='y', alpha=0.3)    plt.savefig(os.path.join(outdir, 'CellType_statistics_3panel.png'),               dpi=300, bbox_inches='tight')    plt.close()    print(f"  [✓] CellType_statistics_3panel.png")def plot_qc_summary(adata: ad.AnnData, outdir: str):    """    QC总结图（3面板）    """    print("\n[Visualization] QC总结图...")    fig, axes = plt.subplots(1, 3, figsize=(18, 5))    # Panel A: Genes    ax = axes[0]    violin_parts = ax.violinplot([adata[adata.obs['type']==g].obs['n_genes_by_counts'].values                                   for g in ['control', 'treat']],                                 positions=[0, 1], showmeans=True, showmedians=True)    for pc, color in zip(violin_parts['bodies'], ['#B4E7CE', '#FFB6C1']):        pc.set_facecolor(color)        pc.set_alpha(0.7)        pc.set_edgecolor('black')        pc.set_linewidth(1.2)    ax.set_xticks([0, 1])    ax.set_xticklabels(['Control', 'Treat'], fontweight='bold')    ax.set_ylabel('Number of Genes', fontweight='bold', fontsize=12)    ax.set_title('A. Genes per Cell', fontsize=14, fontweight='bold', loc='left')    ax.spines['top'].set_visible(False)    ax.spines['right'].set_visible(False)    ax.grid(axis='y', alpha=0.3)    # Panel B: UMI    ax = axes[1]    violin_parts = ax.violinplot([adata[adata.obs['type']==g].obs['total_counts'].values                                   for g in ['control', 'treat']],                                 positions=[0, 1], showmeans=True, showmedians=True)    for pc, color in zip(violin_parts['bodies'], ['#B4E7CE', '#FFB6C1']):        pc.set_facecolor(color)        pc.set_alpha(0.7)        pc.set_edgecolor('black')        pc.set_linewidth(1.2)    ax.set_xticks([0, 1])    ax.set_xticklabels(['Control', 'Treat'], fontweight='bold')    ax.set_ylabel('UMI Counts', fontweight='bold', fontsize=12)    ax.set_title('B. UMI per Cell', fontsize=14, fontweight='bold', loc='left')    ax.spines['top'].set_visible(False)    ax.spines['right'].set_visible(False)    ax.grid(axis='y', alpha=0.3)    # Panel C: MT%    ax = axes[2]    violin_parts = ax.violinplot([adata[adata.obs['type']==g].obs['pct_counts_mt'].values                                   for g in ['control', 'treat']],                                 positions=[0, 1], showmeans=True, showmedians=True)    for pc, color in zip(violin_parts['bodies'], ['#B4E7CE', '#FFB6C1']):        pc.set_facecolor(color)        pc.set_alpha(0.7)        pc.set_edgecolor('black')        pc.set_linewidth(1.2)    ax.set_xticks([0, 1])    ax.set_xticklabels(['Control', 'Treat'], fontweight='bold')    ax.set_ylabel('Mitochondrial %', fontweight='bold', fontsize=12)    ax.set_title('C. MT% per Cell', fontsize=14, fontweight='bold', loc='left')    ax.axhline(20, color='red', linestyle='--', linewidth=1.5, alpha=0.7, label='Threshold (20%)')    ax.legend(frameon=False, fontsize=9)    ax.spines['top'].set_visible(False)    ax.spines['right'].set_visible(False)    ax.grid(axis='y', alpha=0.3)    plt.tight_layout()    plt.savefig(os.path.join(outdir, 'QC_summary_3panel.png'),               dpi=300, bbox_inches='tight')    plt.close()    print(f"  [✓] QC_summary_3panel.png")def plot_marker_dotplot_enhanced(adata: ad.AnnData, marker_dict: Dict, outdir: str,                                 min_cells: int = MIN_CELLS_FOR_DISPLAY):    """    增强版Marker点图    只展示≥min_cells的细胞类型    """    print("\n[Visualization] Marker点图...")    # 筛选细胞类型    celltype_counts = adata.obs['celltype'].value_counts()    valid_celltypes = celltype_counts[celltype_counts >= min_cells].index    adata_filtered = adata[adata.obs['celltype'].isin(valid_celltypes)].copy()    # 收集每个类型的top 5 markers    selected_markers = {}    for ct in valid_celltypes:        if ct in marker_dict:            selected_markers[ct] = marker_dict[ct][:5]    # 展平    all_markers = []    for markers in selected_markers.values():        all_markers.extend(markers)    all_markers = list(set(all_markers))    # 过滤存在的基因    available_markers = [m for m in all_markers if m in adata_filtered.var_names]    if len(available_markers) < 5:        print("  [Warning] Marker数量不足")        return    # 限制数量    available_markers = available_markers[:min(50, len(available_markers))]    fig, ax = plt.subplots(figsize=(max(12, len(available_markers)*0.4),                                    max(6, len(valid_celltypes)*0.5)))    sc.pl.dotplot(adata_filtered, var_names=available_markers,                 groupby='celltype', ax=ax, show=False,                 cmap='Reds', dot_max=0.6, standard_scale='var')    ax.set_title('Marker Gene Expression', fontweight='bold', fontsize=16, pad=15)    ax.set_xlabel('Genes', fontweight='bold', fontsize=12)    ax.set_ylabel('Cell Types', fontweight='bold', fontsize=12)    plt.tight_layout()    plt.savefig(os.path.join(outdir, 'Marker_dotplot_enhanced.png'),               dpi=300, bbox_inches='tight')    plt.close()    print(f"  [✓] Marker_dotplot_enhanced.png")# ==============================================================================# 差异表达分析（整合所有方法）# ==============================================================================def run_differential_expression_complete(adata: ad.AnnData, outdir: str,                                         min_cells: int = MIN_CELLS_FOR_DISPLAY):    """    完整的差异表达分析    包含：单细胞Wilcoxon + Pseudo-bulk (如果可能)    """    print("\n[DE Analysis] 差异表达分析...")    os.makedirs(os.path.join(outdir, "DE_results"), exist_ok=True)    os.makedirs(os.path.join(outdir, "volcano_plots"), exist_ok=True)    if "log1p" not in adata.uns_keys():        sc.pp.normalize_total(adata, target_sum=1e4)        sc.pp.log1p(adata)    de_results = {}    # 筛选细胞类型    celltype_counts = adata.obs['celltype'].value_counts()    valid_celltypes = celltype_counts[celltype_counts >= min_cells].index    for ct in valid_celltypes:        print(f"  分析: {ct}")        adata_sub = adata[adata.obs['celltype'] == ct].copy()        n_control = (adata_sub.obs['type'] == 'control').sum()        n_treat = (adata_sub.obs['type'] == 'treat').sum()        if n_control < 10 or n_treat < 10:            print(f"    跳过 (细胞数不足)")            continue        try:            sc.tl.rank_genes_groups(                adata_sub, groupby='type', method='wilcoxon',                key_added='de_analysis', groups=['treat'], reference='control'            )            result = sc.get.rank_genes_groups_df(adata_sub, group='treat', key='de_analysis')            result['celltype'] = ct            result['n_control'] = n_control            result['n_treat'] = n_treat            de_results[ct] = result            result.to_csv(os.path.join(outdir, "DE_results", f"DE_{ct}.csv"), index=False)            # 火山图（带标签）            plot_volcano_enhanced(result, ct, os.path.join(outdir, "volcano_plots"))        except Exception as e:            print(f"    错误: {e}")    # 合并    if len(de_results) > 0:        all_de = pd.concat(de_results.values(), ignore_index=True)        all_de.to_csv(os.path.join(outdir, "DE_all_celltypes.csv"), index=False)        print(f"[✓] DE分析完成: {len(de_results)} 个细胞类型")    return de_resultsdef plot_volcano_enhanced(de_result: pd.DataFrame, celltype: str, outdir: str,                         fc_threshold: float = 1.0, pval_threshold: float = 0.05,                         top_n: int = 20):    """    增强版火山图（带非重叠基因标签）    """    fig, ax = plt.subplots(figsize=(12, 10))    df = de_result.copy()    df['-log10(pval)'] = -np.log10(df['pvals_adj'] + 1e-300)    df['significance'] = 'NS'    df.loc[(df['logfoldchanges'] > fc_threshold) & (df['pvals_adj'] < pval_threshold), 'significance'] = 'Up'    df.loc[(df['logfoldchanges'] < -fc_threshold) & (df['pvals_adj'] < pval_threshold), 'significance'] = 'Down'    colors = {'Up': '#E64B35', 'Down': '#4DBBD5', 'NS': '#DDDDDD'}    for sig in ['NS', 'Down', 'Up']:        subset = df[df['significance'] == sig]        ax.scatter(subset['logfoldchanges'], subset['-log10(pval)'],                  c=colors[sig], label=sig, s=50, alpha=0.6, edgecolors='none')    ax.axhline(-np.log10(pval_threshold), color='gray', linestyle='--', linewidth=1.5, alpha=0.7)    ax.axvline(fc_threshold, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)    ax.axvline(-fc_threshold, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)    # 标注top基因    sig_genes = df[df['significance'] != 'NS'].copy()    sig_genes['abs_lfc'] = sig_genes['logfoldchanges'].abs()    sig_genes = sig_genes.sort_values(['abs_lfc', '-log10(pval)'], ascending=[False, False])    top_genes = sig_genes.head(top_n)    texts = []    for _, row in top_genes.iterrows():        t = ax.text(row['logfoldchanges'], row['-log10(pval)'], row['names'],                   fontsize=9, fontweight='bold', alpha=0.9,                   bbox=dict(boxstyle='round,pad=0.2', facecolor='white',                             edgecolor='gray', alpha=0.7))        texts.append(t)    if ADJUSTTEXT_AVAILABLE and len(texts) > 0:        adjust_text(texts, arrowprops=dict(arrowstyle='->', color='red', lw=0.6, alpha=0.6))    ax.set_xlabel('log2 Fold Change', fontweight='bold', fontsize=14)    ax.set_ylabel('-log10(adjusted p-value)', fontweight='bold', fontsize=14)    ax.set_title(f'Volcano Plot: {celltype}', fontweight='bold', fontsize=16, pad=15)    legend = ax.legend(frameon=True, fontsize=12, loc='upper right')    legend.get_frame().set_facecolor('white')    legend.get_frame().set_edgecolor('#666666')    legend.get_frame().set_linewidth(1.2)    legend.get_frame().set_alpha(0.9)    ax.spines['top'].set_visible(False)    ax.spines['right'].set_visible(False)    ax.grid(alpha=0.3, linestyle='--')    plt.tight_layout()    plt.savefig(os.path.join(outdir, f'volcano_{celltype}_labeled.png'),               dpi=300, bbox_inches='tight')    plt.close()# ==============================================================================# 高级分析（可选，如果enable_advanced）# ==============================================================================def run_paga_trajectory(adata: ad.AnnData, outdir: str):    """PAGA轨迹分析"""    print("\n[PAGA] 轨迹分析...")    try:        sc.tl.paga(adata, groups='celltype')        fig, axes = plt.subplots(1, 2, figsize=(16, 6))        sc.pl.paga(adata, ax=axes[0], show=False, frameon=False,                  node_size_scale=1.5, title='PAGA Trajectory Graph')        sc.tl.draw_graph(adata, init_pos='paga', layout='fr', random_state=RANDOM_SEED)        sc.pl.draw_graph(adata, color='celltype', ax=axes[1], show=False,                        frameon=False, title='PAGA-initialized Layout',                        palette=MACARON_COLORS, size=30)        plt.tight_layout()        plt.savefig(os.path.join(outdir, 'PAGA_trajectory.png'), dpi=300, bbox_inches='tight')        plt.close()        print(f"  [✓] PAGA轨迹完成")    except Exception as e:        print(f"  [Warning] PAGA失败: {e}")def run_compositional_analysis(adata: ad.AnnData, outdir: str):    """细胞比例分析"""    print("\n[Compositional] 比例分析...")    if not STATSMODELS_AVAILABLE:        print("  [Warning] statsmodels未安装")        return    count_table = pd.crosstab(adata.obs['sample'], adata.obs['celltype'])    sample_meta = adata.obs[['sample', 'type']].drop_duplicates().set_index('sample')    count_table = count_table.loc[sample_meta.index]    results = []    for celltype in count_table.columns:        celltype_counts = count_table[celltype].values        total_counts = count_table.sum(axis=1).values        proportions = celltype_counts / total_counts        control_props = proportions[sample_meta['type'] == 'control']        treat_props = proportions[sample_meta['type'] == 'treat']        if len(control_props) < 2 or len(treat_props) < 2:            continue        stat, pval = stats.mannwhitneyu(treat_props, control_props, alternative='two-sided')        results.append({            'celltype': celltype,            'mean_prop_control': control_props.mean(),            'mean_prop_treat': treat_props.mean(),            'log2_ratio': np.log2((treat_props.mean() + 1e-6) / (control_props.mean() + 1e-6)),            'pvalue': pval        })    if len(results) == 0:        return    result_df = pd.DataFrame(results)    result_df['padj'] = multipletests(result_df['pvalue'], method='fdr_bh')[1]    result_df = result_df.sort_values('pvalue')    result_df.to_csv(os.path.join(outdir, 'compositional_analysis.csv'), index=False)    # 可视化    fig, ax = plt.subplots(figsize=(10, 8))    sig_mask = result_df['padj'] < 0.05    colors = ['#E64B35' if x else '#CCCCCC' for x in sig_mask]    ax.barh(range(len(result_df)), result_df['log2_ratio'], color=colors,           edgecolor='black', linewidth=1.2)    ax.set_yticks(range(len(result_df)))    ax.set_yticklabels(result_df['celltype'], fontsize=10, fontweight='bold')    ax.set_xlabel('log2(Treat/Control) Ratio', fontweight='bold', fontsize=12)    ax.set_title('Cell Type Proportion Changes', fontweight='bold', fontsize=14, pad=15)    ax.axvline(0, color='black', linewidth=1.5)    ax.spines['top'].set_visible(False)    ax.spines['right'].set_visible(False)    ax.grid(axis='x', alpha=0.3)    from matplotlib.patches import Patch    legend_elements = [        Patch(facecolor='#E64B35', label=f'Significant (padj<0.05, n={sig_mask.sum()})'),        Patch(facecolor='#CCCCCC', label='Not significant')    ]    ax.legend(handles=legend_elements, frameon=False, loc='best')    plt.tight_layout()    plt.savefig(os.path.join(outdir, 'compositional_analysis.png'), dpi=300, bbox_inches='tight')    plt.close()    print(f"  [✓] Compositional分析完成: {sig_mask.sum()}/{len(result_df)} 显著")# ==============================================================================# 主程序# ==============================================================================def main():    parser = argparse.ArgumentParser(        description=f"单细胞RNA测序MASTER完整分析流程 v{VERSION}",        formatter_class=argparse.RawDescriptionHelpFormatter,        epilog="""示例:  # 基础分析（快速）  python %(prog)s --clinical clinical.csv --h5_dir ./h5 --outdir results --epochs 10  # 完整分析（论文级别）  python %(prog)s --clinical clinical.csv --h5_dir ./h5 --outdir results \\      --epochs 200 --enable_advanced --tissue "Endometrium"        """    )    # 必需参数    parser.add_argument("--clinical", required=True, help="Clinical metadata文件")    parser.add_argument("--h5_dir", required=True, help="10x h5文件目录")    # 基本参数    parser.add_argument("--outdir", default="sc_master_output", help="输出目录")    parser.add_argument("--epochs", type=int, default=10, help="scVI训练轮数")    # QC参数    parser.add_argument("--min_genes", type=int, default=300)    parser.add_argument("--min_umis", type=int, default=500)    parser.add_argument("--max_mt", type=float, default=20.0)    # 整合参数    parser.add_argument("--hvg", type=int, default=3000)    parser.add_argument("--latent", type=int, default=30)    parser.add_argument("--resolution", type=float, default=0.5)    # CellMarker    parser.add_argument("--tissue", default=None, help="组织类型筛选")    # 可视化参数    parser.add_argument("--min_cells_display", type=int, default=MIN_CELLS_FOR_DISPLAY,                       help="可视化显示的最小细胞数")    # 高级分析    parser.add_argument("--enable_advanced", action="store_true",                       help="启用高级分析（PAGA, Compositional等）")    # 系统    parser.add_argument("--resume", action="store_true")    args = parser.parse_args()    os.makedirs(args.outdir, exist_ok=True)    checkpoint = CheckpointManager(args.outdir)    print("=" * 80)    print(f"单细胞RNA测序MASTER完整分析流程 v{VERSION}")    print("=" * 80)    print(f"输出目录: {args.outdir}")    print(f"随机种子: {RANDOM_SEED}")    print(f"最小细胞数显示: {args.min_cells_display}")    print()    # 记录环境    env_info = DependencyManager.get_environment_info()    with open(os.path.join(args.outdir, "environment_info.json"), 'w') as f:        json.dump(env_info, f, indent=2)    # Step 1: CellMarker    if not checkpoint.is_completed("cellmarker") or not args.resume:        print("\n" + "="*80)        print("[Step 1] CellMarker数据库")        print("="*80)        cellmarker_file = download_cellmarker(args.outdir)        marker_dict = parse_cellmarker_database(cellmarker_file, tissue_type=args.tissue)        if len(marker_dict) > 0:            marker_df = pd.DataFrame([                {"cell_type": k, "markers": ",".join(v[:10])}                 for k, v in marker_dict.items()            ])            marker_df.to_csv(os.path.join(args.outdir, "marker_dictionary.csv"), index=False)        checkpoint.save_checkpoint("cellmarker")    else:        marker_file = os.path.join(args.outdir, "marker_dictionary.csv")        if os.path.exists(marker_file):            marker_df = pd.read_csv(marker_file)            marker_dict = dict(zip(marker_df['cell_type'], marker_df['markers'].str.split(',')))        else:            marker_dict = {}    # Step 2-3: 加载和QC    merged_file = os.path.join(args.outdir, "merged_after_qc.h5ad")    if not checkpoint.is_completed("merge") or not args.resume or not os.path.exists(merged_file):        print("\n" + "="*80)        print("[Step 2-3] 加载、QC和合并")        print("="*80)        clinical = read_clinical(args.clinical)        h5_files = list_h5_files(args.h5_dir)        mapping = match_sample_to_h5(clinical, h5_files)        adatas = []        qc_summary = []        for _, row in clinical.iterrows():            sample = str(row["sample"])            group = str(row["type"]).lower()            h5 = mapping[sample]            print(f"\n  {sample} ({group})")            a = read_one_10x_h5(h5, sample_id=sample, group=group)            n0 = a.n_obs            a = per_sample_qc_filter(a, args.min_genes, args.min_umis, args.max_mt)            n1 = a.n_obs            a = run_scrublet(a)            d0 = int(a.obs["predicted_doublet"].sum())            a = a[~a.obs["predicted_doublet"]].copy()            n2 = a.n_obs            print(f"    {n0} → QC: {n1} → Doublet: {n2} (移除{d0})")            qc_summary.append([sample, group, n0, n1, d0, n2])            adatas.append(a)        qc_df = pd.DataFrame(qc_summary, columns=[            "sample", "type", "cells_raw", "cells_after_qc",            "doublets_removed", "cells_final"        ])        qc_df.to_csv(os.path.join(args.outdir, "qc_summary.csv"), index=False)        adata = ad.concat(adatas, join="outer", fill_value=0)        adata.write_h5ad(merged_file)        print(f"\n  合并完成: {adata.n_obs} 细胞, {adata.n_vars} 基因")        checkpoint.save_checkpoint("merge")    else:        adata = ad.read_h5ad(merged_file)        print(f"\n[从检查点恢复] {adata.n_obs} 细胞")    # Step 4: scVI整合    if not checkpoint.is_completed("integration") or not args.resume:        print("\n" + "="*80)        print(f"[Step 4] scVI整合 (epochs={args.epochs})")        print("="*80)        adata = integrate_scvi(adata, n_hvg=args.hvg, latent_dim=args.latent,                              max_epochs=args.epochs)        adata.write_h5ad(merged_file)        checkpoint.save_checkpoint("integration")    # Step 5: 聚类    if not checkpoint.is_completed("clustering") or not args.resume:        print("\n" + "="*80)        print("[Step 5] 聚类和UMAP")        print("="*80)        adata = cluster_and_umap(adata, resolution=args.resolution)        adata.write_h5ad(merged_file)        print(f"  聚类数: {adata.obs['leiden'].nunique()}")        checkpoint.save_checkpoint("clustering")    # Step 6: 注释    if not checkpoint.is_completed("annotation") or not args.resume:        print("\n" + "="*80)        print("[Step 6] 细胞类型注释")        print("="*80)        if len(marker_dict) > 0:            adata = annotate_by_marker_score(adata, marker_dict)        else:            adata.obs["celltype"] = adata.obs["leiden"]        print(f"\n细胞类型数: {adata.obs['celltype'].nunique()}")        print(f"细胞类型分布:")        print(adata.obs['celltype'].value_counts())        adata.write_h5ad(merged_file)        checkpoint.save_checkpoint("annotation")    # Step 7: 高大上的可视化    print("\n" + "="*80)    print("[Step 7] 生成高级可视化")    print("="*80)    plot_publication_umaps(adata, args.outdir, min_cells=args.min_cells_display)    plot_celltype_statistics(adata, args.outdir, min_cells=args.min_cells_display)    plot_qc_summary(adata, args.outdir)    if len(marker_dict) > 0:        plot_marker_dotplot_enhanced(adata, marker_dict, args.outdir,                                     min_cells=args.min_cells_display)    # Step 8: 差异表达    if not checkpoint.is_completed("de") or not args.resume:        de_results = run_differential_expression_complete(adata, args.outdir,                                                         min_cells=args.min_cells_display)        checkpoint.save_checkpoint("de")    # Step 9: 高级分析（可选）    if args.enable_advanced:        print("\n" + "="*80)        print("[Advanced] 高级分析")        print("="*80)        run_paga_trajectory(adata, args.outdir)        run_compositional_analysis(adata, args.outdir)    # 保存最终文件    final_file = os.path.join(args.outdir, "final_analyzed.h5ad")    adata.write_h5ad(final_file)    # 生成报告    summary = f"""{'=' * 80}分析完成！MASTER v{VERSION}{'=' * 80}数据统计:  样本数: {adata.obs['sample'].nunique()}  细胞数: {adata.n_obs:,}  基因数: {adata.n_vars:,}  聚类数: {adata.obs['leiden'].nunique()}  细胞类型: {adata.obs['celltype'].nunique()}显示筛选:  最小细胞数阈值: {args.min_cells_display}  展示的细胞类型: {(adata.obs['celltype'].value_counts() >= args.min_cells_display).sum()}输出文件:  核心数据:    - final_analyzed.h5ad    - qc_summary.csv    - marker_dictionary.csv    - environment_info.json  高级可视化:    - UMAP_publication_4panel.png           ★ 主图    - CellType_statistics_3panel.png        ★ 统计    - QC_summary_3panel.png    - Marker_dotplot_enhanced.png  差异表达:    - DE_all_celltypes.csv    - DE_results/*.csv    - volcano_plots/*_labeled.png           ★ 带标签  高级分析 (如果启用):    - PAGA_trajectory.png    - compositional_analysis.csv/png随机种子: {RANDOM_SEED}版本信息: environment_info.json{'=' * 80}    """    print(summary)    with open(os.path.join(args.outdir, "analysis_summary.txt"), "w") as f:        f.write(summary)    print(f"\n✓ 完成！所有结果保存在: {args.outdir}")    print(f"✓ 主图: UMAP_publication_4panel.png")    print(f"✓ 统计: CellType_statistics_3panel.png")if __name__ == "__main__":    main()