microRNA deep-sequencing数据分析手段 27

什么是microRNA呢?MicroRNAs(miRNAs)是在真核生物中长约20-22个碱基的对基因表达起负调控作用的RNA分子。它对于发育,细胞分化,分形,凋亡,染色体结构及病毒抗性都起着极其重要的作用。

最近也有研究表明在原核生物及病毒中也发现了类似机制的RNA。

microRNA是高度保守的。受这些高度保守的RNAs调控的基因比例非常高,约占基因组中1-5%的预测基因,10%(多的可至30%)的已知蛋白编码基因受到miRNAs的调控。

miRNAs的工作机制如下图所示:

miRNA分子是从核内的较大的转录子剪切而来。首先会被剪切成一个带发卡环的前体(hairpin precursors)。之后在两种RNaseIII内切酶(RNase III endonucleases)(Drosha和Dicer)的作用下切出发卡环的两段尾巴(发卡环的双链部分,长约~22碱基)以及发卡环本身。两段尾巴此时并未解链,被称为miRNA复合体(duplex)。解链后和目标DNA链结合的那段miRNA被称为成熟miRNA(mature miRNA),与之(不完全)互补另一段RNA被称为星型miRNA(miRNA*)。成熟miRNA与被调控单元形成RNA诱导的沉默复合体(RNA-induced silencing complex, RISC)对目标RNA起负调控作用。需要注意的是mature miRNA和miRNA*并不一定完全互补,下图显示的就是第一个在线虫中发现的miRNA的前体。

随着第二代测序技术(next-generation sequencing (NGS))的发展,人们开发了很多用于microRNAs深度测序分析的软件。这其中包括:

  • mireap: http://sourceforge.net/mailarchive/forum.php?forum_name=mireap-user
  • miRDeep: Friedlander,M.R., Chen,W., Adamidi,C., Maaskola,J., Einspanier,R., Knespel,S. and Rajewsky,N. (2008) Discovering microRNAs from deep sequencing data using miRDeep. Nat. Biotechnol., 26, 407–415.
  • miRanalyzer: Hackenberg,M., Sturm,M., Langenberger,D., Falcon-Perez,J.M. and Aransay,A.M. (2009) miRanalyzer: a microRNA detection and analysis tool for next-generation sequencing experiments. Nucleic Acids Res., 37, W68–W76.
  • miRExpress: Wang,W.C., Lin,F.M., Chang,W.C., Lin,K.Y., Huang,H.D. and Lin,N.S. (2009) miRExpress: analyzing high-throughput sequencing data for profiling microRNA expression. BMC Bioinformatics, 10, 328.
  • miRTRAP: Hendrix,D., Levine,M. and Shi,W. (2010) miRTRAP, a computational method for the systematic identification of miRNAs from high throughput sequencing data. Genome Biol., 11, R39.
  • DSAP: Huang,P.J., Liu,Y.C., Lee,C.C., Lin,W.C., Gan,R.R., Lyu,P.C. and Tang,P. (2010) DSAP: deep-sequencing small RNA analysis pipeline. Nucleic Acids Res., 38, W385–W391.
  • mirTools: Zhu,E., Zhao,F., Xu,G., Hou,H., Zhou,L., Li,X., Sun,Z. and Wu,J. (2010) mirTools: microRNA profiling and discovery based on high-throughput sequencing. Nucleic Acids Res., 38, W392–W397.
  • MIReNA: Mathelier,A. and Carbone,A. (2010) MIReNA: finding microRNAs with high accuracy and no learning at genome scale and from deep sequencing data. Bioinformatics, 26, 2226–2234.
  • miRNAkey: Ronen,R., Gan,I., Modai,S., Sukacheov,A., Dror,G., Halperin,E. and Shomron,N. (2010) miRNAkey: a software for microRNA deep sequencing analysis. Bioinformatics, 26, 2615–2616.

从历史的角度来讲,所有这些软件的先驱应该是miRDeep和mireap。其中miRDeep更是后面多种软件的基础而被广泛使用。本文出于篇辐的考虑也只介绍两个软件,miRDeep及miRanalyzer.

miRDeep是用perl写成的,具有跨平台的优势,容易为生物信息工作者使用。但是因为其缺少用户界面而不受普通生物科研人员的喜爱。而miRanalyzer是具有良好用户界面的服务器端应用,普通生物科研人员可以通过简单易用的几个选项得到较为清晰的分析结果。但是它需要事先将fastq文件转换成unique sequence count文件。当然这一步转换已经有很多工具可以使用,其实也未必就是难事。

但是要了解miRNA Deep-sequencing的分析机理,我们还需要深入地从miRDeep谈起。

在收集的总RNA样品中,用于测序的可能会是多种多样的小片段RNA,上图a所示的就是miRNA的深度测序的理想状态,miRNA前体被切割成三段,成熟miRNA, miRNA*以及发卡环,这时测序的话,会得到多个成熟miRNA的reads,较少量star sequence以及loop的reads,这是因为单链RNA更容易被降解。但是往往实验中是不可能如此理想的,还会有许多miRNA前体被测序,还有如图中b所示的非miRNA但具有与miRNA类似结构的RNA小分子的干扰。幸运的是这些小分子片段干扰因为缺少保护会如同start sequence及loop一样很少被测序到。由此,我们可以得到假阳性较低的结果。


在安装miRDeep时,我们需要同时安装

a     bowtie short read aligner (http://bowtie-bio.sourceforge.net/index.shtml)
b     Vienna package with RNAfold (http://www.tbi.univie.ac.at/~ivo/RNA/)
c.i   SQUID library (http://selab.janelia.org/software.html) goto Squid and download it
c.ii  randfold (http://bioinformatics.psb.ugent.be/software/details/Randfold)
d     Perl package PDF::API2 (http://search.cpan.org/search?query=PDF%3A%3AAPI2&mode=all)

miRDeep首先使用mapper.pl将fastq文件转换成它所需要的格式,其中hg19为bowtie主页上下载的index文件, 0hr.fastq为需要分析的文件。

mapper.pl 0hr.fastq -e -h -i -j -l 18 -m -p hg19 -s reads_collapsed.fa -t reads_collapsed_vs_genome.arf -v -n

而后使用quantifier.pl来分析样品的质量。这里的hsa_precursors.fa及hsa_mature.fa是从miRBase上下载而来。

quantifier.pl -p hsa_precursors.fa -m hsa_mature.fa -r reads_collapsed.fa -y now -t hsa

之后需要将rna模板库转换成dna

rna2dna.pl hsa_precursors.fa >  hsa_precursors.rna2dna.fa
rna2dna.pl hsa_mature.fa >  hsa_mature.rna2dna.fa

最后使用miRDeep2.pl分析数据。其中hg19.fa是从UCSC下载而来。

miRDeep2.pl reads_collapsed.fa hg19.fa reads_collapsed_vs_genome.arf hsa_mature.rna2dna.fa none hsa_precursors.rna2dna.fa -t Human -q expression_analyses/expression_analyses_now/miRBase.mrd 2>report.log

至此我们就可以得到结果了。我们需要注意的是,有一些counts很高的unique sequence并不一定会在结果中,有以下几种可能:
1,污染
2,miRDeep会抛弃在基因组上有超过5个比对匹配结果的序列
3,RNAfold预测的结构不符合标准的miRNA前体
4,其它。

题外话,在本文的写作过程中,阅读了一些文献,发现许多关于miRDeep的文献在讲述其机理时会有诸多谬误。所以本文在讲述的过程中都列有最初的参考文献,以避免因本人能力不足而误导了读者。

27 thoughts on “microRNA deep-sequencing数据分析手段

  1. Reply Emily 8月 21,2013 4:40 下午

    Hello, I am recently interested in tools in miRNA discovery, but I am not major in biology. Could I ask some simple questions? Thank you

  2. Reply zh 8月 25,2013 7:19 上午

    你好 我运行mapper.pl没几秒就报错了。mywork.err内容如下:
    parsing fastq to fasta format
    sh:dir_mapper_seq_/××××.fastq_0547599055_25_08_2013_t_19_49_49/reads.fa: No such file or directory converting rna to dna alphabet
    sh:dir_mapper_seq_/××××.fastq_0547599055_25_08_2013_t_19_49_49/reads_dna.fa: No such file or directory discarding sequences with non-canonical letters
    还有的省略。。。
    请问哈 这是咋回事啊

  3. Reply zh 8月 27,2013 4:24 上午

    非常感谢,我看了下mapper.pl的脚本,它的目录生成写的有点问题,我把测序数据和脚本放在同一目录下问题就解决了。

    还有个问题, hsa_precursors.fa 和 hsa_mature.fa这两个文件我在miRBase中怎么没有找到啊?

    谢谢

  4. Reply xutian 12月 23,2013 1:30 上午

    quantifier.pl 是用来做miRNA表达定量的。
    之后需要将rna模板库转换成dna: 这一步也不想要把

  5. Reply VictorZ 2月 13,2014 2:07 上午

    你好,我想问一下这篇文章似乎还有一半miRanalyzer的部分没有写嘛,博主准备写么?

    • Reply admin 2月 13,2014 9:15 上午

      miRanalyzer部分应该不会写了,因为明白了mirdeep,理解它就没什么问题了。
      关于fastq转换成unique sequence count的角本,之前miRanalyzer提供一个,但现在链接出现了问题。你可以试试使用:
      awk ‘NR%4==2’ in.fastq|sort |uniq -c|awk ‘{print $1″\t”$2}’ > out

  6. Reply VictorZ 2月 13,2014 2:12 上午

    另外,还有个位问题,就是有哪些软件能将fastq转换成unique sequence count文件?求推荐下,谢谢啦

  7. Reply VictorZ 2月 18,2014 9:17 下午

    谢谢博主前面的解答。
    另外,有个问题,就是如果我只想定量已知成熟miRNA的表达情况,quantifier.pl这步其实已经得到了一个结果,而后面miRDeep2.pl这一部又会得到一个成熟miRNA的表达结果,我对比了下两个结果,发现得到的成熟miRNA的read count(rc)基本相同,但normalize后的结果有所区别,我想问博主知道是什么原因造成的么? 还有miRDeep是怎么normalize原始rc的?

  8. Reply WS 3月 17,2014 10:26 下午

    Hi, 博主后续分析成熟的miRNA差异表达,用的是miRDeep哪部份的结果?是expressed_mature的还是最后results表格里的mature detected by miRDeep?前者可能有500–600个成熟miRNA,后者可能只有300个左右的miRNA(小鼠样本里)。results里给出的mature not detected by miRDeep是什么意思?是指这些成熟的miRNA虽然有reads比对上,但是根据miRDeep检测不认为这些符合cutoff要求么? 还有成熟的表达结果中有些mature miRNA的名字是相同的,可能来自不同的前体(不同的染色体),那么计算时如何做?把两部分的reads相加么?谢谢啦~

  9. Reply admin 3月 18,2014 9:12 上午

    结果应该以html文件中显示的为主。你的情况,应该以mature detected by miRDeep的为准。你可以试着把detected by miRDeep的结果加上not detected by miRDeep加起来看看,是否与expressed_mature的结果类似。如果map至不同的染色体,miRDeep会分别都报告一次。这个时候如果你一定要做统计,你需要加上权重。你也可以将mapping的严格度加强,设置g值为0也许可以去除大部分多重mapping的结果。但是也不会工作得特别完美。

  10. Reply obmldy 5月 9,2014 5:18 上午

    博主,请问miRDeep有何局限性

  11. Reply obmldy 5月 9,2014 5:23 上午

    还有,我看到有篇文献是index中的基因组才是需要分析的序列,而你说的0hr.fastq的作用是已经deep-sequencing的miRNA只是需要分解为read,然后mapping到index上去,再找到index上的miRNA序列,求解释到底实在分析哪个

  12. Reply obmldy 5月 11,2014 12:43 上午

    我的意思是说这里的是从hg19中检测已知和未知的miRNA

  13. Reply 丁丁 5月 11,2014 6:01 上午

    博主,你好~我跑完miRDeep2后得到的csv结果文件中最后not detected部分有大片的空白行存在,这个情况你碰到过吗?

  14. Reply obmldy 5月 13,2014 2:27 上午

    博主就是mirdeep对基因组和read的序列有什么要求吗,比如容错率之类的

    • Reply admin 5月 13,2014 10:28 上午

      miRdeep对基因组和read序列没有特殊要求。它设定的容错率为2。看原文:The retained reads were mapped to the dog genome (Canis familiaris version canFam2, from UCSC) using the custom mapping tool, allowing for up to two edits. For each read, mappings of suboptimal edit distance were discarded (if the best mapping was edit distance 1, all edit distance 2 mappings were discarded).

  15. Reply monkeyjun 7月 30,2014 9:45 上午

    博主老师好 ,我是一名生物信息的本科生,现在正在做miRNA-seq分析,想向您请教下这mirdeep2的结果因该怎么看,希望得到您的指导。谢谢!

Leave a Reply

  

  

  

%d 博主赞过: