circRNA_finder是一款环状RNA预测软件，在对果蝇的研究中采用该软件进行了环状RNA的预测，该软件的源代码托管在github上，网址如下
软件的安装比较简便，直接下载解压就可以了，需要注意的是，该软件依赖以下几个软件
- perl
- awk
- STAR
- samtools
其中samtools的版本必须是低于1.0的版本，因为两个版本samtools sort的用法有变化。软件的使用包含如下两个步骤
虽然软件提供了一个名为runStar.pl的脚本(如下)，但是由于STAR的版本问题，使用起来并不方便。该脚本本质上是对STAR的封装，直接用STAR就好了，参数设置可以参考脚本中的设置。
./runStar.pl --inFile1 [R1 fastq] --inFile2 [R2 fastq] --genomeDir [path to STAR genome] --maxMismatch [max mismatches realtive to read length, default 0.02] --outPrefix [output directory and prefix]
1. STAR构建index
STAR --runMode genomeGenerate \
--runThreadN <# cpus> \
--genomeDir <genome output directory> \
--genomeFastaFiles <input Genome FASTA file>
- --runThreadN是指你要用几个cpu来运行；
- --genomeDir构建索引输出文件的目录；
- --genomeFastaFiles你的基因组fasta文件所在的目录
2. STAR比对
首先要明确的是，STAR软件的参数实在是太多了，那么经常用的就那么几个。
下面的代码里，其中--twopassMode参数凭个人爱好是否使用，这个参数是非常非常耗时的，加上这个参数以后的运行时间是正常的两倍。而且非常消耗内存（star的twopassMode问题）。但是比对结果更为精准。
$ srun --mem=128G --cpus-per-task=20 --time=2:00:00 --pty STAR \
--runThreadN 40 \ #线程数
--runMode alignReads \#比对模式
--readFilesCommand zcat \#说明你的fastq文件是压缩形式的，就是.gz结尾的，不加的话会报错
--quantMode TranscriptomeSAM GeneCounts \ #将reads比对至转录本序列
--twopassMode Basic \#先按索引进行第一次比对，而后把第一次比对发现的新剪切位点信息加入到索引中进行第二次比对。这个参数可以保证更精准的比对情况，但是费时也费内存。
--outSAMtype BAM Unsorted \ #输出BAM文件，不进行排序。如果不加这一行，只输出SAM文件。
--outSAMunmapped None \
--genomeDir /gpfs/home/fangy04/downloads/STAR_index/GRCh38/ \ #索引文件目录
--readFilesIn /gpfs/home/fangy04/downloads/SRR8112732_1.fastq.gz /gpfs/home/fangy04/downloads/SRR8112732_2.fastq.gz \ #两个fastq文件目录
--outFileNamePrefix DRB_TT_seq_SRR8112732 #输出文件前缀
因为后续运行ciriRNA_finder 除了正常比对结果还需要以下三个文件：
所以在STAR 比对中最少需要增加两个参数：--chimOutType Junctions SeparateSAMold --chimSegmentMin 10;--chimSegmentMin默认为0，如果为0着不会输出junction信息文件
用了服务器里128G只比对了一对fastq文件，用时20分钟，嗯，大家可以体验一下~
STAR比对后会输出好几个文件。上面的代码里，如果你没有加--outSAMtype BAM Unsorted这一个参数的设置，只会输出SAM。在这些文件里，最重要的就是“*.Aligned.out.sam”，你可以使用samtools对其进行bam文件的转换以及排序。但是如果你设置了输出文件为bam，那么直接用samtools进行排序就好了。
输出的文件：
9216920116 Jun 28 17:06 DRB_TT_seq_SRR8112732Aligned.out.bam #这个文件是最重要的，用来后续进行remove duplicates和sort
1166235552 Jun 28 17:06 DRB_TT_seq_SRR8112732Aligned.toTranscriptome.out.bam #这个文件是那些比对到转录本上的reads组成的bam文件
2034 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.final.out
20188 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.out
2571 Jun 28 17:06 DRB_TT_seq_SRR8112732Log.progress.out
1585521 Jun 28 17:06 DRB_TT_seq_SRR8112732ReadsPerGene.out.tab
6732305 Jun 28 17:06 DRB_TT_seq_SRR8112732SJ.out.tab #剪切的信息
8192 Jun 28 16:51 DRB_TT_seq_SRR8112732_STARgenome
8192 Jun 28 16:51 DRB_TT_seq_SRR8112732_STARpass1
如果你想查看比对率之类的信息，可以查看*.Log.final.out文件：
$ cat DRB_TT_seq_SRR8112732Log.final.out
Started job on | Jun 28 16:44:01
Started mapping on | Jun 28 16:54:06
Finished on | Jun 28 17:06:52
Mapping speed, Million of reads per hour | 369.36
Number of input reads | 78590839 #fastq文件的信息
Average input read length | 150 #read长度
UNIQUE READS:
Uniquely mapped reads number | 58096443 #唯一比对上的reads数量
Uniquely mapped reads % | 73.92%
Average mapped length | 144.37
Number of splices: Total | 1880228
Number of splices: Annotated (sjdb) | 1554548 #剪切数
Number of splices: GT/AG | 1513028
Number of splices: GC/AG | 50306
Number of splices: AT/AC | 4877
Number of splices: Non-canonical | 312017 #非典型剪切数
Mismatch rate per base, % | 0.73%
Deletion rate per base | 0.01%
Deletion average length | 1.61
Insertion rate per base | 0.01%
Insertion average length | 1.63
MULTI-MAPPING READS: #多重比对数
Number of reads mapped to multiple loci | 4439909
% of reads mapped to multiple loci | 5.65%
Number of reads mapped to too many loci | 405437
% of reads mapped to too many loci | 0.52%
UNMAPPED READS: #未比对上的reads
Number of reads unmapped: too many mismatches | 0
% of reads unmapped: too many mismatches | 0.00%
Number of reads unmapped: too short | 15355117
% of reads unmapped: too short | 19.54%
Number of reads unmapped: other | 293933
% of reads unmapped: other | 0.37%
CHIMERIC READS: #嵌合的reads数
Number of chimeric reads | 0
% of chimeric reads | 0.00%
而输出文件里的*SJ.out.tab文件，则是splice junctions的一些信息，其中需要注意的是：对于junction的位置信息，STAR则是按照intron的起始和终止位置来定，而其他的一些软件则是按照exon的位置来决定的（参考:比对软件STAR的简单使用）：
$ head DRB_TT_seq_SRR8112732SJ.out.tab
chr1 10078 10178 2 2 0 0 1 26
chr1 10084 10178 2 2 0 0 1 32
chr1 10090 10178 2 2 0 0 1 32
chr1 10102 10178 2 2 1 0 4 37
chr1 10113 10183 2 2 1 0 1 37
chr1 10137 10282 2 2 1 0 1 30
chr1 10137 10288 2 2 1 1 1 37
chr1 10137 10294 2 2 1 0 1 30
chr1 10162 10397 2 2 1 0 1 13
chr1 10179 10225 2 2 1 3 2 27
#第一列：染色体名称
#第二列：内含子起始位点
#第三列：内含子结束位点
#第四列：链。1是正链，2是负链。
#第五列：内含子的motif: 0: 非典型; 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5:AT/AC, 6: GT/AT
#第六列：0: 未注释的； 1：注释的
#第七列：横跨这个junction上的唯一比对上的reads数
#第八列：横跨这个junction上的muoti-mapping的reads数
#第九列：maximum spliced alignment overhang
3. 运行circRNA_finder
perl \
postProcessStarAlignment.pl \
--starDir star_out_dir \
--minLen 100
--outDir output_dir
运行完成之后，会三个文件,对应的后缀如下所示
- _filteredJunctions.bed
- _s_filteredJunctions.bed
- _s_filteredJunctions_fw.bed
第一个文件为所有环状RNA的结果文件；第二个文件为剪切位点符合GT-AG剪切信号的环状RNA；第三个文件和第二个文件的环状RNA相同，只不过新增了环状RNA连接点附近的线性RNA平均测序深度信息。通常情况下，我们选择第二个文件的结果作为最终的环状RNA预测结果，该文件内容示意如下
每一行代表一个环状RNA,第一列代表染色体编号，第二列和第三列分别代表起始和终止位置，第四列代表name,这里用数字编号加上字母s表示，第五列代表环状RNA的junction reads数目,也就是表达量，第六列代表正负链信息。
参考资料：
