160_碱基编辑检测分析¶

设计sgRNA(靶向目标碱基在编辑窗口内)
    ↓
转染碱基编辑器+sgRNA到细胞
    ↓
培养48-72h
    ↓
提取基因组DNA
    ↓
PCR扩增靶位点(~200-300bp)
    ↓
深度测序(NGS, 通常>10000x)
    ↓
CRISPResso2/BE-Analyzer分析
    ↓
编辑效率、纯度、bystander评估

# === 安装CRISPResso2 ===
pip install CRISPResso2
# 或 conda install -c bioconda crispresso2

# === 基本碱基编辑分析 ===
CRISPResso --fastq_r1 sample_R1.fq.gz \
  --fastq_r2 sample_R2.fq.gz \
  --amplicon_seq ATCGATCGATCGATCGATCGATCGATCGATCGATCGATCG \
  --guide_seq ATCGATCGATCGATCGATCG \
  --base_editor_output \
  --conversion_nuc_from C \
  --conversion_nuc_to T \
  --output_folder crispresso_output/ \
  --name sample1 \
  --min_frequency_alleles_around_cut_to_plot 0.01 \
  --plot_window_size 20

# 参数说明:
# --amplicon_seq: 扩增子的参考序列(无PAM)
# --guide_seq: sgRNA序列(不含PAM, 与amplicon_seq方向一致)
# --base_editor_output: 开启碱基编辑模式
# --conversion_nuc_from/to: CBE用C→T, ABE用A→G
# --plot_window_size: 分析窗口大小

# === ABE分析 ===
CRISPResso --fastq_r1 sample_R1.fq.gz \
  --fastq_r2 sample_R2.fq.gz \
  --amplicon_seq REFERENCE_SEQUENCE \
  --guide_seq GUIDE_20NT \
  --base_editor_output \
  --conversion_nuc_from A \
  --conversion_nuc_to G \
  --output_folder abe_output/ \
  --name abe_sample1

# === 批量分析 ===
CRISPRessoBatch --batch_settings batch.txt \
  --output_folder batch_output/ \
  --base_editor_output \
  --conversion_nuc_from C \
  --conversion_nuc_to T

# batch.txt格式:
# name  fastq_r1  fastq_r2  amplicon_seq  guide_seq
# sample1  s1_R1.fq.gz  s1_R2.fq.gz  AMPLICON_SEQ  GUIDE_SEQ
# sample2  s2_R1.fq.gz  s2_R2.fq.gz  AMPLICON_SEQ  GUIDE_SEQ

# BE-Analyzer: http://www.rgenome.net/be-analyzer/
# 也有命令行版本

# === 本地Python实现碱基编辑定量 ===

import pandas as pd
import numpy as np
from collections import Counter
from Bio import SeqIO, pairwise2

def analyze_base_editing(fastq_file, reference, guide_seq, 
                          edit_from='C', edit_to='T', window=(4, 8)):
    """碱基编辑效率分析"""

    # 确定编辑窗口位置
    guide_start = reference.find(guide_seq)
    if guide_start == -1:
        guide_start = reference.find(str(Seq(guide_seq).reverse_complement()))

    edit_positions = list(range(guide_start + window[0] - 1, guide_start + window[1]))

    # 统计各位置的编辑
    total_reads = 0
    edited_reads = 0
    position_edits = {pos: 0 for pos in edit_positions}

    for record in SeqIO.parse(fastq_file, "fastq"):
        seq = str(record.seq)
        # 简化: 假设已与reference对齐
        alignment = pairwise2.align.globalms(reference, seq, 2, -1, -5, -1, 
                                              one_alignment_only=True)
        if not alignment:
            continue

        total_reads += 1
        read_edited = False

        ref_aligned, query_aligned = alignment[0][0], alignment[0][1]
        ref_pos = 0
        for i, (r, q) in enumerate(zip(ref_aligned, query_aligned)):
            if r != '-':
                if ref_pos in edit_positions:
                    if r == edit_from and q == edit_to:
                        position_edits[ref_pos] += 1
                        read_edited = True
                ref_pos += 1

        if read_edited:
            edited_reads += 1

    # 计算效率
    results = {
        'total_reads': total_reads,
        'edited_reads': edited_reads,
        'overall_efficiency': edited_reads / total_reads if total_reads > 0 else 0,
        'position_efficiency': {pos: count/total_reads for pos, count in position_edits.items()}
    }

    return results

# 使用示例
results = analyze_base_editing(
    "sample_merged.fq",
    reference="ATCGATCGATCGATCGATCGATCGATCG",
    guide_seq="ATCGATCGATCGATCGATCG",
    edit_from='C', edit_to='T',
    window=(4, 8)
)

print(f"总reads: {results['total_reads']}")
print(f"编辑效率: {results['overall_efficiency']:.1%}")
print("各位置效率:")
for pos, eff in results['position_efficiency'].items():
    print(f"  Position {pos}: {eff:.1%}")

# === 编辑产物分类 ===
# 碱基编辑产物类型:
# 1. 完美编辑(desired edit only)
# 2. Bystander编辑(编辑窗口内其他碱基也被编辑)
# 3. Indel(切口导致的插入缺失，碱基编辑的副作用)
# 4. 未编辑

def classify_editing_outcomes(crispresso_output_dir):
    """分类CRISPResso2输出的编辑结果"""

    # 读取CRISPResso2的等位基因频率表
    alleles = pd.read_csv(
        f"{crispresso_output_dir}/Alleles_frequency_table.zip",
        sep='\t', compression='zip'
    )

    outcomes = {
        'perfect_edit': 0,
        'bystander_edit': 0,
        'indel': 0,
        'unedited': 0
    }

    target_position = 5  # 目标编辑位点(相对于guide起始)

    for _, row in alleles.iterrows():
        aligned_seq = row['Aligned_Sequence']
        ref_seq = row['Reference_Sequence']
        freq = row['%Reads'] / 100

        has_indel = '-' in aligned_seq or '-' in ref_seq

        if has_indel:
            outcomes['indel'] += freq
        else:
            # 检查目标位点和bystander
            edits = []
            for i, (r, a) in enumerate(zip(ref_seq, aligned_seq)):
                if r != a and r != '-' and a != '-':
                    edits.append(i)

            if len(edits) == 0:
                outcomes['unedited'] += freq
            elif len(edits) == 1 and target_position in edits:
                outcomes['perfect_edit'] += freq
            else:
                outcomes['bystander_edit'] += freq

    return outcomes

# 可视化编辑产物
import matplotlib.pyplot as plt

outcomes = classify_editing_outcomes("crispresso_output/")
labels = list(outcomes.keys())
sizes = list(outcomes.values())

fig, ax = plt.subplots(figsize=(8, 6))
colors = ['#2ecc71', '#f39c12', '#e74c3c', '#95a5a6']
ax.pie(sizes, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90)
ax.set_title("Base Editing Outcome Classification")
plt.savefig("editing_outcomes.png", dpi=150)

# === Cas-OFFinder: 预测脱靶位点 ===
# 下载: http://www.rgenome.net/cas-offinder/
# 输入文件格式:
# /path/to/genome_dir
# NNNNNNNNNNNNNNNNNNNNNGG    (PAM pattern)
# ATCGATCGATCGATCGATCGNNN 3  (guide + 允许的错配数)

cat > input.txt << 'EOF'
/path/to/hg38/
NNNNNNNNNNNNNNNNNNNNNGG
ATCGATCGATCGATCGATCGNNN 4
EOF

cas-offinder input.txt C output.txt
# 输出: 预测的脱靶位点列表

# === 脱靶位点的深度测序验证 ===
# 对top 10预测脱靶位点设计引物
# PCR扩增 → 深度测序 → CRISPResso2分析

# 批量分析多个脱靶位点
for site in OT1 OT2 OT3 OT4 OT5; do
  CRISPResso --fastq_r1 ${site}_R1.fq.gz \
    --fastq_r2 ${site}_R2.fq.gz \
    --amplicon_seq ${site}_ref \
    --guide_seq ${site}_guide \
    --base_editor_output \
    --conversion_nuc_from C \
    --conversion_nuc_to T \
    --output_folder offtarget_results/${site}/ \
    --name ${site}
done

import matplotlib.pyplot as plt
import numpy as np

# === 编辑窗口热图 ===
# 多个sgRNA在不同位置的编辑效率
guides = ['sgRNA_1', 'sgRNA_2', 'sgRNA_3', 'sgRNA_4']
positions = list(range(1, 21))  # protospacer 20个位置

# 模拟数据(实际从CRISPResso2结果提取)
efficiency_matrix = np.random.rand(len(guides), 20) * 0.5
# 编辑窗口(位置4-8)效率更高
efficiency_matrix[:, 3:8] *= 3

fig, ax = plt.subplots(figsize=(12, 4))
im = ax.imshow(efficiency_matrix, cmap='YlOrRd', aspect='auto', vmin=0, vmax=1)
ax.set_xticks(range(20))
ax.set_xticklabels(positions)
ax.set_yticks(range(len(guides)))
ax.set_yticklabels(guides)
ax.set_xlabel("Position in protospacer")
ax.set_ylabel("sgRNA")
plt.colorbar(im, label="Editing efficiency")
ax.set_title("Base Editing Window Profile")

# 标注编辑窗口
ax.axvline(x=2.5, color='blue', linestyle='--', alpha=0.5)
ax.axvline(x=7.5, color='blue', linestyle='--', alpha=0.5)
ax.text(5, -0.7, 'Editing Window', ha='center', color='blue')

plt.tight_layout()
plt.savefig("editing_window.png", dpi=150)

# === 完整碱基编辑分析流程 ===

# 1. 质控
fastp -i R1.fq.gz -I R2.fq.gz -o clean_R1.fq.gz -O clean_R2.fq.gz

# 2. CRISPResso2分析
CRISPResso --fastq_r1 clean_R1.fq.gz --fastq_r2 clean_R2.fq.gz \
  --amplicon_seq AMPLICON_SEQUENCE \
  --guide_seq GUIDE_20NT \
  --base_editor_output \
  --conversion_nuc_from C --conversion_nuc_to T \
  --output_folder results/ --name experiment1 \
  --min_average_read_quality 30 \
  --min_single_bp_quality 20

# 3. 批量处理
CRISPRessoBatch --batch_settings samples.txt \
  --base_editor_output --conversion_nuc_from C --conversion_nuc_to T \
  --output_folder batch_results/

# 4. 可视化结果
# CRISPResso2自动生成HTML报告和PNG图片
# 关键图: nucleotide_frequency, base_editor_substitution