crispr-screen-analysis

--- name: tooluniverse-crispr-screen-analysis description: Comprehensive CRISPR screen analysis for functional genomics. Analyze pooled or arrayed CRISPR screens (knockout, activation, interference) to identify essential genes, synthetic lethal interactions, and drug targets. Perform sgRNA count processing, gene-level scoring (MAGeCK, BAGEL), quality control, pathway enrichment, and drug target prioritization. Use for CRISPR screen analysis, gene essentiality studies, synthetic lethality detection, functional genomics, drug target validation, or identifying genetic vulnerabilities. --- # ToolUniverse CRISPR Screen Analysis Comprehensive skill for analyzing CRISPR-Cas9 genetic screens to identify essential genes, synthetic lethal interactions, and therapeutic targets through robust statistical analysis and pathway enrichment. ## Overview CRISPR screens enable genome-wide functional genomics by systematically perturbing genes and measuring fitness effects. This skill provides an 8-phase workflow for: - Processing sgRNA count matrices - Quality control and normalization - Gene-level essentiality scoring (MAGeCK-like and BAGEL-like approaches) - Synthetic lethality detection - Pathway enrichment analysis - Drug target prioritization with DepMap integration - Integration with expression and mutation data ## Core Workflow ### Phase 1: Data Import & sgRNA Count Processing **Load sgRNA Count Matrix** ```python import pandas as pd import numpy as np def load_sgrna_counts(counts_file): """ Load sgRNA count matrix from MAGeCK format or generic TSV. Expected format: sgRNA | Gene | Sample1 | Sample2 | Sample3 | ... sgRNA_1 | BRCA1 | 1500 | 1200 | 1100 | ... sgRNA_2 | BRCA1 | 1800 | 1500 | 1400 | ... """ counts = pd.read_csv(counts_file, sep='\t') # Validate required columns required_cols = ['sgRNA', 'Gene'] if not all(col in counts.columns for col in required_cols): raise ValueError(f"Missing required columns: {required_cols}") # Extract sample columns sample_cols = [col for col in counts.columns if col not in ['sgRNA', 'Gene']] # Create count matrix count_matrix = counts[sample_cols].copy() count_matrix.index = counts['sgRNA'] # Gene mapping sgrna_to_gene = dict(zip(counts['sgRNA'], counts['Gene'])) metadata = { 'n_sgrnas': len(counts), 'n_genes': counts['Gene'].nunique(), 'n_samples': len(sample_cols), 'sample_names': sample_cols, 'sgrna_to_gene': sgrna_to_gene } return count_matrix, metadata # Load counts counts, meta = load_sgrna_counts("sgrna_counts.txt") print(f"Loaded {meta['n_sgrnas']} sgRNAs targeting {meta['n_genes']} genes across {meta['n_samples']} samples") ``` **Create Experimental Design Table** ```python def create_design_matrix(sample_names, conditions, timepoints=None): """ Create experimental design linking samples to conditions. Example: Sample | Condition | Timepoint | Replicate T0_rep1 | baseline | 0 | 1 T14_rep1 | treatment | 14 | 1 """ design = pd.DataFrame({ 'Sample': sample_names, 'Condition': conditions }) if timepoints is not None: design['Timepoint'] = timepoints # Auto-detect replicates design['Replicate'] = design.groupby('Condition').cumcount() + 1 return design # Example usage sample_names = ['T0_rep1', 'T0_rep2', 'T14_rep1', 'T14_rep2', 'T14_rep3'] conditions = ['baseline', 'baseline', 'treatment', 'treatment', 'treatment'] design = create_design_matrix(sample_names, conditions) ``` ### Phase 2: Quality Control & Filtering **Assess sgRNA Distribution** ```python def qc_sgrna_distribution(count_matrix, min_reads=30, min_samples=2): """ Quality control for sgRNA distribution. - Remove sgRNAs with low read counts - Check for outlier samples - Assess library representation """ results = {} # 1. Library size per sample library_sizes = count_matrix.sum(axis=0) results['library_sizes'] = library_sizes results['median_library_size'] = library_sizes.median() # 2. Zero-count sgRNAs zero_counts = (count_matrix == 0).sum(axis=1) results['zero_counts'] = zero_counts results['sgrnas_with_zeros'] = (zero_counts > 0).sum() # 3. Low-count sgRNAs (< min_reads in > min_samples) low_count_mask = (count_matrix < min_reads).sum(axis=1) > (len(count_matrix.columns) - min_samples) results['low_count_sgrnas'] = low_count_mask.sum() # 4. Gini coefficient (library skewness) def gini_coefficient(counts): sorted_counts = np.sort(counts) n = len(counts) cumsum = np.cumsum(sorted_counts) return (2 * np.sum((np.arange(1, n+1)) * sorted_counts)) / (n * cumsum[-1]) - (n + 1) / n results['gini_per_sample'] = {col: gini_coefficient(count_matrix[col].values) for col in count_matrix.columns} # 5. Recommend filtering results['filter_recommendation'] = { 'min_reads': min_reads, 'min_samples_above_threshold': min_samples, 'sgrnas_to_remove': low_count_mask.sum() } return results # Run QC qc_results = qc_sgrna_distribution(counts, min_reads=30, min_samples=2) print(f"Library sizes: {qc_results['library_sizes']}") print(f"Low-count sgRNAs to remove: {qc_results['filter_recommendation']['sgrnas_to_remove']}") ``` **Filter Low-Count sgRNAs** ```python def filter_low_count_sgrnas(count_matrix, sgrna_to_gene, min_reads=30, min_samples=2): """ Remove sgRNAs with insufficient read counts. """ # Keep sgRNAs with >= min_reads in >= min_samples keep_mask = (count_matrix >= min_reads).sum(axis=1) >= min_samples filtered_counts = count_matrix[keep_mask].copy() filtered_mapping = {k: v for k, v in sgrna_to_gene.items() if k in filtered_counts.index} print(f"Filtered: {(~keep_mask).sum()} sgRNAs removed, {keep_mask.sum()} retained") return filtered_counts, filtered_mapping # Apply filtering filtered_counts, filtered_mapping = filter_low_count_sgrnas(counts, meta['sgrna_to_gene']) ``` ### Phase 3: Normalization **Library Size Normalization** ```python def normalize_counts(count_matrix, method='median'): """ Normalize sgRNA counts to account for library size differences. Methods: - 'median': Median ratio normalization (like DESeq2) - 'total': Total count normalization (CPM-like) """ if method == 'median': # Calculate geometric mean for each sgRNA across samples pseudo_ref = np.exp(np.log(count_matrix + 1).mean(axis=1)) - 1 # Calculate size factors for each sample size_factors = {} for col in count_matrix.columns: ratios = count_matrix[col] / pseudo_ref ratios = ratios[ratios > 0] # Remove zeros size_factors[col] = ratios.median() # Normalize normalized = count_matrix.div(pd.Series(size_factors), axis=1) elif method == 'total': # CPM-like normalization size_factors = count_matrix.sum(axis=0) / 1e6 normalized = count_matrix.div(size_factors, axis=1) else: raise ValueError(f"Unknown normalization method: {method}") return normalized, size_factors # Normalize norm_counts, size_factors = normalize_counts(filtered_counts, method='median') ``` **Log-Fold Change Calculation** ```python def calculate_lfc(norm_counts, design, control_condition='baseline', treatment_condition='treatment'): """ Calculate log2 fold changes between treatment and control. """ # Get sample names for each condition control_samples = design[design['Condition'] == control_condition]['Sample'].tolist() treatment_samples = design[design['Condition'] == treatment_condition]['Sample'].tolist() # Calculate mean counts control_mean = norm_counts[control_samples].mean(axis=1) treatment_mean = norm_counts[treatment_samples].mean(axis=1) # Log2 fold change (add pseudocount to avoid log(0)) lfc = np.log2((treatment_mean + 1) / (control_mean + 1)) return lfc, control_mean, treatment_mean # Calculate LFC lfc, control_mean, treatment_mean = calculate_lfc(norm_counts, design) ``` ### Phase 4: Gene-Level Scoring (MAGeCK-like) **Aggregate sgRNA Scores to Gene Level** ```python def mageck_gene_scoring(lfc, sgrna_to_gene, method='rra'): """ Gene-level essentiality scoring using MAGeCK-like approach. Methods: - 'rra': Robust Rank Aggregation (identify genes with consistently low-ranking sgRNAs) - 'mean': Simple mean LFC across sgRNAs """ # Create gene-level aggregation gene_lfc = {} for sgrna, gene in sgrna_to_gene.items(): if sgrna in lfc.index: if gene not in gene_lfc: gene_lfc[gene] = [] gene_lfc[gene].append(lfc[sgrna]) if method == 'rra': # Simplified RRA: rank sgRNAs, calculate p-value for each gene # based on whether its sgRNAs are enriched at the top (negative selection) # or bottom (positive selection) # Rank all sgRNAs by LFC ranked_sgrnas = lfc.sort_values() ranks = {sgrna: rank for rank, sgrna in enumerate(ranked_sgrnas.index, 1)} gene_scores = {} for gene, sgrna_list in gene_lfc.items(): # Get ranks for this gene's sgRNAs gene_ranks = [ranks[sgrna] for sgrna in sgrna_list if sgrna in ranks] if len(gene_ranks) > 0: # Use mean rank as score (lower = more essential) gene_scores[gene] = { 'score': np.mean(gene_ranks), 'n_sgrnas': len(gene_ranks), 'mean_lfc': np.mean([lfc[sg] for sg in sgrna_list if sg in lfc.index]) } # Convert to DataFrame gene_df = pd.DataFrame(gene_scores).T gene_df['rank'] = gene_df['score'].rank() elif method == 'mean': # Simple mean LFC gene_df = pd.DataFrame({ gene: { 'mean_lfc': np.mean(sgrna_lfcs), 'n_sgrnas': len(sgrna_lfcs), 'score': np.mean(sgrna_lfcs) } for gene, sgrna_lfcs in gene_lfc.items() }).T # Sort by essentiality (negative LFC = essential) gene_df = gene_df.sort_values('mean_lfc') return gene_df # Gene-level scoring gene_scores = mageck_gene_scoring(lfc, filtered_mapping, method='rra') print(f"Top 10 essential genes:\n{gene_scores.head(10)[['mean_lfc', 'n_sgrnas']]}") ``` **Bayes Factor Scoring (BAGEL-like)** ```python def bagel_bayes_factor(lfc, sgrna_to_gene, essential_genes=None, nonessential_genes=None): """ BAGEL-like Bayes Factor calculation for gene essentiality. Uses reference sets of known essential and non-essential genes to calculate likelihood ratios. """ # Default reference gene sets (core essential genes) if essential_genes is None: essential_genes = ['RPL5', 'RPS6', 'POLR2A', 'PSMC2', 'PSMD14'] # Example if nonessential_genes is None: nonessential_genes = ['AAVS1', 'ROSA26', 'HPRT1'] # Example # Get LFC distributions for reference sets essential_lfc = [lfc[sg] for sg, g in sgrna_to_gene.items() if g in essential_genes and sg in lfc.index] nonessential_lfc = [lfc[sg] for sg, g in sgrna_to_gene.items() if g in nonessential_genes and sg in lfc.index] if len(essential_lfc) < 3 or len(nonessential_lfc) < 3: print("Warning: Insufficient reference genes for BAGEL scoring") return None # Estimate distributions (simplified) essential_mean, essential_std = np.mean(essential_lfc), np.std(essential_lfc) nonessential_mean, nonessential_std = np.mean(nonessential_lfc), np.std(nonessential_lfc) # Calculate Bayes Factor for each gene gene_bf = {} gene_lfc_map = {} for sgrna, gene in sgrna_to_gene.items(): if sgrna in lfc.index: if gene not in gene_lfc_map: gene_lfc_map[gene] = [] gene_lfc_map[gene].append(lfc[sgrna]) for gene, sgrna_lfcs in gene_lfc_map.items(): mean_lfc = np.mean(sgrna_lfcs) # Likelihood under essential distribution from scipy.stats import norm l_essential = norm.pdf(mean_lfc, essential_mean, essential_std) # Likelihood under non-essential distribution l_nonessential = norm.pdf(mean_lfc, nonessential_mean, nonessential_std) # Bayes Factor (avoid division by zero) bf = l_essential / (l_nonessential + 1e-10) gene_bf[gene] = { 'bayes_factor': bf, 'mean_lfc': mean_lfc,