一篇文章学会miRNA-seq分析

第一讲：文献选择与解读

前阵子逛BioStar论坛的时候看到了一个关于miRNA分析的问题，提问者从NCBI的SRA中下载文献提供的原始数据，然后处理的时候出现了问题。我看到他列出的数据来自iron torrent测序仪，而且我以前也没有做过miRNA-seq的数据分析， 就自学了一下。因为我有RNA-seq的基础，所以理解学习起来比较简单。

在这里记录自己的学习过程，希望对需要的朋友有帮助。

这里选择的文章是2014年发表的，作者用ET-1刺激human iPSCs (hiPSC-CMs) 细胞前后，观察miRNA和mRNA表达量的变化，我并没有细看文章的生物学意义，仅仅从数据分析的角度解读一下这篇文章，mRNA表达量用的是Affymetrix Human Genome U133 Plus 2.0 Array，分析起来特别容易。得到表达矩阵，然后用limma这个包找差异表达基因即可。

但是miRNA分析起来就有点麻烦了，作者用的是iron torrent测序仪。不过本次从SRA数据中心下载的是已经去掉接头的fastq格式测序数据，所以这里其实并不需要考虑测序仪的特异性。

关于该文章的几个资料

  • paper : http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0108051

  • Aggarwal P, Turner A, Matter A, Kattman SJ et al. RNA expression profiling of human iPSC-derived cardiomyocytes in a cardiac hypertrophy model. PLoS One 2014;9(9):e108051. PMID: 25255322

  • The accession numbers are 1. SuperSeries (mRNA+miRNA) - GSE60293

  • mRNA expression array - GSE60291 (Affymetrix Human Genome U133 Plus 2.0 Array)

  • miRNA-Seq - GSE60292 (Ion Torrent)

  • GEO : http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60292

  • FTP : ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP045/SRP045420

接下来我们要知道文章做了哪些分析，然后才能自己模仿看是否可以得到同样的分析结果。

文章数据处理流程

  • Ion Torrent's Torrent Suite version 3.6 was used for basecalling

  • Raw sequencing reads were aligned using the SHRiMP2 aligner and were aligned against the human reference genome (hg19) for novel miRNA prediction and then against a custom reference sequence file containing miRBase v.20 known human miRNA hairpins, tRNA, rRNA, adapter sequences and predicted novel miRNA sequences.(Genome_build: hg19, miRBase v.20 human miRNA hairpins)

  • The miRDeep2 package (default parameters) was used to predict novel (as yet undescribed) miRNAs

  • Alignments with less than 17 bp matches and a custom 3′ end phred q-score threshold of 17 were filtered out.

  • miRNA quanitification was done using HTSeq v0.5.3p3 using the default union parameter. Differential miRNA expression was analyzed using the DESeq (v.1.12.1) R/Bioconductor package

  • In this study, differentially expressed genes that had a false discovery rate cutoff at 10% (FDR< = 0.1), a log2 fold change greater than 1.5 and less than −1.5 were considered significant.

  • Target gene prediction was performed using the TargetScan (version 6.2) database

  • We also used miRTarBase (version 4.3), to identify targets that have been experimentally validated

  • miR-Deep2 and miReap predict exact precursor sequence according from mature sequence

    文章提到了fastq数据质量控制标准，数据比对工具，比对的参考基因组（两条比对线路），获得miRNA表达量，miRNA预测，miRNA靶基因预测。

这也是我们学习miRNA-seq的数据分析的标准套路， 而且作者给出了所有的分析结果，我们完全可以通过自己的学习来重现他的分析过程。

  • Supplementaryfilesformatandcontent: tab-delimited text files containing raw read counts for known mature human miRNAs.（表达矩阵）

  • We detected 836 known human mature miRNAs in the control-CMs and 769 in the ET1-CMs

  • Based on our miRNA-Seq data, we predicted 506 sequences to be potentially novel, as yet undescribed miRNAs.

  • In order to validate the expression profiles of the miRNAs detected, we performed RT-qPCR on a subset of five known human mature and five of our predicted novel miRNAs.

  • we obtained a total of 1,922 predicted miRNA-mRNA pairs represented by 309 genes and 174 known mature human miRNAs.

当然仅仅是套路分析还不够，所以他进行了 miRNA和mRNA 进行网络分析并做了少量湿实验来验证，最后还扯了一些生物学意义。

第二讲：搜集学习资料

因为我是完全从零开始入门miRNA-seq分析，所以收集的资料比较齐全。

首先看了部分中文资料，了解miRNA测序是怎么回事，该分析什么，然后主要围绕着第一讲文献里的分析步骤来搜索资料。

miRNA定义

我首先找到了miRNA定义：http://nar.oxfordjournals.org/content/34/suppl_1/D135.full

MicroRNAs (miRNAs) are small RNA molecules, which are ∼22 nt sequences that have an important role in the translational regulation and degradation of mRNA by the base's pairing to the 3′-untranslated regions (3′-UTR) of the mRNAs. The miRNAs are derived from the precursor transcripts of ∼70–120 nt sequences, which fold to form as stem–loop structures, which are thought to be highly conserved in the evolution of genomes. Previous analyses have suggested that ∼1% of all human genes are miRNA genes, which regulate the production of protein for 10% or more of all human coding genes。

选择参考序列

然后我比较纠结的问题是参考序列如何选择，因为miRNA序列很少，把它map到3G大小的人类基因组有点浪费计算资源，正好我的服务器又坏了，不想太麻烦，想用自己的个人电脑搞定这个学习过程。

我看到很多帖子提到的都是用bowtie比对到参考miRNA数据库(miRNA count: 28645 entries)http://www.mirbase.org/ ，从这个数据库，我明白了前体miRNA和成熟的miRNA的区别，前体miRNA长度一般是70–120 碱基，一般是茎环结果，也就是发夹结构，所以叫做hairpin。成熟之后，一般22 个碱基，在miRNA数据库很容易下载到这些数据，目前人类这个物种已知的成熟miRNA共有2588条序列，而前体miRNA共有1881条序列，我下载（下载时间2016年6月 ）的代码如下所示：

  1. wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.gz   ##　28645　reads

  2. wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.zip   ##   35828 reads

  3. wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.zip

  4. wget ftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3 ##

  5. wget ftp://mirbase.org/pub/mirbase/CURRENT/miFam.dat.zip

  6. grep sapiens mature.fa |wc  　# 2588

  7. grep sapiens hairpin.fa |wc       # 1881

  8. ## Homo sapiens

  9. perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' hairpin.fa >hairpin.human.fa

  10. perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' mature.fa >mature.human.fa

  11. # 这里值得一提的是miRBase数据库下载的序列，居然都是用U表示的，也就是说就是miRNA序列，而不是转录成该miRNA的基因序列，而我们测序的都是基因序列。

通过这个代码制作的hairpin.human.fa 和 mature.human.fa 就是本次数据分析的参考基因组。

在搜集资料的过程中，我看到了一篇文献讲挖掘1000genomes的数据找到位于miRNA的snp位点

https://genomemedicine.biomedcentral.com/articles/10.1186/gm363

看起来比较新奇，不过跟本次学习过程没有关系，我就是记录一下，有空回来学习学习。

博客讲解如何分析miRNA数据

  • http://genomespot.blogspot.com/2013/08/quick-alignment-of-microrna-seq-data-to.html

公司数据分析流程：

  • http://bioinfo5.ugr.es/miRanalyzer/miRanalyzer_tutorial.html

  • http://www.partek.com/sites/default/files/Assets/UserGuideMicroRNAPipeline.pdf

  • http://partek.com/Tutorials/microarray/microRNA/miRNA_tutorial.pdf

  • http://www.arraystar.com/reviews/microrna-sequencing-data-analysis-guideline/

  • http://bioinfo5.ugr.es/sRNAbench/sRNAbench_tutorial.pdf

  • http://seqcluster.readthedocs.io/mirna_annotation.html

耶鲁大学

  • http://www.yale.edu/giraldezlab/miRNA.html

南方基因

  • http://www.southgene.com/newsshow.php?cid=55&id=73

miRNA研究整套方案

  • http://wenku.baidu.com/view/5f38577a31b765ce05081429.html?re=view

Biostar 讨论帖子

  • https://www.biostars.org/p/3344/

  • https://www.biostars.org/p/98486/

miRNA-seq数据处理实战指南

  • http://bib.oxfordjournals.org/content/early/2015/04/17/bib.bbv019.full

直接用一个包搞定

  • http://bioconductor.org/packages/release/bioc/html/easyRNASeq.html

github流程：miRNA Analysis Pipeline v0.2.7

  • https://github.com/bcgsc/mirna/tree/master/v0.2.7

  • https://tools.thermofisher.com/content/sfs/manuals/CO25176_0512.pdf

miRNA annotation

  • http://seqcluster.readthedocs.io/mirna_annotation.html

网页版分析工具

  • https://wiki.uio.no/projects/clsi/images/2/2f/HTS2014miRNAanalysisLifeportal14final.pdf

  • http://www.training.prace-ri.eu/uploads/tx_pracetmo/NGSdataAnalysisWithChipster.pdf

可视化IGV User Guide

  • http://www.broadinstitute.org/igv/book/export/html/6

比较特殊的是新的miRNA预测，miRNA靶基因预测，这块软件太多并没有成型的流程和标准。

第三讲：下载公共数据

前面已经讲到了该文章的数据已经上传到NCBI的SRA数据中心，所以直接根据索引号下载，然后用SRAtoolkit转出我们想要的fastq测序数据即可。

下载的数据一般要进行质量控制，可视化展现一下质量如何，然后根据大题测序质量进行简单过滤。所以需要提前安装一些软件来完成这些任务，包括：sratoolkit /fastx_toolkit /fastqc/bowtie2/hg19/miRBase/SHRiMP

下面是我用新服务器下载安装软件的一些代码记录，因为fastx_toolkit /fastqc我已经安装过，就不列代码了

  1. ## pre-step: download sratoolkit /fastx_toolkit_0.0.13/fastqc/bowtie2/hg19/miRBase/SHRiMP

  2. ## http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software

  3. ## http://www.ncbi.nlm.nih.gov/books/NBK158900/

  4. ## 我这里特意挑选的二进制版本程序下载的，这样直接解压就可以用，但是需要挑选适合自己的操作系统的程序。

  5. cd ~/biosoft

  6. mkdir sratoolkit &&  cd sratoolkit

  7. wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.6.3/sratoolkit.2.6.3-centos_linux64.tar.gz

  8. ##

  9. ##  Length: 63453761 (61M) [application/x-gzip]

  10. ##  Saving to: "sratoolkit.2.6.3-centos_linux64.tar.gz"

  11. tar zxvf sratoolkit.2.6.3-centos_linux64.tar.gz

  12. cd ~/biosoft

  13. mkdir bowtie &&  cd bowtie

  14. wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip/download

  15. #Length: 27073243 (26M) [application/octet-stream]

  16. #Saving to: "download"

  17. mv download  bowtie2-2.2.9-linux-x86_64.zip

  18. unzip bowtie2-2.2.9-linux-x86_64.zip

  19. ## http://compbio.cs.toronto.edu/shrimp/

  20. mkdir SHRiMP &&  cd SHRiMP

  21. wget http://compbio.cs.toronto.edu/shrimp/releases/SHRiMP_2_2_3.lx26.x86_64.tar.gz

  22. tar zxvf SHRiMP_2_2_3.lx26.x86_64.tar.gz

  23. cd SHRiMP_2_2_3

  24. export SHRIMP_FOLDER=$PWD  ## 这个软件使用的时候比较奇葩，需要设置到环境变量，不能简单的调用全路径

SHRiMP这个软件比较小众，我也是第一次听说过。

本来我计划是能用bowtie搞定，但是第一次比对出了一个bug，就是下载的miRNA序列里面的U没有转换成T，所以导致比对率非常之低。于是我不得不根据文章里面记录的软件SHRiMP 来做比对，最后发现比对率完全没有改善，搞得我都在怀疑是不是作者乱来了。

下面是下载数据，质量控制的代码，希望大家可以照着运行一下。

  1. ## step1 : download raw data

  2. mkdir miRNA_test && cd miRNA_test

  3. echo {14..19} |sed 's/ /\n/g' |while read id; \

  4. do  wget "ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP045/SRP045420/SRR15427$id/SRR15427$id.sra"  ;\

  5. done

  6. ## step2 :  change sra data to fastq files.

  7. ## 主要是用shell脚本来批量下载

  8. ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done

  9. rm *sra

  10. ##  33M --> 247M

  11. #Read 1866654 spots for SRR1542714.sra

  12. #Written 1866654 spots for SRR1542714.sra

  13. ## step3 : download the results from paper

  14. ## http://www.bio-info-trainee.com/1571.html

  15. ## ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1009/suppl/GSE1009_RAW.tar

  16. mkdir paper_results && cd paper_results

  17. wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60292/suppl/GSE60292_RAW.tar

  18. ## tar xvf GSE60292_RAW.tar

  19. ls *gz |while read id ; do (echo $id;zcat $id | cut -f 2 |perl -alne '{$t+=$_;}END{print $t}');done

  20. ls *gz |xargs gunzip

  21. ## step4 : quality assessment

  22. ls *fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done

  23. ## Sequence length 8-109

  24. ## %GC 52

  25. ## Adapter Content passed

  26. ## write a script : :: cat >filter.sh

  27. ls *fastq |while read id

  28. do

  29. echo $id

  30. ~/biosoft/fastx_toolkit_0.0.13/bin/fastq_quality_filter -v -q 20 -p 80 -Q33  -i $id -o tmp ;

  31. ~/biosoft/fastx_toolkit_0.0.13/bin/fastx_trimmer -v -f 1 -l 27 -i tmp  -Q33 -z -o ${id%%.*}_clean.fq.gz ;

  32. done

  33. rm tmp

  34. ## discarded 12%~~49%%

  35. ls *_clean.fq.gz | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done

  36. mkdir QC_results

  37. mv *zip *html QC_results

这个代码是我自己根据文章的理解写出的，因为我本身不擅长miRNA数据分析，所以在进行QC的时候参数选择可能并不是那么友好.

  1. ~/biosoft/fastx_toolkit_0.0.13/bin/fastq_quality_filter -v -q 20 -p 80 -Q33  -i $id -o tmp ;

  2. ~/biosoft/fastx_toolkit_0.0.13/bin/fastx_trimmer -v -f 1 -l 27 -i tmp  -Q33 -z -o ${id%%.*}_clean.fq.gz ;

最后得到的clean.fq.gz系列文件，就是我需要进行比对的序列。

第四讲：测序数据比对

序列比对是大多数类型数据分析的核心，如果要利用好测序数据，比对细节非常重要，我这里只是研读一篇文章也就没有对比对细节过多考虑，只是列出自己的代码和自己的几点思考，力求重现文章作者的分析结果。

对miRNA-seq数据有两条比对策略

  1. 下载miRBase数据库里面的已知miRNA序列来进行比对

  2. 直接比对到参考基因组(比如人类的是hg19/hg38)

前面的比对非常简单，而且很容易就可以数出已经的所以miRNA序列的表达量；后面的比对有点耗时，而且算表达量的时候也不是很方便，但是它的优点是可以来预测新的miRNA，所以大多数文章都会把这两条路给走一下。

本文选择的是SHRiMP这个小众软件，起初我并没有在意，就用的bowtie2而已，参考基因组就用了miRBase数据库下载的人类的参考序列。

  1. ## step5 : alignment to miRBase v21 (hairpin.human.fa/mature.human.fa )

  2. #### step5.1 using bowtie2 to do alignment

  3. mkdir  bowtie2_index &&  cd bowtie2_index

  4. ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ../hairpin.human.fa hairpin_human

  5. ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ../mature.human.fa  mature_human

  6. ls *_clean.fq.gz | while read id ; do  ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -x miRBase/bowtie2_index/hairpin_human -U $id   -S ${id%%.*}.hairpin.sam ; done

  7. ## overall alignment rate:  10.20% / 5.71%/ 10.18%/ 4.36% / 10.02% / 4.95%  (before convert U to T )

  8. ## overall alignment rate:  51.77% / 70.38%/51.45% /61.14%/ 52.20% / 65.85% (after convert U to T )

  9. ls *_clean.fq.gz | while read id ; do  ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -x miRBase/bowtie2_index/mature_human  -U $id   -S ${id%%.*}.mature.sam ; done

  10. ## overall alignment rate:  6.67% / 3.78% / 6.70% / 2.80%/ 6.55% / 3.23%    (before convert U to T )

  11. ## overall alignment rate:  34.94% / 46.16%/ 35.00%/ 38.50% / 35.46% /42.41%(after convert U to T )

  12. #### step5.2 using SHRiMP to do alignment

  13. ##    http://compbio.cs.toronto.edu/shrimp/README

  14. ##    3.5 Mapping cDNA reads against a miRNA database

  15. cd ~/biosoft/SHRiMP/SHRiMP_2_2_3

  16. export SHRIMP_FOLDER=$PWD

  17. cd -

  18. ##　　We project the database with:

  19. $SHRIMP_FOLDER/utils/project-db.py --seed 00111111001111111100,00111111110011111100,00111111111100111100,00111111111111001100,00111111111111110000 \

  20. --h-flag --shrimp-mode ls miRBase/hairpin.human.fa

  21. ##

  22. $SHRIMP_FOLDER/bin/gmapper-ls -L  hairpin.human-ls SRR1542716.fastq  --qv-offset 33   \

  23. -o 1 -H -E -a -1 -q -30 -g -30 --qv-offset 33 --strata -N 8  >map.out 2>map.log

大家可以看到我们把测序reads比对到前体miRNA和成熟的miRNA结果是有略微区别的，因为一个前体miRNA可以形成多个成熟的miRNA，而并不是所有的成熟的miRNA形式都被记录在数据库，所以一般推荐比对到前体miRNA数据库，这样还可以预测新的成熟miRNA，也是非常有意义的。

另外非常重要的一点是，把U变成T前后比对率差异非常大，这其实是一个非常蠢的错误，我就不多说了。但是做到这一步，其实可以跟文章来做验证，文章有提到比对率，比对的序列。

我也是在博客里面看到这个信息的：

Thank you so much!. Yes I contacted the lab-guy and he just said that trimmed the first 4 bp and last 4bp. ( as you found) So I firstly trimmed the adapter sequences(TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC) And then, trimmed the first 4bp and last 4bp from reads, which leads to the 22bp peak of read-length distribution(instead of 24bp) Anyhow, I tried to map with bowtie2 again.

  1. > bowtie2 --local -N 1 -L 16

  2. > -x ../miRNA_reference/hairpin_UtoT.fa

  3. > -U first4bptrimmed_A1-SmallRNA_S1_L001_R1_001_Illuminaadpatertrim.fastq

  4. > -S f4_trimmed.sam

I also changed hairpin.fa file (U to T) Oh.. thank you David, Finallly, I got

  1. >  2565353 reads; of these:

  2. >  2565353 (100.00%) were unpaired; of these:

  3. >  479292 (18.68%) aligned 0 times

  4. >  11959 (0.47%) aligned exactly 1 time

  5. >  2074102 (80.85%) aligned >1 times

  6. >  81.32% overall alignment rate

第五讲：获取miRNA表达量

得到比对后的sam/bam文件只能算是level2的数据，一般我们给他人share的结果也是直接给表达矩阵的， miRNA分析跟mRNA分析类似，但是它的表达矩阵更好获取一点。

如果是mRNA，我们一般会跟基因组来比较，而基因组就是24条参考染色体，想知道具体比对到了哪个基因，需要根据基因组注释文件来写程序提取表达量信息，现在比较流行的是htseq这个软件，我前面也写过教程如何安装和使用，这里就不啰嗦了。但是对于miRNA，因为我比对的就是那1881条前体miRNA序列，所以直接分析比对的sam/bam文件就可以知道每条参考miRNA序列的表达量了。

  1. ## step6: counts the reads which mapping to each miRNA reference.

  2. ## we need to exclude unmapped as well as multiple-mapped  reads

  3. ## XS:i:<n> Alignment score for second-best alignment. Can be negative. Can be greater than 0 in --local mode

  4. ## NM:i:1   ## NM i Edit distance to the reference, including ambiguous bases but excluding clipping

  5. #The following command exclude unmapped (-F 4) as well as multiple-mapped (grep -v “XS:”) reads

  6. #samtools view -F 4 input.bam | grep -v "XS:" | wc -l

  7. ## 180466//1520320

  8. ##cat >count.hairpin.sh

  9. ls *hairpin.sam  | while read id

  10. do

  11. samtools view  -SF 4 $id |perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h }'  > ${id%%_*}.hairpin.counts

  12. done

  13. ## bash count.hairpin.sh

  14. ##cat >count.mature.sh

  15. ls *mature.sam  | while read id

  16. do

  17. samtools view  -SF 4 $id |perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h }'  > ${id%%_*}.mature.counts

  18. done

  19. ## bash count.mature.sh

上面的代码，是我自己写的脚本来算表达量，非常简单，因为我没有考虑细节，直接想得到各个样本测序数据的表达量而已。如果是比对到了参考基因组，就要根据miRNA的gff注释文件用htseq等软件来计算表达量。

得到了表达量，就可以跟文献来做比较

  1. ### step7: compare the results with paper's

  2. GSM1470353: control-CM, experiment1; Homo sapiens; miRNA-Seq   SRR1542714

  3. GSM1470354: ET1-CM, experiment1; Homo sapiens; miRNA-Seq  SRR1542715

  4. GSM1470355: control-CM, experiment2; Homo sapiens; miRNA-SeqSRR1542716

  5. GSM1470356: ET1-CM, experiment2; Homo sapiens; miRNA-Seq SRR1542717

  6. GSM1470357: control-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542718

  7. GSM1470358: ET1-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542719

  8. ### 下面我用R语言来检验一下，我得到的分析结果跟文章发表的结果的区别。

  9. a=read.table("bowtie_bam/SRR1542714.mature.counts")

  10. b=read.table("paper_results/GSM1470353_iPS_010313_Unstim_known_miRNA_counts.txt")

  11. plot(log(tmp[,2]),log(tmp[,3]))

  12. cor(tmp[,2],tmp[,3])

  13. ##[1] 0.8413439

表达量相关性还不错，至少说明我跟作者的分析结果非常一致，应该是没有分析错。

这个代码是我自己根据文章的理解写出的，因为我本身不擅长miRNA数据分析，所以在进行alignment的时候参数选择可能并不是那么友好。

第六讲：miRNA表达量差异分析

这一讲是miRNA-seq数据分析的分水岭，前面的5讲说的是读文献下载数据比对，然后计算表达量，属于常规的流程分析，一般在公司测序之后都可以拿到分析结果，或者文献也会给出下载结果。

但是单纯的分析一个样本意义不大，一般来说，我们做研究都是针对于不同状态下的miRNA表达量差异分析，然后做注释，功能分析，网络分析，这才是重点和难点。

我这里就直接拿文献处理好的miRNA表达量来展示如何做下游分析，首先就是差异分析。

根据文献，我们可以知道样本的分类情况是

GSM1470353: control-CM, experiment1; Homo sapiens; miRNA-Seq SRR1542714

GSM1470354: ET1-CM, experiment1; Homo sapiens; miRNA-Seq SRR1542715

GSM1470355: control-CM, experiment2; Homo sapiens; miRNA-SeqSRR1542716

GSM1470356: ET1-CM, experiment2; Homo sapiens; miRNA-Seq SRR1542717

GSM1470357: control-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542718

GSM1470358: ET1-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542719

可以看到是6个样本的测序数据，分成两组，就是ET1刺激了CM细胞系前后对比而已！

同时，我们也拿到了这6个样本的表达矩阵，计量单位是counts的reads数，所以我们一般会选用DESeq2，edgeR这样的常用包来做差异分析，当然，做差异分析的工具还有十几个，我这里只是拿一根最顺手的举例子，就是DESeq2。

下面的代码有点长，因为我在bioconductor系列教程里面多次提到了DESeq2使用方法，这里就只贴出代码，反正我要说的重点是，我们通过差异分析得到了差异miRNA列表

  1. ### step8: differential expression analysis by R package for miRNA expression patterns:

  2. ## 文章里面提到的结果是：

  3. MicroRNA sequencing revealed over 250 known and 34 predicted novel miRNAs to be differentially expressed between ET-1 stimulated and unstimulated control hiPSC-CMs.

  4. ## (FDR < 0.1 and 1.5 fold change)

  5. rm(list=ls())

  6. setwd('J:\\miRNA_test\\paper_results')  ##把从GEO里面下载的文献结果放在这里

  7. sampleIDs=c()

  8. groupList=c()

  9. allFiles=list.files(pattern = '.txt')

  10. i=allFiles[1]

  11. sampleID=strsplit(i,"_")[[1]][1]

  12. treat=strsplit(i,"_")[[1]][4]

  13. dat=read.table(i,stringsAsFactors = F)

  14. colnames(dat)=c('miRNA',sampleID)

  15. groupList=c(groupList,treat)

  16. for (i in allFiles[-1]){

  17. sampleID=strsplit(i,"_")[[1]][1]

  18. treat=strsplit(i,"_")[[1]][4]

  19. a=read.table(i,stringsAsFactors = F)

  20. colnames(a)=c('miRNA',sampleID)

  21. dat=merge(dat,a,by='miRNA')

  22. groupList=c(groupList,treat)

  23. }

  24. ### 上面的代码只是为了把6个独立的表达文件给合并成一个表达矩阵

  25. ## we need to filter the low expression level miRNA

  26. exprSet=dat[,-1]

  27. rownames(exprSet)=dat[,1]

  28. suppressMessages(library(DESeq2))

  29. exprSet=ceiling(exprSet)

  30. (colData <- data.frame(row.names=colnames(exprSet), groupList=groupList))

  31. ## DESeq2就是这么简单的用

  32. dds <- DESeqDataSetFromMatrix(countData = exprSet,

  33. colData = colData,

  34. design = ~ groupList)

  35. dds <- DESeq(dds)

  36. png("qc_dispersions.png", 1000, 1000, pointsize=20)

  37. plotDispEsts(dds, main="Dispersion plot")

  38. dev.off()

  39. res <- results(dds)

  40. ## 画一些图，相当于做QC吧

  41. png("RAWvsNORM.png")

  42. rld <- rlogTransformation(dds)

  43. exprSet_new=assay(rld)

  44. par(cex = 0.7)

  45. n.sample=ncol(exprSet)

  46. if(n.sample>40) par(cex = 0.5)

  47. cols <- rainbow(n.sample*1.2)

  48. par(mfrow=c(2,2))

  49. boxplot(exprSet,  col = cols,main="expression value",las=2)

  50. boxplot(exprSet_new, col = cols,main="expression value",las=2)

  51. hist(exprSet[,1])

  52. hist(exprSet_new[,1])

  53. dev.off()library(RColorBrewer)

  54. (mycols <- brewer.pal(8, "Dark2")[1:length(unique(groupList))])

  55. # Sample distance heatmap

  56. sampleDists <- as.matrix(dist(t(exprSet_new)))

  57. #install.packages("gplots",repos = "http://cran.us.r-project.org")

  58. library(gplots)

  59. png("qc-heatmap-samples.png", w=1000, h=1000, pointsize=20)

  60. heatmap.2(as.matrix(sampleDists), key=F, trace="none",

  61. col=colorpanel(100, "black", "white"),

  62. ColSideColors=mycols[groupList], RowSideColors=mycols[groupList],

  63. margin=c(10, 10), main="Sample Distance Matrix")

  64. dev.off()

  65. png("MA.png")

  66. DESeq2::plotMA(res, main="DESeq2", ylim=c(-2,2))

  67. dev.off()

  68. ## 重点就是这里啦，得到了差异分析的结果

  69. resOrdered <- res[order(res$padj),]

  70. resOrdered=as.data.frame(resOrdered)

  71. write.csv(resOrdered,"deseq2.results.csv",quote = F)

  72. ##下面也是一些图，主要是看看样本之间的差异情况

  73. library(limma)

  74. plotMDS(log(counts(dds, normalized=TRUE) + 1))

  75. plotMDS(log(counts(dds, normalized=TRUE) + 1) - log(t( t(assays(dds)[["mu"]]) / sizeFactors(dds) ) + 1))

  76. plotMDS( assays(dds)[["counts"]] )  ## raw count

  77. plotMDS( assays(dds)[["mu"]] ) ##- fitted values.

最后我们得到的差异分析结果：deseq2.results.csv，就可以根据FDR和fold change来挑选符合要求的差异miRNA。

第七讲：miRNA样本配对mRNA表达量获取

这一讲其实算不上自学miRNA-seq分析，本质是affymetrix的mRNA表达芯片数据分析，而且还是最常用的那种GPL570 HG-U133Plus2，但因为是跟miRNA样本配对检测而且后面会利用到这两个数据分析结果来做共表达网络分析等等，所以就贴出对该芯片数据的分析结果。

文章里面也提到了 Messenger RNA expression analysis identified 731 probe sets with significant differential expression，作者挑选的差异分析结果的显著基因列表如下 http://journals.plos.org/plosone/article/asset?unique&id=info:doi/10.1371/journal.pone.0108051.s002 mRNA expression array - GSE60291 (Affymetrix Human Genome U133 Plus 2.0 Array)

hgu133plus2 芯片数据很常见，可以从GEO里面下载该study的原始测序数据，然后用affy,limma包来分析，也可以直接用GEOquery包来下载作者分析好的表达矩阵，然后直接做差异分析。我这里选择的是后者，而且我跟作者分析方法有一点区别是，我先把探针都注释好了基因，然后只挑最大表达量的基因。而作者是直接对探针为单位的的表达矩阵进行差异分析，对分析结果里面的探针进行基因注释。我这里无法给出哪种方法好的绝对评价。代码如下

  1. m(list=ls())

  2. library(GEOquery)

  3. library(limma)

  4. GSE60291 <- getGEO('GSE60291', destdir=".",getGPL = F)

  5. #下面是表达矩阵

  6. exprSet=exprs(GSE60291[[1]])

  7. library("annotate")

  8. GSE60291[[1]]

  9. ## 下面是分组信息

  10. pdata=pData(GSE60291[[1]])

  11. treatment=factor(unlist(lapply(pdata$title,function(x) strsplit(as.character(x),"-")[[1]][1])))

  12. #treatment=relevel(treatment,'control')

  13. ## 下面做基因注释

  14. platformDB='hgu133plus2.db'

  15. library(platformDB, character.only=TRUE)

  16. probeset <- featureNames(GSE60291[[1]])

  17. #EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

  18. SYMBOL <-  lookUp(probeset, platformDB, "SYMBOL")

  19. ## 下面对每个基因挑选最大表达量探针

  20. a=cbind(SYMBOL,exprSet)

  21. ## remove the duplicated probeset

  22. rmDupID <-function(a=matrix(c(1,1:5,2,2:6,2,3:7),ncol=6)){

  23. exprSet=a[,-1]

  24. rowMeans=apply(exprSet,1,function(x) mean(as.numeric(x),na.rm=T))

  25. a=a[order(rowMeans,decreasing=T),]

  26. exprSet=a[!duplicated(a[,1]),]

  27. #

  28. exprSet=exprSet[!is.na(exprSet[,1]),]

  29. rownames(exprSet)=exprSet[,1]

  30. exprSet=exprSet[,-1]

  31. return(exprSet)

  32. }

  33. exprSet=rmDupID(a)

  34. rn=rownames(exprSet)

  35. exprSet=apply(exprSet,2,as.numeric)

  36. rownames(exprSet)=rn

  37. exprSet[1:4,1:4]

  38. #exprSet=log(exprSet) ## based on e

  39. boxplot(exprSet,las=2)

  40. ## 下面用limma包来进行芯片数据差异分析

  41. design=model.matrix(~ treatment)

  42. fit=lmFit(exprSet,design)

  43. fit=eBayes(fit)

  44. #vennDiagram(decideTests(fit))

  45. DEG=topTable(fit,coef=2,n=Inf,adjust='BH')

  46. dim(DEG[abs(DEG[,1])>1.2 & DEG[,5]<0.05,])  ## 806 genes

  47. write.csv(DEG,"ET1-normal.DEG.csv")

得到的ET1-normal.DEG.csv 文件就是我们的差异分析结果，可以跟文章提供的差异结果做比较，几乎一模一样。

如果根据logFC:1.2 和pValue:0.05来挑选，可以拿到806个基因。

第八讲：miRNA-mRNA表达相关下游分析

通过前面的分析，我们已经量化了ET1刺激前后的细胞的miRNA和mRNA表达水平，也通过成熟的统计学分析分别得到了差异miRNA和mRNA，这时候我们就需要换一个参考文献了，因为前面提到的那篇文章分析的不够细致，我这里选择了浙江大学的一篇TCGA数据挖掘分析文章：

Identifying miRNA/mRNA negative regulation pairs in colorectal cancer

里面首先查找miRNA-mRNA基因对，因为miRNA主要还是负向调控mRNA表达，所以根据我们得到的两个表达矩阵做相关性分析，很容易得到符合统计学意义的miRNA-mRNA基因对，具体分析内容如下

  • 把得到的差异miRNA的表达量画一个热图，看看它是否能显著的分类

  • 用miRWalk2.0等数据库或者根据来获取这些差异miRNA的validated target genes

  • 然后看看这些pairs of miRNA- target genes的表达量相关系数，选取显著正相关或者负相关的pairs

  • 这些被选取的pairs of miRNA- target genes拿去做富集分析

  • 最后这些pairs of miRNA- target genes做PPI网络分析

首先我们看第一个热图的实现

  1. resOrdered=na.omit(resOrdered)

  2. DEmiRNA=resOrdered[abs(resOrdered$log2FoldChange)>log2(1.5) & resOrdered$padj <0.01 ,]

  3. write.csv(resOrdered,"deseq2.results.csv",quote = F)

  4. DEmiRNAexprSet=exprSet[rownames(DEmiRNA),]

  5. write.csv(DEmiRNAexprSet,'DEmiRNAexprSet.csv')

  6. DEmiRNAexprSet=read.csv('DEmiRNAexprSet.csv',stringsAsFactors = F)

  7. exprSet=as.matrix(DEmiRNAexprSet[,2:7])

  8. rownames(exprSet)=rownames(DEmiRNAexprSet)

  9. heatmap(exprSet)

  10. gplots::heatmap.2(exprSet)

  11. library(pheatmap)

  12. ## http://biit.cs.ut.ee/clustvis/

因为我前面保存的表达量就基于counts的，所以画热图还需要进行normalization，我这里懒得弄了，就用了一个网页版工具可以自动生成热图：http://biit.cs.ut.ee/clustvis/(并不是懒，其实就是想推荐这个工具而已)

miRNA-heatmap

感觉还不错，可以很清楚的看到ET1刺激前后细胞中miRNA表达量变化

然后就是检验我们感兴趣的有显著差异的miRNA的target genes，这时候有两种方法：一个是先由数据库得到已经被检验的miRNA的target genes；另一种是根据miRNA和mRNA表达量的相关性来预测。

用数据库来查找MiRNA的作用基因，非常多的工具，比较常用的有TargetScan/miRTarBase

  • http://nar.oxfordjournals.org/content/early/2015/11/19/nar.gkv1258.full

  • http://mirtarbase.mbc.nctu.edu.tw/

  • http://mirtarbase.mbc.nctu.edu.tw/cache/download/6.1/hsa_MTI.xlsx

  • http://www.targetscan.org/vert_71/ (version 7.1 (June 2016))

我还看到过一个整合工具： miRecords (DIANA-microT, MicroInspector, miRanda, MirTarget2, miTarget, NBmiRTar, PicTar, PITA, RNA22, RNAhybrid and TargetScan/TargertScanS)里面提到了查找miRNA的作用基因这一过程，高假阳性，至少被5种工具支持，才算是真的。

还有很多类似的工具，miRWalk2，psRNATarget 网页版工具。

最后值得一提的是中山大学的： starBase

Pan-Cancer Analysis Platform is designed for deciphering Pan-Cancer Networks of lncRNAs, miRNAs, ceRNAs and RNA-binding proteins (RBPs) by mining clinical and expression profiles of 14 cancer types (>6000 samples) from The Cancer Genome Atlas (TCGA) Data Portal (all data available without limitations).

虽然我没有仔细用，但是看介绍好牛的样子。

还有一个R包：miRLAB，它是先通过算所有配对的miRNA- genes的表达量相关系数，选取显著正相关或者负相关的pairs，然后反过来通过已知数据库来验证。

后面我就不讲了，主要看你得到miRNA的时候其它生物学数据是否充分，如果是癌症病人，有生存相关数据，可以做生存分析，如果你同时测了甲基化数据，可以做甲基化相关分析。

如果只是单纯的miRNA测序数据，可以回过头去研究一下de novo的miRNA预测的步骤，也是研究重点。

编辑校对：思考问题的熊

赞 (0)

相关推荐