20260910 143546 Crispr 筛选 Hits 鉴定全流程实现详解

 sgRNA库设计 → FASTQ → 计数 → QC → 拷贝数校正 → Hit Calling → 通路富集 → 可视化
     |           |        |       |         |            |           |          |
  Step 1      Step 2   Step 3  Step 4    Step 5       Step 6      Step 7    Step 8

# 检查文库均匀度
import pandas as pd
import numpy as np

def library_qc(counts_df):
    """文库质量控制指标"""
    out = pd.DataFrame(index=counts_df.columns)
    out['n_detected'] = (counts_df > 0).sum()
    out['pct_zero'] = (counts_df == 0).sum() / len(counts_df) * 100
    out['pct_lowcount'] = (counts_df < 30).sum() / len(counts_df) * 100
    out['gini'] = counts_df.apply(gini_coefficient)
    out['median_count'] = counts_df.median()
    out['p90_p10_ratio'] = counts_df.quantile(0.9) / counts_df.quantile(0.1).replace(0, np.nan)
    return out

def gini_coefficient(x):
    """计算 Gini 指数"""
    x = np.sort(x[x > 0])
    n = len(x)
    cumsum = np.cumsum(x)
    return (2 * np.sum((np.arange(1, n+1)) * x) - (n + 1) * cumsum[-1]) / (n * cumsum[-1])

# fastp 预处理
fastp \
    --input_plasmid_R1.fq.gz --input_plasmid_R2.fq.gz \
    --output_Day0_R1.fq.gz --output_Day0_R2.fq.gz \
    --detect_adapter_for_pe \
    --correction \
    --cut_front --cut_tail \
    --thread 16

# 生成质量报告
fastqc *.fq.gz -t 16 -o qc_reports/
multiqc qc_reports/ -o multiqc_report/

mageck count \
    --list-seq brunello_library.csv \
    --sample-label Plasmid,Day0,Veh_r1,Veh_r2,Drug_r1,Drug_r2 \
    --fastq Plasmid.fq.gz Day0.fq.gz Veh_r1.fq.gz Veh_r2.fq.gz Drug_r1.fq.gz Drug_r2.fq.gz \
    --norm-method median \
    --output-prefix screen \
    --trim-5 AUTO

sgRNA,Sequence,Gene
BRCA1_1,ATGGATTTATCTGCTCTTCG,BRCA1
BRCA1_2,CAGCAGATACTTGATGCATC,BRCA1
NTC_0001,GACGCATCGAATCAATAGCC,NonTargeting_0001

# Replicate 一致性检查
import pandas as pd
import numpy as np
from scipy.stats import pearsonr, spearmanr

def replicate_concordance(counts_df, groups):
    """检查生物学重复的一致性"""
    results = []
    for group_name, samples in groups.items():
        if len(samples) >= 2:
            log_counts = np.log2(counts_df[samples] + 1)
            for i in range(len(samples)):
                for j in range(i+1, len(samples)):
                    r_pearson, p_p = pearsonr(
                        log_counts[samples[i]], log_counts[samples[j]])
                    r_spearman, p_s = spearmanr(
                        log_counts[samples[i]], log_counts[samples[j]])
                    results.append({
                        'group': group_name,
                        'pair': f"{samples[i]} vs {samples[j]}",
                        'pearson_r': r_pearson,
                        'spearman_r': r_spearman
                    })
    return pd.DataFrame(results)

# 使用
groups = {
    'Vehicle': ['Veh_r1', 'Veh_r2'],
    'Drug': ['Drug_r1', 'Drug_r2']
}
concordance = replicate_concordance(counts, groups)
print(concordance)
# 要求: Pearson >= 0.8, Spearman > 0.7

library(CRISPRcleanR)

# 加载文库注释
data(KY_Library_v1.0)

# 归一化 + 计算折变化
norm_counts <- ccr.NormfoldChanges(
    filename = 'counts.txt',
    min_reads = 30,
    EXPname = 'hct116_screen',
    libraryAnnotation = KY_Library_v1.0
)

# 基因组位置校正
cleaned <- ccr.GWclean(
    norm_counts,
    libraryAnnotation = KY_Library_v1.0,
    saveToPath = 'corrected/'
)

# 校正后的计数可用于 MAGeCK 下游分析

python drugz.py \
    -i screen.count.txt \
    -o drugz_results.txt \
    -c Veh_r1,Veh_r2 \
    -x Drug_r1,Drug_r2

mageck test \
    --count-table screen.count.txt \
    --treatment-id Drug_r1,Drug_r2 \
    --control-id Veh_r1,Veh_r2 \
    --norm-method median \
    --output-prefix drug_vs_veh \
    --gene-lfc-method median \
    --sort-criteria pos

# 1. 计算折叠变化
BAGEL.py fc \
    -i screen.count_normalized.txt \
    -o bagel_fc.txt \
    -c Day0 \
    -e Veh_r1,Veh_r2,Drug_r1,Drug_r2

# 2. 计算 Bayes Factor
BAGEL.py bf \
    -i bagel_fc.txt \
    -o bagel_bf.txt \
    -e CEGv2.txt \
    -n NEGv1.txt

# 3. 精确率-召回率曲线
BAGEL.py pr \
    -i bagel_bf.txt \
    -o bagel_pr \
    -e CEGv2.txt \
    -n NEGv1.txt

import pandas as pd

def consensus_hits(drugz_path, mageck_path, fdr_threshold=0.05):
    """跨方法共识 hit 调用"""
    # drugZ hits
    drugz = pd.read_csv(drugz_path, sep='\t')
    drugz_synth = drugz[drugz['synth_fdr'] < fdr_threshold]['gene'].tolist()
    drugz_supp = drugz[drugz['supp_fdr'] < fdr_threshold]['gene'].tolist()

    # MAGeCK hits
    mageck = pd.read_csv(mageck_path, sep='\t')
    mageck_neg = mageck[mageck['neg|fdr'] < fdr_threshold]['id'].tolist()
    mageck_pos = mageck[mageck['pos|fdr'] < fdr_threshold]['id'].tolist()

    # 共识
    consensus_sensitizer = set(drugz_synth) & set(mageck_neg)  # 敏感基因
    consensus_resistance = set(drugz_supp) & set(mageck_pos)   # 抗性基因

    # 单方法 hits（需正交验证）
    drugz_only_synth = set(drugz_synth) - set(mageck_neg)
    mageck_only_neg = set(mageck_neg) - set(drugz_synth)

    return {
        'consensus_sensitizers': consensus_sensitizer,
        'consensus_resistance': consensus_resistance,
        'drugz_only': drugz_only_synth,
        'mageck_only': mageck_only_neg
    }

library(clusterProfiler)

# 对 hit 基因做 GO 富集
go_result <- enrichGO(
    gene = hit_genes,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.2
)

# KEGG 通路富集
kegg_result <- enrichKEGG(
    gene = hit_genes,
    organism = "hsa",
    pAdjustMethod = "BH"
)

# GSEA：用 MAGeCK 排名（neg|score 或 neg|lfc）
library(fgsea)
ranks <- mageck_results$`neg|lfc`
names(ranks) <- mageck_results$id
fgsea_result <- fgsea(
    pathways = kegg_pathways,
    stats = ranks,
    minSize = 15,
    maxSize = 500
)

import matplotlib.pyplot as plt
import numpy as np

def crispr_volcano(gene_summary_path, direction='neg', fdr_threshold=0.05, lfc_threshold=1.0):
    """CRISPR 筛选火山图"""
    df = pd.read_csv(gene_summary_path, sep='\t')
    lfc_col = f'{direction}|lfc'
    fdr_col = f'{direction}|fdr'

    fig, ax = plt.subplots(figsize=(10, 8))
    sig = (df[fdr_col] < fdr_threshold) & (np.abs(df[lfc_col]) > lfc_threshold)

    ax.scatter(df.loc[~sig, lfc_col],
               -np.log10(df.loc[~sig, fdr_col].clip(lower=1e-10)),
               c='lightgray', alpha=0.3, s=8, label='n.s.')
    ax.scatter(df.loc[sig, lfc_col],
               -np.log10(df.loc[sig, fdr_col].clip(lower=1e-10)),
               c='red' if direction == 'neg' else 'blue',
               alpha=0.7, s=20, label=f'FDR<{fdr_threshold}, |LFC|>{lfc_threshold}')

    top = df.loc[sig].nsmallest(20, fdr_col)
    for _, r in top.iterrows():
        ax.annotate(r['id'], (r[lfc_col], -np.log10(max(r[fdr_col], 1e-10))),
                    fontsize=7, ha='center', va='bottom')

    ax.axhline(-np.log10(fdr_threshold), ls='--', c='black', lw=0.5)
    ax.axvline(0, c='gray', lw=0.5)
    ax.axvline(-lfc_threshold, ls=':', c='gray', lw=0.5)
    ax.axvline(lfc_threshold, ls=':', c='gray', lw=0.5)
    ax.set_xlabel(f'{direction} LFC (log2 fold change)')
    ax.set_ylabel('-log10(FDR)')
    ax.legend()
    return fig

 ┌────────────┐   ┌────────────┐   ┌─────────────────┐
 │  Step 1    │   │  Step 2    │   │   Step 3        │
 │ 文库设计    │──►│ FASTQ 预处理│──►│  MAGeCK count   │
 │ Brunello   │   │ fastp      │   │  FASTQ→计数矩阵  │
 └────────────┘   └────────────┘   └────────┬────────┘
                                            │
                    ┌───────────────────────▼──────────────┐
                    │         Step 4: QC                   │
                    │  Gini / replicate / CEGv2 PR-AUC     │
                    └───────────────┬──────────────────────┘
                                    │
                    ┌───────────────▼──────────────────────┐
                    │    Step 5: 拷贝数校正 (癌症细胞系)    │
                    │    CRISPRcleanR / Chronos             │
                    └───────────────┬──────────────────────┘
                                    │
              ┌─────────────────────▼──────────────────────┐
              │            Step 6: Hit Calling              │
              │  ┌─────────┐  ┌──────────┐  ┌───────────┐  │
              │  │ drugZ   │  │ MAGeCK   │  │ BAGEL2    │  │
              │  │ (药物)  │  │ RRA/MLE  │  │ (必需性)  │  │
              │  └────┬────┘  └────┬─────┘  └─────┬─────┘  │
              │       └─────┬──────┘              │        │
              │         共识 hits                │        │
              └─────────────────┬─────────────────────────┘
                                │
              ┌─────────────────▼──────────────┐
              │     Step 7: 通路富集            │
              │  GO / KEGG / GSEA               │
              └─────────────────┬──────────────┘
                                │
              ┌─────────────────▼──────────────┐
              │     Step 8: 可视化              │
              │  Volcano / Heatmap / Rank plot  │
              └────────────────────────────────┘