方便快捷的smallRNA数据分析流程

生信店小二 · 简书 · · 2020-11-22 11:30

文章预览

前言

今天来跟大家分享一个smallRNA的分析流程—— sRNAnalyzer 。这个流程是由perl语言写的，对于老生信人来说应该很熟悉，毕竟perl以前也是很受欢迎的编程语言，不过对于近几年入行生信的来说就有点陌生了，因为现在做生信分析大多使用的是python和R语言了。所以对于这个流程的源码我也是看的似懂非懂的状态，因为我对perl也是一点不懂。那为什么要分享这个流程呢？首先是这个流程可扩展性比较好，就是你可以自己定义所使用的参考基因组的数量和顺序；再者就是定量miRNA时直接定量为前体或者成熟的miRNA，相当方便；还有就是在定量的时候可以考虑错配的碱基数目，可以分别统计错配0、1、2个碱基的定量情况。该流程做了很多工作还有一些其他额外的功能，大家可以根据需要自己去探索一下，我只是使用了其中定量的功能，下面也就定量的功能来做一个介绍。

准备工作

在使用这个流程之前先保证你的工作环境中已经安装好了 fastx_toolkit 、 cutadapt 、 bowtie 这三个软件，流程中的数据预处理和比对会用到这些软件。下面还需要准备参考基因组，如果你没有特别的要求可以直接去软件的官网直接下载：
http://srnanalyzer.systemsbiology.net/ 。到官网后你可以看到如下图所示的表格，第一个链接是流程的代码，剩下的是配套的参考基因组链接，如果只是做简单的定量，只需下载smallRNA数据库文件 sRNA_DBs 就可以了(该文件包含了GtRNAdb 、miRBase 、MirGeneDB 、piRBase snoRNABase参考基因组)，其他的数据库文件挺大的下载起来需要很长时间。如果你觉得参考基因组不够或不是自己想要的，你也可以直接用bowtile软件自己构建，注意是 bowtile 不是bowtile2。

数据预处理

准备好软件和参考基因组，下面可以开始分析流程了，首先就是数据预处理，做这个步骤时，需要填写一个流程的参数配置文件：
pipeline_config.conf：

preprocess:
  kit:        NEB
  gzip:       true
  stop-oligo: false
        
alignment:
  type: single
  human_miRNA:     2
  human_miRNA_sub: 2
  human_piRNA:     2
  human_snoRNA:    2

可以看出来该配置文件包含了数据预处理和比对两个步骤的参数设置，preprocess设置的是数据预处理的参数：
kit: NEB，设置建库试剂盒，据此来设置不同的接头序列，还有两个参数为Illumina、4N
gzip: true，设置fastq是否为压缩格式
stop-oligo: false，设置是否有结束序列

alignment设置比对时的参数:
type: single，设置数据是单端还是双端
human_miRNA: 2 #设置比对时最大的错配碱基数
human_miRNA_sub: 2
human_piRNA: 2
human_snoRNA: 2
设置好配置文件，可以使用脚本做预处理了，将所有需要处理的fastq文件放置在一个文件夹里面，然后进入到该文件夹里面，运行如下的代码，脚本就会将所有的文件都处理好，会在该文件加下生成处理好的文件：

preprocess.pl --config pipeline_config.conf

结果类似如下：

GSM3593646_Cutadapt_Prinseq.report
GSM3593646_Cutadapt.report
GSM3593646.fastq.gz
GSM3593646_Processed.fa

比对及定量

当做好数据预处理步骤后，就可以做比对及定量了，不过在此之前需要两个配置文件，一个是上面准备的配置文件，另一个是数据库配置文件，如下所示：
DB_config.conf:

base: path/to/database
human_miRNA:     miRBase/hairpin_hsa_anno
human_miRNA_sub: miRBase/hairpin_hsa_sub_anno
human_piRNA:     piRBase/piR_human_v1.0
human_snoRNA:    snoRNABase/snoRNABase
virus_miRNA:     miRBase/hairpin_virus_anno
plant_miRNA:     miRBase/hairpin_plant_anno
all_miRNA:       miRBase/hairpin_anno
all_miRNA_sub:   miRBase/hairpin_sub_anno
MirGeneDB_human_miRNA:  MirGeneDB/MirGeneDB_hsa_precursor_anno
GtRNAdb_human_tRNA:     GtRNAdb/hg38-tRNAs
GtRNAdb_bacteria_tRNA:  GtRNAdb/bacterial-tRNAs
GtRNAdb_all_tRNA:       GtRNAdb/GtRNAdb-all-tRNAs

要注意path/to/database是所有数据库的共同路径，例如上面所示的human_miRNA参考基因组的绝对路径就是path/to/database/miRBase/hairpin_hsa_anno。还有一点需要注意的是，数据库的顺序也是有意义的，因为在比对时，会按照这里面的顺序逐个去比对，没有比对上的read会比对下一个参考基因组，直到所有的参考基因组都比对一遍。这样有一个好处就是read只会比对到一个参考基因组，不存在同时比对到多个基因组的情况。
准备好配置文件，就可以做比对及定量分析了，代码如下：

align.pl fastadir pipeline_config.conf DB_config.conf

fastadir：数据预处理的文件夹
pipeline_config.conf：数据预处理时的流程配置文件
DB_config.conf：数据库配置文件
做完比对和定量后，会在文件加下面多生成一些文件，类似如下所示：

GSM3593663_Cutadapt_Prinseq.report
GSM3593663.fastq.gz        
GSM3593663_Processed.fa              #处理后数据，fasta格式
GSM3593663_Processed.profile       #记录了每一条read比对结果
GSM3593663_Cutadapt.report
GSM3593663_Processed.dist           #reads长度分布文件
GSM3593663_Processed.feature     #定量结果
GSM3593663_Processed_unMatch.fa   #没有比对上参考基因组的reads

到此定量就完成了，但是流程为了方便大家，还准备了一个汇总结果的脚本，也就是说如果有很多样本，该脚本会将所有样本的每一种结果都汇总到单独的一个文件里面，代码如下：

summarize.pl DB_config.conf --project smallrna

会在文件夹下面生成以下的文件：

smallrna_des.feature
smallrna_distCount.sum
smallrna.feature
smallrna_matchRate.sum
smallrna_stp.sum
smallrna_des.profile
smallrna.distRate.sum
smallrna_matchCount.sum
smallrna.profile
smallrna_trimRate.sum

这样所有的样本的结果都汇总在一个文件中，方便后续的使用。到此，定量就全部完成了。其实这个流程还有一个特别之处，就是miRNA的参考基因组，使用的全部是前体的序列，只不过在head部分保留了成熟miRNA序列在前体中的位置信息，然后根据比对的结果落在成熟miRNA区域在整个read长度占比情况来判断属于前体或者成熟的miRNA，也就是说如果与成熟miRNA的区域重叠超过60%则属于成熟的miRNA，否则属于前体。
下面展示一下经过处理的miRNA参考基因组的fasta文件：

>hsa-let-7a-1|0:79||hsa-let-7a-1-5p|5:26||hsa-let-7a-1-3p|56:76| MI0000060 Homo sapiens let-7a-1 stem-loop
TGGGATGAGGTAGTAGGTTGTATAGTTTTAGGGTCACACCCACCACTGGGAGATAACTATACAATCTACTGTCTTTCCTA
>hsa-let-7b|0:82||hsa-let-7b-5p|5:26||hsa-let-7b-3p|59:80| MI0000063 Homo sapiens let-7b stem-loop
CGGGGTGAGGTAGTAGGTTGTGTGGTTTCAGGGCAGTGATGTTGCCCCTCGGAAGATAACTATACAACCTACTGCCTTCCCTG
>hsa-let-7c|0:83||hsa-let-7c-5p|10:31||hsa-let-7c-3p|55:76| MI0000064 Homo sapiens let-7c stem-loop
GCATCCGGGTTGAGGTAGTAGGTTGTATGGTTTAGAGTTACACCCTGGGAGTTAACTGTACAACCTTCTAGCTTTCCTTGGAGC

经过处理miRNA参考基因组中虽然保留了成熟miRNA的信息，但是需要注意一点，就是结果中成熟miRNA的ID可能与 miRBase 数据库中不完全一致，例如“hsa-mir-101-1”前体的成熟miRNA在miRBase数据库中的ID是“hsa-miR-101-3p”，但在sRNAnalyzer流程的结果中成熟miRNA的ID是“hsa-miR-101-1-3p”。可能作者在处理miRNA参考基因组时，成熟miRNA的ID是根据直接前体ID来生成的而不是用的数据库中的名称，这一点大家一定要留意一下！！！

最后

这个流程用来定量还是很方便的，如果想要做novel miRNA的预测可以结合使用 miRDeep2 。今天就分享到这里了~~~

………………………………

原文地址：访问原文地址
快照地址：访问文章快照
总结与预览地址：访问总结与预览

分享到微博