环状RNA检测CIRI¶

circRNA的产生机制——反向剪接（back-splicing）：

正常线性剪接：  外显子1 → 外显子2 → 外显子3（5'→3'方向）
反向剪接：      外显子2的3'端 ← 连接 → 外显子2的5'端（形成环）

               ┌──────────────┐
               │   circRNA    │
               │  外显子2-3   │
               └──────────────┘

BSJ（Back-Spliced Junction）就是这个连接点
测序reads跨越BSJ就是circRNA的证据

circRNA的特点：
1. 共价闭合环 → 没有5'帽和3'poly(A)尾
2. 抗RNase R降解 → 比线性RNA更稳定
3. 组织/发育阶段特异性表达
4. 可作为miRNA海绵（sponge）吸附miRNA
5. 少数可以翻译成蛋白质

从RNA-seq中检测circRNA的核心思路：
  找到跨越BSJ的reads → 这些reads不能映射到正常基因组上

两类方法：
  ① 基于比对的方法（alignment-based）：
     BWA-MEM → CIRI2/CIRI3
     STAR → CIRCexplorer2
     Segemehl → 自带嵌合比对

  ② 基于伪参考的方法：
     构建所有可能的BSJ参考序列 → 直接比对

挑战：
  - BSJ reads数量少（通常只有几条）
  - 线性RNA的假嵌合比对会产生假阳性
  - 重复序列区域比对错误
  - 需要高覆盖深度

# ===== 安装CIRI2 =====
# CIRI2基于Perl，需要BWA和SAMtools
conda install -c bioconda bwa         # 安装BWA比对工具
conda install -c bioconda samtools    # 安装SAMtools
conda install -c bioconda ciri2       # 安装CIRI2（如果conda有的话）

# 或手动安装
wget https://sourceforge.net/projects/ciri/files/CIRI2/CIRI_v2.0.6.zip  # 下载
unzip CIRI_v2.0.6.zip       # 解压
chmod +x CIRI2.pl            # 赋予执行权限

# ===== 安装CIRI3（2026新版） =====
# CIRI3基于Java，一站式整合工具
# 从GitHub下载最新版
wget https://github.com/bioinfo-biols/CIRI3/releases/latest/download/CIRI3.jar  # 下载
java -jar CIRI3.jar --help   # 测试安装

# ===== 安装CIRIquant =====
pip install CIRIquant          # 安装定量工具

# 1. 建立BWA索引（CIRI2需要BWA-MEM比对）
bwa index -a bwtsw hg38.fa   # 建立BWA索引（-a bwtsw适合大基因组）

# 2. BWA-MEM比对（必须用BWA-MEM，不能用bowtie2）
bwa mem -T 19 \               # 最小分数阈值19（允许嵌合比对）
  -t 8 \                      # 8个线程
  hg38.fa \                   # 参考基因组
  reads_R1.fq.gz \            # 正向reads
  reads_R2.fq.gz \            # 反向reads
  > aligned.sam                # 输出SAM文件（不要转BAM！CIRI2需要SAM）

# 注意：CIRI2需要SAM格式输入，不要排序或转换为BAM
# 这是和其他工具不同的地方

# 3. 准备基因注释文件（GTF格式）
# 从GENCODE下载
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.annotation.gtf.gz  # 下载注释
gunzip gencode.v46.annotation.gtf.gz  # 解压

# 运行CIRI2检测circRNA
perl CIRI2.pl \
  -I aligned.sam \            # BWA-MEM比对的SAM文件
  -O circRNA_output.txt \     # 输出文件
  -F hg38.fa \                # 参考基因组FASTA
  -A gencode.v46.annotation.gtf \  # 基因注释GTF
  -T 8 \                      # 线程数
  -G circRNA.log              # 日志文件

# CIRI2输出列说明：
# circRNA_ID:    circRNA唯一标识（chr:start|end）
# chr:           染色体
# circRNA_start: 起始位置
# circRNA_end:   终止位置
# #junction_reads: BSJ支持reads数（越多越可靠）
# SM_MS_SMS:     BSJ reads的比对模式分类
# #non_junction_reads: 非BSJ reads数
# junction_reads_ratio: BSJ reads占比
# circRNA_type:  类型（exon/intron/intergenic）
# gene_id:       宿主基因ID
# strand:        链方向
# junction_reads_ID: BSJ reads的ID列表

# CIRI3整合了检测+定量+差异分析
# 步骤1：检测circRNA
java -jar CIRI3.jar detect \
  -i aligned.sam \             # SAM输入
  -o sample_circRNA \          # 输出前缀
  -r hg38.fa \                 # 参考基因组
  -a gencode.v46.annotation.gtf \  # 注释文件
  -t 8                         # 线程数

# 步骤2：精确定量（类似CIRIquant）
java -jar CIRI3.jar quant \
  -i sample_circRNA.txt \      # 步骤1的输出
  -r hg38.fa \                 # 参考基因组
  -b aligned.bam \             # BAM文件
  -o sample_quant               # 输出前缀

# 步骤3：差异分析（多样本比较）
java -jar CIRI3.jar de \
  -i sample_list.txt \         # 样本列表文件
  -o diff_circRNA \            # 输出前缀
  -g group_info.txt            # 分组信息

# CIRI3相比CIRI2的改进：
# - 更高的灵敏度和准确性
# - 整合检测+定量+差异分析
# - 支持更多物种
# - 改进的假阳性过滤

# 安装CIRCexplorer2
pip install circexplorer2       # pip安装

# CIRCexplorer2有5个模块：
# align   — 比对模块
# parse   — 解析嵌合比对
# annotate — 注释circRNA
# assemble — 组装circRNA全长序列
# denovo   — 从头检测未注释的circRNA

# ===== 方法1：使用STAR比对 =====
# 先用STAR比对（需要支持嵌合比对的参数）
STAR --runThreadN 8 \                           # 线程数
  --genomeDir star_index/ \                      # STAR索引目录
  --readFilesIn reads_R1.fq.gz reads_R2.fq.gz \ # 输入文件
  --readFilesCommand zcat \                      # gzip文件用zcat
  --chimSegmentMin 10 \                          # 嵌合比对最小片段长度
  --chimOutType Junctions \                      # 输出嵌合连接信息
  --outSAMtype BAM SortedByCoordinate \          # 输出排序BAM
  --outFileNamePrefix sample_                    # 输出前缀

# 解析STAR的嵌合比对结果
CIRCexplorer2 parse \
  -t STAR \                    # 比对工具类型
  sample_Chimeric.out.junction \  # STAR的嵌合连接文件
  > CIRCexplorer2_parse.log    # 日志

# 注释circRNA
CIRCexplorer2 annotate \
  -r ref_genePred.txt \        # 基因注释（genePred格式）
  -g hg38.fa \                 # 参考基因组
  -b back_spliced_junction.bed \  # parse步骤的输出
  -o circRNA_known.txt         # 输出已知circRNA

# genePred格式转换（从GTF转）
# gtfToGenePred gencode.v46.annotation.gtf ref_genePred.txt

# 从头检测（找未注释的新circRNA）
CIRCexplorer2 denovo \
  -r ref_genePred.txt \        # 基因注释
  -g hg38.fa \                 # 参考基因组
  -b back_spliced_junction.bed \  # BSJ文件
  --abs abs_folder/ \          # TopHat比对结果目录
  -o circRNA_denovo.txt        # 输出文件

# 输出文件列说明：
# chrom: 染色体
# start: 起始位置
# end: 终止位置
# name: circRNA名称
# score: BSJ reads数
# strand: 链方向
# thickStart/thickEnd: 编码区范围
# itemRgb: 颜色
# exonCount: 外显子数
# exonSizes: 各外显子大小
# exonOffsets: 各外显子偏移
# readNumber: BSJ reads数
# circType: circRNA类型
# geneName: 基因名
# isoformName: 转录本名
# flankIntron: 侧翼内含子

# CIRIquant对circRNA做精确定量
# 原理：构建伪参考序列 → 重新比对 → 消除线性RNA干扰

CIRIquant \
  -t 8 \                        # 线程数
  -1 reads_R1.fq.gz \           # 正向reads
  -2 reads_R2.fq.gz \           # 反向reads
  --config config.yml \          # 配置文件（指定工具路径和参考基因组）
  -o output_dir/ \               # 输出目录
  --circ circRNA_output.txt \    # CIRI2的检测结果
  -p sample_prefix               # 输出前缀

# 配置文件config.yml内容：
# name: hg38
# tools:
#   bwa: /path/to/bwa
#   hisat2: /path/to/hisat2
#   stringtie: /path/to/stringtie
#   samtools: /path/to/samtools
# reference:
#   fasta: /path/to/hg38.fa
#   gtf: /path/to/gencode.v46.annotation.gtf
#   bwa_index: /path/to/bwa_index/hg38

# CIRIquant输出列：
# circ_id: circRNA标识
# circ_type: 类型
# bsj: BSJ reads数
# ratio: circRNA/(circRNA+线性RNA)的比值
# exp: 表达量（CPM）

# 读取CIRI2输出
library(tidyverse)

circ_data <- read.delim("circRNA_output.txt")  # 读取CIRI2结果
cat("检测到circRNA数量:", nrow(circ_data), "\n")  # 总数

# 过滤低质量circRNA
circ_filtered <- circ_data %>%
  filter(`X.junction_reads` >= 2) %>%   # 至少2条BSJ reads
  filter(circRNA_type == "exon")        # 只保留外显子型circRNA
cat("过滤后circRNA数量:", nrow(circ_filtered), "\n")

# circRNA类型分布
type_counts <- table(circ_data$circRNA_type)  # 统计各类型数量
barplot(type_counts,
        main = "circRNA Type Distribution",  # 标题
        xlab = "Type",                        # X轴标签
        ylab = "Count",                       # Y轴标签
        col = "steelblue")                    # 颜色

# circRNA长度分布
circ_data$length <- circ_data$circRNA_end - circ_data$circRNA_start  # 计算长度
ggplot(circ_data, aes(x = log10(length))) +
  geom_histogram(binwidth = 0.1, fill = "steelblue") +  # 直方图
  labs(x = "log10(circRNA Length)", y = "Count",
       title = "circRNA Length Distribution") +
  theme_minimal()

# 宿主基因分析
host_genes <- circ_data %>%
  group_by(gene_id) %>%                    # 按宿主基因分组
  summarise(n_circ = n()) %>%              # 统计每个基因产生的circRNA数
  arrange(desc(n_circ))                    # 按数量排序
head(host_genes, 10)  # 查看产生最多circRNA的基因

# circRNA与miRNA相互作用预测（需要circBank或circInteractome数据库）
# 大部分circRNA通过miRNA海绵机制发挥功能

# 推荐使用CirComPara2整合多个工具的结果
# CirComPara2整合了CIRI2、CIRCexplorer2、find_circ、DCC、Segemehl

# 安装CirComPara2
git clone https://github.com/egaffo/CirComPara2.git  # 克隆
cd CirComPara2
pip install -r requirements.txt  # 安装依赖

# 运行CirComPara2
# 配置文件指定样本信息和工具路径
python CirComPara2.py \
  --config config.yaml \       # 配置文件
  --threads 8                  # 线程数

# CirComPara2的策略：
# 1. 同时运行多个检测工具
# 2. 只保留至少2个工具都检测到的circRNA
# 3. 大大降低假阳性率

# 2024年大规模Benchmark结果：
# - 16个工具检测>315,000个circRNA
# - 工具间精确度相似（中位数>95%）
# - 灵敏度差异巨大（1,372到58,032个）
# - 推荐使用多工具取交集

# circRNA检测流程
BWA-MEM比对(SAM) → CIRI2/CIRI3检测BSJ
STAR比对(嵌合) → CIRCexplorer2解析+注释
多工具整合 → CirComPara2取交集（推荐）

# 工具选择
快速检测: CIRI2（BWA-MEM → SAM → CIRI2）
全面分析: CIRCexplorer2（5个模块）
精确定量: CIRIquant（伪参考序列法）
一站式:   CIRI3（2026年，检测+定量+差异）
多工具:   CirComPara2（≥2工具交集）

# 过滤标准
BSJ reads: ≥2条（严格≥5条）
circRNA类型: exon型最可靠
RNase R验证: 实验验证金标准
多工具验证: ≥2个工具检测到

# circRNA命名规则
格式: chr:start|end
示例: chr1:12345|67890
含义: 染色体:BSJ起始位置|BSJ终止位置

# 关键参数
BWA-MEM -T: 19（允许短片段嵌合比对）
STAR chimSegmentMin: 10（嵌合片段最小长度）
CIRI2线程: -T 8（根据CPU调整）
CIRCexplorer2: 需要genePred格式注释