希望所有的学徒,实习生,学员都能在生信技能树的舞台上绽放自己的光彩!
本次测评CircRNA-seq上游分析的两大最新工具CIRCexplorer3
及CIRIquant
。CIRCexplorer3
是2019年发表在Genomics Proteomics Bioinformatics(2020 IF=7.69)上,目前引用量是22次;CIRIquant
2020年发表在nature communications上,目前引用量是54次。
不考虑算法的前提下比较这两款软件:两款软件运行均比较慢,40个线程下双端测序的一个样本约需2小时。其中CIRCexplorer3
运行更慢一些,且需要安装非常多的依赖包。但是,CIRCexplorer3
github官网并没有列出这些依赖包,且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经验的同学来说可以略过。
软件下载使用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文件修剪
笔者使用的数据为项目数据,这里我们以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"
###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
###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 ./
原始数据已去接头。这里贴一下用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
官网列出的主要软件包括:
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
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 |
### 批量处理
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
该软件还提供了差异分析,总共三步走。
方便起见,这里仅简单做两组差异分析:
(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原始数据CIRCexplorer2
参考:https://circexplorer2.readthedocs.io/en/latest/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
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
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分析软件。
联系客服