miRNA-seq数据分析
miRNA-seq Data Analysis   

仪器设备

1. 服务器 (型号：I620-G20；操作系统：centos7；CPU：Xeon(R) CPU E5-2620 v4 @ 2.10GHz，8核，32G内存)
   注：学习者最好是有mac或者linux系统，需要8G+的内存，500G的存储，如果是Windows，则需要安装git，notepad++，虚拟机 (安装linux系统) 等,建议初学者使用mac。
2. 普通个人电脑 (需要已经安装好R，能调用个人终端)
   

软件

1. FASTX-Toolkit (v0.014)
2. FastQC (v0.11.8)
3. mulitiQC(v1.6)
4. miRDeep2(v2.0.0) (Anders和Huber，2010)
5. DESeq(v1.34.0) (Friedlander等，2012)
6. R(v3.5.1)
   

实验步骤

注：

1. miRNA-seq数据获得
   本实验方案以2017年发表在Oncotarget杂志上的文章：“Systematic analysis reveals tumor-enhancing and -suppressing microRNAs in Drosophila epithelial tumors.” (Shu等，2017) 中miRNA-seq的部分数据作为测试数据。
2. 去掉adapter
   在终端输入命令：$ fastx_clipper -a adapter -c –i输入文件名.fastq -z –o; 输出文件名.fastq.gz,操作如下图所示：
   
   
3. 通过FastQC对数据进行质量评估
   $ fastqc所在目录/fastqc -t 4输入文件
   注：若想指定输出目录，可以增加参数“-o路径”设置。
   操作如下图所示：
   
   在输出文件夹中每个输入文件都有两个输出结果，一个结果集合压缩包 (.zip文件)和一个网页结果 (.html文件)。一般通过查看网页结果即可知道基本的测序质量。
   结果如图1：
   
   
   图1. fastqc质控结果. A. 表示每个碱基的测序质量，在绿色区域表示测序质量很高；B. 展示了测序测到的长度分布，主要分布在21~23 nt，这与miRNA的大小符合。
   
   当具有多个样品时，查看多个fastqc的结果时可以使用multiQC软件，将所有的结果综合到一起查看，输入命令：
   $ multiqc目标文件夹
   查看方式与结果与fastqc是一致的，在此就不进行赘述。
4.  对测序文件进行比对
   目前比对miRNA的程序很多，比如miRDeep2、miRExpress、miRNAkey以及sRNAbench等，同时还有一些在线的分析软件。我们采用综合指标比较好的miRDeep2进行比对 (Bisgin等，2018)。
   
   1)

   建立index文件，把fastq文件转换为fasta文件，代码如下：
   $ bowtie-build基因组文件输出文件名
   注：这里的使用的是从flybase上下载的果蝇参考基因组文件。
   操作如图所示：
   
   

   2)

   将测序文件与参考基因组进行比对，比对命令如下：
   $ mapper.pl 输入文件-e-h -i-j -m -k adapter -l 18 -p参考基因组的index -s处理过的reads输出的文件名-t输出文件名.arf
   操作如图所示：
   
   其中所涉及到的参数具体含义入下：
   -c输入文件是fasta格式
   -h解析为fasta格式
   -i将rna转换为dna字母表 (以映射基因组)
   -j删除所有包含字母序列的条目，除了a，c，g，t，u，n，A，C，G，T，U，N之外
   -m collapse reads
   -k 3'适配器序列
   -l n忽视长度低于n的序列
   -p将处理过的reads map到之前建立过索引的基因组上，注意输入的是index文件的前缀
   -s 指出将处理过的reads输出到某个文件，自己命名
   -t 指出将mapping的结果输出到某个文件，自己命名，必须是.arf文件
   -o设置线程数
   运行结束后，在终端屏幕上会显示出一个比对结果的summary，如图2所示：
   
   
   图2. miRDeep2比对输出结果
   
   

   3)

   计算比对到已知miRNA的counts数
   $ quantifier.pl -p前体序列参考文件.fa -m成熟序列参考文件.fa -r上一步处理过的reads输出文件-s从miRBase上下载的star序列文件-t物种名称-y文件名后缀
   操作如图所示：
   
   运行结束后，会生成一个miRNAs_expressed_all_samples_now.csv文件。所有测到的已知miRNA的counts等信息都在这个文件中。我们后续的分析都是基于这个文件进行的。

   4)

   鉴定测序得到的未知miRNA
   一般情况下，鉴定新的基因不常用，但是如果需要也可以进行以下的操作：miRDeep2.pl处理过的reads输出文件基因组文件处理过的输出文件名.arf成熟miRNA文件其他物种的成熟miRNA文件研究物种miRNA前体的文件-t物种名称2>report.log
   如果只有reads，arf文件和genome文件，需要用none表示成熟miRNA文件、其他物种的成熟miRNA文件和研究物种miRNA前体的文件。
   操作如图所示：
   
   运行结束后，结果会保存在.html文件中。
5. miRNA表达矩阵以及差异表达分析
   miRNA可以使用DESeq (DESeq2)，edgeR以及limma等R包进行差异表达分析。由于示例数据是没有重复的数据，所以我们采用DESeq对数据进行分析 (DESeq2不支持没有重复组的分析，建议读者使用具有重复的实验数据)。将R的工作环境设置到含有输出的miRNAs_expression的文件夹中，运行如下脚本 (图3)：
   
   
   图3. DESEq分析R脚本
   
   运行脚本后我们可以获得miRNA差异基因表格以及差异最显著的前30个miRNA的热图 (如图4)。
   
   
   图4. miRNA差异基因热图
   
   
6. 下游分析
   目前进行miRNA测序数据分析，大多是要和mRNA测序数据关联来看。这里就涉及到了miRNA寻找对应的targets问题。目前寻找miRNA靶点的方法主要分为计算机预测和实验的方法。计算机预测靶点，可以参考该网站的介绍https://www.plob.org/article/1157.html。
   而实验方法，目前在果蝇方面，我们实验室首次实现了一种可以检测个体水平上的miRNA靶标的实验Ago1-RIP-seq，通过这种方法我们可以探究已知以及潜在的miRNA靶标，为miRNA的后续分析提供了非常坚实的实验指导。目前我们正在开发Ago1-RIP-seq的数据库，很快可以上线，为广大的科研人员提供服务。

参考文献

1. Anders, S. and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol 11(10): R106.
2. Bisgin, H., Gong, B., Wang, Y. and Tong, W. (2018). Evaluation of bioinformatics approaches for next-generation sequencing analysis of microRNAs with a toxicogenomics study design. Front Genet 9: 22.
3. Friedlander, M. R., Mackowiak, S. D., Li, N., Chen, W. and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res 40(1): 37-52.
4. Shu, Z., Huang, Y. C., Palmer, W. H., Tamori, Y., Xie, G., Wang, H., Liu, N. and Deng, W. M. (2017). Systematic analysis reveals tumor-enhancing and -suppressing microRNAs in Drosophila epithelial tumors. Oncotarget 8(65): 108825-108839.