CircRNA-seq上游分析工具测评：CIRIquant VS. CIRCexplorer3

生信技能树

发布于 2021-12-17 19:55:36

3.1K0

发布于 2021-12-17 19:55:36

文章被收录于专栏：生信技能树

前言

本次测评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分析软件。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-12-16，如有侵权请联系 cloudcommunity@tencent.com 删除

node.js

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度