CircRNA-seq上游分析工具测评：CIRIquant VS. CIRCexplorer3
前言
本次测评CircRNA-seq上游分析的两大最新工具CIRCexplorer3及CIRIquant。CIRCexplorer3是2019年发表在Genomics Proteomics Bioinformatics（2020 IF=7.69）上，目前引用量是22次；CIRIquant2020年发表在nature communications上，目前引用量是54次。
不考虑算法的前提下比较这两款软件：两款软件运行均比较慢，40个线程下双端测序的一个样本约需2小时。其中CIRCexplorer3运行更慢一些，且需要安装非常多的依赖包。但是，CIRCexplorer3github官网并没有列出这些依赖包，且git clone安装后，也没有提供依赖包的list文件；尽管这款软件可用一句代码即可完成所有模块，且得到最终计算出的表达矩阵，但是在缺少部分依赖包的情况下，这款软件把能运行的部分运行完毕，但一遇到当前模块缺少对应依赖包时，便报错停止运行，且报错信息给的非常不明确。这样的话，造成了软件给出什么报错，我们苦苦debug，对应的装当前模块缺少的依赖包。但由于软件本身运行非常耗时，哪怕漏一个依赖包，都得从头再来......最后需要注意的是，这款软件对依赖包的版本有限制（笔者血泪教训得到的经验），这一点github官网也没有明确说清楚，我是看了issue部分才知道的。
ps：笔者前期并未试用过
CIRCexplorer1和CIRCexplorer2，笔者推测使用过这两个旧版本的用户更适合上手CIRCexplorer3。
而相比之下，CIRIquant非常贴心，不仅提供一键式下载，还有示例数据及官方的文档，见https://ciri-cookbook.readthedocs.io/en/latest/index.html。
因此，笔者推荐Linux熟练度不够的同学，不建议尝试CIRCexplorer3。另一个角度考虑的话，CIRIquant发表在NC上，且是最新的软件。当然，想要折腾一下，尝试不同算法的同学，两个软件均可做尝试。
接下来，是两款软件的测评情况。我们先从示例数据的下载开始。
一. 示例数据下载
此模块对于有过RNA-seq经验的同学来说可以略过。
Step0. 数据下载前的准备：
软件下载使用conda和mamba：
注：mamba是conda的一个超级加速软件
conda create -n circrna python=2
mamba install -y -c bioconda sra-tools #GEO官方推荐下载工具
mamba install -y -c bioconda fastq-dump #SRR转fq工具
mamba install -y -c bioconda fastp #质控与修剪
mamba install -y -c bioconda fastqc #fq文件质控
mamba install -y -c bioconda trim-galore #fq文件修剪
Step1. 首先是获取SRA ID：
笔者使用的数据为项目数据，这里我们以GSE108505为例：
点击上图的SRA Run Selector按钮，到达下图界面，点击下载Metadata即为表型信息，点击Accession List即可获得SRA ID。
根据表型数据可知，9个数据为双端数据，且9个样分为3组处理。开始下载数据：
mkdir 1.sra_data
cd 1.sra_data
###1.Raw_data: 下载SRR数据
cat >cirRNA.id #注：随后键入Enter键，输入如下内容：
SRR6417989
SRR6417990
SRR6417991
SRR6417992
SRR6417993
SRR6417994
SRR6417995
SRR6417996
SRR6417997
#注：快捷键ctrl+C即可结束当前的输入
##批量下载
cat cirRNA.id |while read id;do (nohup prefetch -c -p $id 1> $id.log &);done # 后台下载
##把所有下载好的SRR移动到当前目录
mv SRR*/* ./
rm -r SRR*/
##检查文件下载的大小和完整性
##cat nohup.out | grep "failed"
Step2. SRR转fastq文件
###2.SRR转fastq
mkdir cd ../2.raw_fastq/
cd ../2.raw_fastq/
##做一个软连接文件
ln -s ../1.sra_data/*.sra ../2.raw_fastq/
##创建文件转换fastq.sh脚本文件，注意示例数据为双端数据
cat >paired_fastq.sh #注：随后键入Enter键，输入如下内容：
ls *.sra | while read id;do ( nohup fastq-dump --gzip --split-3 -O ./ $id & );done
#注：快捷键ctrl+C即可结束当前的输入
##运行，即可将单个SRR文件转为_1.fastq.gz和_2.fastq.gz
bash paired_fastq.sh
Step3. fastp质控
###3.fastq_results
mkdir ../3.fastq_results/
cd ../3.fastq_results/
ln -s ../2.raw_fastq/*.fastq.gz ./
##生成的fastqc放到3.fastq_results中,-t指定30个线程数。
nohup fastqc -t 30 -o ./ ./*.fastq.gz &
##使用MutliQC整合FastQC结果，将后缀为.html的文件进行multiqc处理
multiqc ./
Step4. clean: trim_galore进行过滤
原始数据已去接头。这里贴一下用trim_galore去接头的代码：
mkdir ../4.trim_galore_clean
cd ../4.trim_galore_clean
##如果9个样本均需要过滤
cat >id.txt #注：随后键入Enter键，输入如下内容：
SRR6417989
SRR6417990
SRR6417991
SRR6417992
SRR6417993
SRR6417994
SRR6417995
SRR6417996
SRR6417997
#注：快捷键ctrl+C即可结束当前的输入
##写入bash命令
cat >trim_galore.bash #注：随后键入Enter键，输入如下内容：
cat id.txt | while read id;
do
trim_galore --quality 25 --phred33
--length 36 -j 20 -e 0.1 --stringency 3
-o ./
--paired ${id}_1.fastq.gz ${id}_2.fastq.gz
> ./${id}_trim.log 2>&1;
done
#注：快捷键ctrl+C即可结束当前的输入
## 开始过滤
nohup bash trim_galore.bash &
拿到clean后的fastq文件即可使用CIRIquant或CIRCexplorer3软件输出一系列的中间文件包括bam数据 ，及最终的表达数据
二. CIRIquant分析circRNA-seq fastq原始数据
原文是Zhang, J., Chen, S., Yang, J. et al. Accurate quantification of circular RNAs identifies extensive circular isoform switching events. Nat Commun 11, 90 (2020) doi:10.1038/s41467-019-13840-9，
官网教程：https://ciri-cookbook.readthedocs.io/en/latest/index.html
Step1. 软件准备
官网列出的主要软件包括：
- Python 2.7
- bwa 0.7.17-r1188
- hisat2 v2.1.0
- stringtie v2.0
- samtools >=v1.9
- Perl v5.26.2
###1.常规的软件下载：
conda create -n ciriquant python=2
conda activate ciriquant
conda install -y bwa stringtie
conda install -y hisat2=2.1.0
conda install -y samtools=1.9
conda install -c conda-forge numpy
pip install --upgrade pip
###2.安装CIRIquant
pip install CIRIquant
Step2. 先对一个样本进行定量
mkdir 5.final_matrix
###1.根据官网说明，需要制作一个chr1.yml文件
##注:以下软件的路径及参考基因组路径需自定义
cat >chr1.yml #注：随后键入Enter键，输入如下内容：
name: mm10
tools:
# 软件路径
bwa: /home/data/vip10t17/miniconda3/envs/ciriquant/bin/bwa
hisat2: /home/data/vip10t17/miniconda3/envs/ciriquant/bin/hisat2
stringtie: /home/data/vip10t17/miniconda3/envs/ciriquant/bin/stringtie
samtools: /home/data/vip10t17/miniconda3/envs/ciriquant/bin/samtools
reference:
# 文件路径
fasta: /home/data/vip10t17/mouse_rnaseq/3.cirRNA/ref_mouse/mm10.fa
gtf: /home/data/vip10t17/mouse_rnaseq/3.cirRNA/ref_mouse/mm10_ref.gtf
# 索引路径
bwa_index: /home/data/refdir/index/bwa/mm10
hisat_index: /home/data/refdir/index/hisat/mm10/genome
#注：快捷键ctrl+C即可结束当前的输入
###2.先试一个样本
mkdir test_out_dir
CIRIquant -t 40
-1 ./SRR6417989_1.fastq.gz
-2 ./SRR6417989_2.fastq.gz
--config ./chr1.yml
-o test_out_dir
-p SRR6417989
**注意：**一个样本最终的占用空间大约为11G，但是中间会生成临时的sam文件，约占90G空间（结束后会自动删除），确保硬盘空间大小足够
40个线程的运行量下，一个样本大概需要耗时1个小时40分钟，输出文件如下：
cd ./test_out_dir
tree -h
查看最终的输出结果：
cat SRR6417989.gtf | head
想要看懂这些结果，需要浏览官网：https://ciri-cookbook.readthedocs.io/en/latest/CIRIquant_2_quantification.html
大致如下：
Output format
The main output of CIRIquant is a GTF file, that contains detailed information of BSJ and FSJ reads of circRNAs and annotation of circRNA back-spliced regions in the attribute columns
Description of each columns’s value
column | name | description |
|---|---|---|
1 | chrom | chromosome / contig name |
2 | source | CIRIquant |
3 | type | circRNA |
4 | start | 5' back-spliced junction site |
5 | end | 3' back-spliced junction site |
6 | score | CPM of circRNAs (#BSJ / #Mapped reads) |
7 | strand | strand information |
8 | . | . |
9 | attributes | attributes seperated by semicolon |
The attributes containing several pre-defined keys and values:
key | description |
|---|---|
circ_id | name of circRNA |
circ_type | circRNA types: exon / intron / intergenic |
bsj | number of bsj reads |
fsj | number of fsj reads |
junc_ratio | circular to linear ratio: 2 * bsj / ( 2 * bsj + fsj) |
rnaser_bsj | number of bsj reads in RNase R data (only when --RNaseR is specificed) |
rnaser_fsj | number of fsj reads in RNase R data (only when --RNaseR is specificed) |
gene_id | ensemble id of host gene |
gene_name | HGNC symbol of host gene |
gene_type | type of host gene in gtf file |
Step3. 然后进行批量处理
### 批量处理
for i in `ls *_1.fastq.gz`
do
id=${i/_1.fastq.gz/}
mkdir ${id}_out
echo "CIRIquant -t 40
-1 ./${id}_1.fastq.gz
-2 ./${id}_2.fastq.gz
--config ./chr1.yml
-o ${id}_out
-p ${id}"
done >ciriquant.sh
nohup bash ciriquant.sh 1>nohup.log &
将输出的gtf文件拷贝到集中到gtf_file文件夹下：
mkdir gtf_file
## 写脚本
ls *_1.fastq.gz | while read id
do
id=${id/_1.fastq.gz/}
file=${id}_out
echo "cp $file/${id}.gtf ./gtf_file/"
done >gather.command
## 检查一下
cat gather.command
## 运行
bash gather.command
Step4. 差异分析
该软件还提供了差异分析，总共三步走。
方便起见，这里仅简单做两组差异分析：
- differentiation 3 days (D3) vs. differentiation 6 days (D6)
- differentiation 3 days (D3) vs. proliferation （P）
（1）先准备一些配置文件：
cd ./gtf_file
## group ("C" for control, "T" for treatment)
##D3 VS. D6
cat >sample_1.lst #注：随后键入Enter键，输入如下内容
D3_1 ./SRR6417989.gtf C 1
D3_2 ./SRR6417990.gtf C 2
D3_3 ./SRR6417991.gtf C 3
D6_1 ./SRR6417992.gtf T 1
D6_2 ./SRR6417993.gtf T 2
D6_3 ./SRR6417994.gtf T 3
#注：快捷键ctrl+C即可结束当前的输入
prep_CIRIquant -i sample_1.lst
--lib library_info_1.csv
--circ circRNA_info_1.csv
--bsj circRNA_bsj_1.csv
--ratio circRNA_ratio_1.csv
#查看生成的文件
ls -lh *_1* | cut -d " " -f 5-
1.6M 11月 28 14:09 circRNA_bsj_1.csv
3.0M 11月 28 14:09 circRNA_info_1.csv
1.7M 11月 28 14:09 circRNA_ratio_1.csv
240 11月 28 14:09 library_info_1.csv
159 11月 28 14:08 sample_1.lst
##D3 VS. P
cat >sample_2.lst #注：随后键入Enter键，输入如下内容
D3_1 ./SRR6417989.gtf C 1
D3_2 ./SRR6417990.gtf C 2
D3_3 ./SRR6417991.gtf C 3
P_1 ./SRR6417995.gtf T 1
P_2 ./SRR6417996.gtf T 2
P_3 ./SRR6417997.gtf T 3
#注：快捷键ctrl+C即可结束当前的输入
#运行
prep_CIRIquant -i sample_2.lst
--lib library_info_2.csv
--circ circRNA_info_2.csv
--bsj circRNA_bsj_2.csv
--ratio circRNA_ratio_2.csv
#查看生成的文件
ls -lh *_2* | cut -d " " -f 5-
3.2M 11月 28 14:11 circRNA_bsj_2.csv
6.2M 11月 28 14:11 circRNA_info_2.csv
3.4M 11月 28 14:11 circRNA_ratio_2.csv
259 11月 28 14:11 library_info_2.csv
168 11月 28 14:08 sample_2.lst
然后可以将这些计数矩阵（csv文件）导入到R中，可以用DESeq2（DESeqDataSetFromMatrix()函数）或edgeR（DGEList()函数）分析。
（2）StringTie的输出文件：
##D3 VS. D6
cat >sample_gene1.lst #注：随后键入Enter键，输入如下内容
D3_1 ../SRR6417989_out/gene/SRR6417989_out.gtf
D3_2 ../SRR6417990_out/gene/SRR6417990_out.gtf
D3_3 ../SRR6417991_out/gene/SRR6417991_out.gtf
D6_1 ../SRR6417992_out/gene/SRR6417992_out.gtf
D6_2 ../SRR6417993_out/gene/SRR6417993_out.gtf
D6_3 ../SRR6417994_out/gene/SRR6417994_out.gtf
#注：快捷键ctrl+C即可结束当前的输入
#产生 gene_count_matrix_1.csv和transcript_count_matrix_1.csv文件
prepDE.py -i sample_gene1.lst -g gene_count_matrix_1.csv -t transcript_count_matrix_1.csv
##D3 VS. P
cat >sample_gene2.lst #注：随后键入Enter键，输入如下内容
D3_1 ../SRR6417989_out/gene/SRR6417989_out.gtf
D3_2 ../SRR6417990_out/gene/SRR6417990_out.gtf
D3_3 ../SRR6417991_out/gene/SRR6417991_out.gtf
P_1 ../SRR6417995_out/gene/SRR6417995_out.gtf
P_2 ../SRR6417996_out/gene/SRR6417996_out.gtf
P_3 ../SRR6417997_out/gene/SRR6417997_out.gtf
#注：快捷键ctrl+C即可结束当前的输入
#产生 gene_count_matrix.csv
prepDE.py -i sample_gene2.lst -g gene_count_matrix_2.csv -t transcript_count_matrix_2.csv
（3）差异分析：
注：需要安装R和edgeR，optparse包
##示例
CIRI_DE_replicate --lib library_info_1.csv #library information by CIRIquant
--bsj circRNA_bsj_1.csv # circRNA expression matrix
--gene gene_count_matrix_1.csv # gene expression matrix
--out circRNA_de_D3_VS_D6.csv #output differential expression result
##PRJNA712946
CIRI_DE_replicate --lib library_info_1.csv --bsj circRNA_bsj_1.csv --gene gene_count_matrix_1.csv --out circRNA_de_D3_VS_D6.csv
##SRP107902
CIRI_DE_replicate --lib library_info_2.csv --bsj circRNA_bsj_2.csv --gene gene_count_matrix_2.csv --out circRNA_de_D3_VS_P.csv
差异分析结果示例：
三. CIRCexplorer3分析circRNA-seq fastq原始数据
Step1. 软件准备
- CIRCexplorer2 (>=2.3.6)，
CIRCexplorer2参考：https://circexplorer2.readthedocs.io/en/latest/ - HISAT2 (>=2.0.5)
- StringTie (>1.3.6)
- Package (python 2.7 +)
- pysam (>=0.8.4)
- pybedtools
- CIRCexplorer3，软件地址：https://github.com/YangLab/CLEAR
conda create -n cicrrna python=2 #此处拼写错误，应该是circrna，软件已经安装好了我就不修改了
##下面的软件带版本号的，严格要求软件的版本
mamba install -y -c bioconda hisat2 pysam stringtie bowtie2 bwa
mamba install -y -c bioconda circexplorer2
conda install -y -c bioconda/label/cf201901 mapsplice
mamba install -y samtools=0.1.18.0
mamba install -y -c bioconda bowtie=1.0.0.0
conda install -y -c bioconda ucsc-genepredtogtf
conda install -y -c bioconda ucsc-gtftogenepred
## 安装最新的CIRCexplorer3
git clone https://github.com/YangLab/CLEAR
cd CLEAR
python ./setup.py install
安装tophat软件:
wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.12.Linux_x86_64.tar.gz
tar -zxvf tophat-2.0.12.Linux_x86_64.tar.gz
cd ./tophat-2.0.12.Linux_x86_64
ln -s /home/data/ssy43/software_install/tophat-2.0.12.Linux_x86_64/tophat ~/miniconda3/envs/cicrrna/bin/
ln -s /home/data/ssy43/software_install/tophat-2.0.12.Linux_x86_64/tophat2 ~/miniconda3/envs/cicrrna/bin/
#还是报错
tophat
# 更改python的版本为2.7
cd ~/software_install/tophat-2.0.12.Linux_x86_64/
vim tophat
vim tophat-fusion-post
Step2. CIRCexplorer3软件所需的基因注释文件及索引构建
CIRCexplorer2 需要基因注释文件和参考基因组序列文件来注释环状 RNA。基因注释文件应为 Gene Predictions 和 RefSeq Genes with Gene Names 格式，参考基因组序列文件包含所有具有各自染色体 ID 的基因组序列。基因注释文件中的所有染色体 ID 都必须包含在参考基因组序列文件中，否则这两个文件之间的不一致可能会导致运行 CIRCexplorer2 时出现不可检测的错误。
作者还提供了一个脚本来下载，可以使用 fetch_ucsc.py 脚本下载所有必需的基因注释和参考基因组序列文件，用于环状 RNA 鉴定。
fetch_ucsc.py 是一个包含在 CIRCexplorer2 中的 Python 小脚本，用于帮助用户为 CIRCexplorer2 准备相关的东西。它可以下载和格式化基因注释文件（RefSeq、KnownGenes 或 Ensembl）和两个物种（人类：hg19、hg38；小鼠：mm9、mm10）的参考基因组序列文件。所有这些文件都将从最新版本的 UCSC 中获取。
如下：
##人类的
#1.下载人类RefSeq数据库基因注释文件
fetch_ucsc.py hg38 ref hg38_ref.txt
#2.下载人类KnownGenes数据库基因注释文件
fetch_ucsc.py hg38 kg hg38_kg.txt
#3.下载人类Ensembl数据库基因注释文件
fetch_ucsc.py hg19 ens hg19_ens.txt
#4.下载人类参考基因组序列文件
fetch_ucsc.py hghg38 fa hg38.fa
##同样小鼠的：
#1.小鼠 RefSeq数据库 基因注释文件
fetch_ucsc.py mm10 ref mm10_ref.txt
#2.小鼠 KnownGenes数据库基因注释文件
fetch_ucsc.py mm10 kg mm10_kg.txt
#3.小鼠 Ensembl数据库基因注释文件
fetch_ucsc.py mm10 ens mm10_ens.txt
#4.小鼠参考基因组序列文件
fetch_ucsc.py mm10 fa mm10.fa
然后转换为 GTF 格式：
# 将基因注释文件转换为GTF格式（需要genePredToGtf）
cut -f2-11 hg38_ref.txt|genePredToGtf file stdin hg38_ref.gtf
cut -f2-11 mm10_ref.txt|genePredToGtf file stdin mm10_ref.gtf
CIRCexploer2 TopHat2/TopHat-Fusion pipeline 需要 Bowtie 和 Bowtie2 索引文件作为参考基因组。可以使用 bowtie-build 和 bowtie2-build 来索引相关基因组。或者可以使用 CIRCexplorer2 align 自动索引基因组文件:
# Bowtie 建索引基因组
bowtie-build --threads 40 mm10.fa bowtie1_index # 输出为.ebwt后缀的索引
# Bowtie2 建索引基因组
# bowtie2-build mm10.fa bowtie2_inde
Step3. 开始分析
mkdir ../6.final_matrix_2/
cd ../6.final_matrix_2/
ln -s ../2.raw_fastq/*.fastq.gz ./
### 1. 先跑一个看看
mm10_hisat_index=/home/data/refdir/server/reference/index/hisat/mm10/genome
mm10_bowtie_index=~/ref_annotation_Geneset/5.circrna/mm10/bowtie1_index
annotation_gtf=/home/data/ssy43/ref_annotation_Geneset/5.circrna/mm10/mm10_ref.gtf
mm10_fa=~/ref_annotation_Geneset/5.circrna/mm10/mm10.fa
nohup clear_quant -p 10 -1 SRR14078153_1.fastq.gz -2 SRR14078153_2.fastq.gz
-g $mm10_fa -i $mm10_hisat_index
-j $mm10_bowtie_index
-G $annotation_gtf
-o ./SRR14078153_output/ &
### 2.批量处理
mm10_hisat_index=/home/data/refdir/server/reference/index/hisat/mm10/genome
mm10_bowtie_index=~/ref_annotation_Geneset/5.circrna/mm10/bowtie1_index
annotation_gtf=/home/data/ssy43/ref_annotation_Geneset/5.circrna/mm10/mm10_ref.gtf
mm10_fa=~/ref_annotation_Geneset/5.circrna/mm10/mm10.fa
for i in `ls *_1.fastq.gz`
do
id=${i/_1.fastq.gz/}
output_dir=~/mouse_brain_rnaseq/3.cirRna/${id}_output
echo "clear_quant -p 40 -1 ${id}_1.fastq.gz -2 ${id}_2.fastq.gz
-g $mm10_fa -i $mm10_hisat_index
-j $mm10_bowtie_index
-G $annotation_gtf
-o $output_dir"
done >CIRCexplorer3.sh
## 开始了漫长的运行
nohup bash CIRCexplorer3.sh 1>nohup.log &
tree -h
最终的输出结果说明：
输出为 output_dir/quant/quant.txt：
Field | Description |
|---|---|
chrom | Chromosome |
start | Start of circular RNA |
end | End of circular RNA |
name | Circular RNA/Junction reads |
score | Flag of fusion junction realignment |
strand | + or - for strand |
thickStart | No meaning |
thickEnd | No meaning |
itemRgb | 0,0,0 |
exonCount | Number of exons |
exonSizes | Exon sizes |
exonOffsets | Exon offsets |
readNumber | Number of junction reads |
circType | Type of circular RNA |
geneName | Name of gene |
isoformName | Name of isoform |
index | Index of exon or intron |
flankIntron | Left intron/Right intron |
FPBcirc | Expression of circRNA |
FPBlinear | Expression of cognate linear RNA |
CIRCscore | Relative expression of circRNA |
结语
如前言部分所述，笔者认为，CIRCexplorer3用户不友好，Linux熟练度不够的同学慎重尝试。当然，想要折腾一下，尝试不同算法的同学，两个软件均可做尝试。此外，还有circ_find，是一款2013发表在nature上的经典Circular RNAs分析软件。
相关文章
- 【技术种草】cdn+轻量服务器+hugo=让博客“云原生”一下
- CLB运维&运营最佳实践 ---访问日志大洞察
- vnc方式登陆服务器
- 轻松学排序算法：眼睛直观感受几种常用排序算法
- 十二个经典的大数据项目
- 为什么使用 CDN 内容分发网络？
- 大数据——大数据默认端口号列表
- Weld 1.1.5.Final，JSR-299 的框架
- JavaFX 2012：彻底开源
- 提升as3程序性能的十大要点
- 通过凸面几何学进行独立于边际的在线多类学习
- 利用行动影响的规律性和部分已知的模型进行离线强化学习
- ModelLight：基于模型的交通信号控制的元强化学习
- 浅谈Visual Source Safe项目分支
- 基于先验知识的递归卡尔曼滤波的代理人联合状态和输入估计
- 结合网络结构和非线性恢复来提高声誉评估的性能
- 最佳实践丨云开发CloudBase多环境管理实践
- TimeVAE：用于生成多变量时间序列的变异自动编码器
- 具有线性阈值激活的神经网络：结构和算法
- 内网渗透之横向移动 -- 从域外向域内进行密码喷洒攻击