CLIP-seq RNA结合蛋白分析¶

# CLIP数据特点
# - 单端测序为主（eCLIP有PE版本）
# - reads长度短且不均一
# - PCR重复率高（需UMI去重）
# - 需要size-matched input作为背景控制

# 数据文件
# IP样本：RBP-bound RNA fragments
# Input/SMInput：size-matched input control

# === eCLIP数据预处理（ENCODE pipeline）===

# 1. 接头去除（cutadapt，eCLIP有两轮接头）
cutadapt -a AGATCGGAAGAGC -O 1 --times 2 -m 18 \
    -o trimmed_R1.fq.gz raw_R1.fq.gz
# 如果有UMI在R2中：
cutadapt -a AGATCGGAAGAGC -O 1 --times 2 -m 18 \
    -o trimmed_R2.fq.gz raw_R2.fq.gz

# 2. UMI提取（如果有）
umi_tools extract --stdin=trimmed_R1.fq.gz --stdout=umi_R1.fq.gz \
    --bc-pattern=NNNNNNNNNN  # 10bp UMI

# 3. 比对（STAR）
STAR --runThreadN 16 \
    --genomeDir star_index \
    --readFilesIn umi_R1.fq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outFilterMultimapNmax 1 \
    --outFilterMismatchNmax 2 \
    --outFileNamePrefix eclip_

# 4. UMI去重
umi_tools dedup -I eclip_Aligned.sortedByCoord.out.bam \
    -S eclip_dedup.bam --method unique

# 或使用坐标去重（无UMI时）
samtools markdup -r eclip_sorted.bam eclip_dedup.bam

# 5. 同样处理SMInput
# ... 相同步骤处理 input 样本 ...

# === PureCLIP（基于HMM的精确交联位点检测）===
pureclip -i eclip_dedup.bam -bai eclip_dedup.bam.bai \
    -g genome.fa \
    -o pureclip_crosslink_sites.bed \
    -or pureclip_regions.bed \
    -nt 16 \
    -iv 'chr1;chr2;chr3' \  # 用于训练参数的染色体
    -dm 8  # merge distance

# === CLIPper（ENCODE eCLIP官方peak caller）===
clipper -b eclip_dedup.bam \
    -s hg38 \
    -o clipper_peaks.bed \
    --processors 16

# 用input做归一化过滤
# 计算IP/Input的fold enrichment
# 保留 log2(IP/Input) > 3 的peaks

# === Piranha ===
Piranha -s eclip_dedup.bam \
    -o piranha_peaks.bed \
    -p 0.01 \
    -b 50  # bin size

# === 使用IDR合并重复（ENCODE标准）===
idr --samples rep1_peaks.bed rep2_peaks.bed \
    --input-file-type bed \
    --output-file idr_peaks.bed \
    --idr-threshold 0.01

# === CTK iCLIP分析 ===
# CTK: https://github.com/chaolinzhanglab/ctk

# 1. 解析和去重
perl ctk/stripBarcode.pl -format fastq \
    -len 5 raw.fq.gz stripped.fq.gz barcode.txt

# 比对后去PCR重复（基于barcode+mapping position）
perl ctk/tag2collapse.pl -v --random-barcode \
    -EM 30 --seq-error-model em-local \
    eclip_sorted.bed eclip_uniq.bed

# 2. 识别交联位点（CITS: CrossLink-Induced Truncation Sites）
perl ctk/tag2peak.pl -big -ss -v \
    --valley-seeking --valley-depth 0.9 \
    -p 0.001 \
    eclip_uniq.bed eclip_peaks.bed \
    --out-boundary eclip_boundaries.bed

# 3. 获取单碱基分辨率的交联位点
perl ctk/CITS.pl -big -ss -v \
    -p 0.001 \
    eclip_uniq.bed eclip_CITS.bed

# === Motif发现 ===

# 提取peak区域序列
bedtools getfasta -fi genome.fa -bed peaks.bed -fo peak_sequences.fa -s

# HOMER motif发现
findMotifsGenome.pl peaks.bed hg38 homer_output/ \
    -rna -size 50 -mis 1 -p 8

# MEME/DREME
dreme -rna -p peak_sequences.fa -o dreme_output/
meme peak_sequences.fa -rna -oc meme_output/ -mod zoops -nmotifs 10 -minw 4 -maxw 10

# === 功能注释 ===
library(GenomicRanges)
library(GenomicFeatures)

# 加载基因注释
txdb <- makeTxDbFromGFF("gencode.v38.gtf")
peaks_gr <- import("peaks.bed")

# 注释peak位置
utr3 <- threeUTRsByTranscript(txdb, use.names = TRUE)
utr5 <- fiveUTRsByTranscript(txdb, use.names = TRUE)
cds <- cdsBy(txdb, by = "tx", use.names = TRUE)
introns <- intronsByTranscript(txdb, use.names = TRUE)

# 统计分布
n_3utr <- sum(countOverlaps(peaks_gr, unlist(utr3)) > 0)
n_5utr <- sum(countOverlaps(peaks_gr, unlist(utr5)) > 0)
n_cds <- sum(countOverlaps(peaks_gr, unlist(cds)) > 0)
n_intron <- sum(countOverlaps(peaks_gr, unlist(introns)) > 0)

cat(sprintf("3'UTR: %d (%.1f%%)\n5'UTR: %d (%.1f%%)\nCDS: %d (%.1f%%)\nIntron: %d (%.1f%%)\n",
    n_3utr, n_3utr/length(peaks_gr)*100,
    n_5utr, n_5utr/length(peaks_gr)*100,
    n_cds, n_cds/length(peaks_gr)*100,
    n_intron, n_intron/length(peaks_gr)*100))

# === CLIP + RNA-seq整合 ===

# 1. 确定RBP的直接靶基因
# 有peak的基因 = 直接靶标
target_genes <- unique(peaks_annotated$gene_name)

# 2. 与RBP KD后的差异基因比较
kd_de_genes <- read.csv("RBP_knockdown_DEGs.csv")

# 直接靶标中有多少差异表达
direct_targets_de <- intersect(target_genes, kd_de_genes$gene[kd_de_genes$padj < 0.05])
cat(sprintf("Direct targets with expression change: %d / %d (%.1f%%)\n",
    length(direct_targets_de), length(target_genes),
    length(direct_targets_de)/length(target_genes)*100))

# 3. Fisher精确检验：RBP靶标是否富集于差异基因
all_genes <- unique(kd_de_genes$gene)
fisher_table <- matrix(c(
  length(intersect(target_genes, kd_de_genes$gene[kd_de_genes$padj < 0.05])),
  length(setdiff(target_genes, kd_de_genes$gene[kd_de_genes$padj < 0.05])),
  sum(kd_de_genes$padj < 0.05) - length(direct_targets_de),
  length(all_genes) - length(target_genes) - sum(kd_de_genes$padj < 0.05) + length(direct_targets_de)
), nrow = 2)
fisher.test(fisher_table)

# 4. 调控方向：3'UTR结合通常促降解
# 靶标在KD后上调 → RBP促进降解
# 靶标在KD后下调 → RBP促进稳定/翻译

# eCLIP标准流程
cutadapt -a ADAPTER -m 18 -o trimmed.fq.gz raw.fq.gz
STAR --genomeDir index --readFilesIn trimmed.fq.gz --outSAMtype BAM SortedByCoordinate
umi_tools dedup -I sorted.bam -S dedup.bam
clipper -b dedup.bam -s hg38 -o peaks.bed
findMotifsGenome.pl peaks.bed hg38 motif_out/ -rna