一篇文章学会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月 )的代码如下所示:

wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.gz   ## 28645 readswget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.zip   ##   35828 reads wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.zipwget ftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3 ##wget ftp://mirbase.org/pub/mirbase/CURRENT/miFam.dat.zipgrep sapiens mature.fa |wc   # 2588 grep sapiens hairpin.fa |wc       # 1881 ## Homo sapiensperl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' hairpin.fa >hairpin.human.faperl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' mature.fa >mature.human.fa# 这里值得一提的是miRBase数据库下载的序列,居然都是用U表示的,也就是说就是miRNA序列,而不是转录成该miRNA的基因序列,而我们测序的都是基因序列。

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

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

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

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

博客讲解如何分析miRNA数据

公司数据分析流程:

耶鲁大学

南方基因

miRNA研究整套方案

Biostar 讨论帖子

miRNA-seq数据处理实战指南

直接用一个包搞定

github流程:miRNA Analysis Pipeline v0.2.7

miRNA annotation

网页版分析工具

可视化IGV User Guide

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

第三讲:下载公共数据

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

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

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

## pre-step: download sratoolkit /fastx_toolkit_0.0.13/fastqc/bowtie2/hg19/miRBase/SHRiMP## http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software## http://www.ncbi.nlm.nih.gov/books/NBK158900/## 我这里特意挑选的二进制版本程序下载的,这样直接解压就可以用,但是需要挑选适合自己的操作系统的程序。cd ~/biosoftmkdir sratoolkit &&  cd sratoolkitwget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.6.3/sratoolkit.2.6.3-centos_linux64.tar.gz####  Length: 63453761 (61M) [application/x-gzip]##  Saving to: "sratoolkit.2.6.3-centos_linux64.tar.gz"tar zxvf sratoolkit.2.6.3-centos_linux64.tar.gzcd ~/biosoftmkdir bowtie &&  cd bowtiewget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip/download#Length: 27073243 (26M) [application/octet-stream]#Saving to: "download"mv download  bowtie2-2.2.9-linux-x86_64.zipunzip bowtie2-2.2.9-linux-x86_64.zip## http://compbio.cs.toronto.edu/shrimp/mkdir SHRiMP &&  cd SHRiMPwget http://compbio.cs.toronto.edu/shrimp/releases/SHRiMP_2_2_3.lx26.x86_64.tar.gztar zxvf SHRiMP_2_2_3.lx26.x86_64.tar.gz cd SHRiMP_2_2_3export SHRIMP_FOLDER=$PWD  ## 这个软件使用的时候比较奇葩,需要设置到环境变量,不能简单的调用全路径

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

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

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

## step1 : download raw datamkdir miRNA_test && cd miRNA_testecho {14..19} |sed 's/ /\n/g' |while read id; \do  wget "ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP045/SRP045420/SRR15427$id/SRR15427$id.sra"  ;\done## step2 :  change sra data to fastq files.## 主要是用shell脚本来批量下载ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;donerm *sra##  33M --> 247M#Read 1866654 spots for SRR1542714.sra#Written 1866654 spots for SRR1542714.sra## step3 : download the results from paper## http://www.bio-info-trainee.com/1571.html## ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1009/suppl/GSE1009_RAW.tarmkdir paper_results && cd paper_resultswget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60292/suppl/GSE60292_RAW.tar## tar xvf GSE60292_RAW.tarls *gz |while read id ; do (echo $id;zcat $id | cut -f 2 |perl -alne '{$t+=$_;}END{print $t}');donels *gz |xargs gunzip## step4 : quality assessmentls *fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done## Sequence length 8-109## %GC 52## Adapter Content passed## write a script : :: cat >filter.shls *fastq |while read iddoecho $id~/biosoft/fastx_toolkit_0.0.13/bin/fastq_quality_filter -v -q 20 -p 80 -Q33  -i $id -o tmp ;~/biosoft/fastx_toolkit_0.0.13/bin/fastx_trimmer -v -f 1 -l 27 -i tmp  -Q33 -z -o ${id%%.*}_clean.fq.gz ;donerm tmp## discarded 12%~~49%%ls *_clean.fq.gz | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;donemkdir QC_resultsmv *zip *html QC_results

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

~/biosoft/fastx_toolkit_0.0.13/bin/fastq_quality_filter -v -q 20 -p 80 -Q33  -i $id -o tmp ;~/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数据库下载的人类的参考序列。

## step5 : alignment to miRBase v21 (hairpin.human.fa/mature.human.fa )#### step5.1 using bowtie2 to do alignmentmkdir  bowtie2_index &&  cd bowtie2_index~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ../hairpin.human.fa hairpin_human~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ../mature.human.fa  mature_humanls *_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## overall alignment rate:  10.20% / 5.71%/ 10.18%/ 4.36% / 10.02% / 4.95%  (before convert U to T )## overall alignment rate:  51.77% / 70.38%/51.45% /61.14%/ 52.20% / 65.85% (after convert U to T )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## overall alignment rate:  6.67% / 3.78% / 6.70% / 2.80%/ 6.55% / 3.23%    (before convert U to T )## overall alignment rate:  34.94% / 46.16%/ 35.00%/ 38.50% / 35.46% /42.41%(after convert U to T )#### step5.2 using SHRiMP to do alignment##    http://compbio.cs.toronto.edu/shrimp/README##    3.5 Mapping cDNA reads against a miRNA databasecd ~/biosoft/SHRiMP/SHRiMP_2_2_3export SHRIMP_FOLDER=$PWDcd -##  We project the database with:$SHRIMP_FOLDER/utils/project-db.py --seed 00111111001111111100,00111111110011111100,00111111111100111100,00111111111111001100,00111111111111110000 \ --h-flag --shrimp-mode ls miRBase/hairpin.human.fa##$SHRIMP_FOLDER/bin/gmapper-ls -L  hairpin.human-ls SRR1542716.fastq  --qv-offset 33   \-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.

> bowtie2 --local -N 1 -L 16> -x ../miRNA_reference/hairpin_UtoT.fa> -U first4bptrimmed_A1-SmallRNA_S1_L001_R1_001_Illuminaadpatertrim.fastq> -S f4_trimmed.sam

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

>  2565353 reads; of these:>  2565353 (100.00%) were unpaired; of these:>  479292 (18.68%) aligned 0 times>  11959 (0.47%) aligned exactly 1 time>  2074102 (80.85%) aligned >1 times>  81.32% overall alignment rate

第五讲:获取miRNA表达量

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

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

## step6: counts the reads which mapping to each miRNA reference.## we need to exclude unmapped as well as multiple-mapped  reads## XS:i:<n> Alignment score for second-best alignment. Can be negative. Can be greater than 0 in --local mode## NM:i:1   ## NM i Edit distance to the reference, including ambiguous bases but excluding clipping#The following command exclude unmapped (-F 4) as well as multiple-mapped (grep -v “XS:”) reads#samtools view -F 4 input.bam | grep -v "XS:" | wc -l## 180466//1520320##cat >count.hairpin.shls *hairpin.sam  | while read iddosamtools view  -SF 4 $id |perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h }'  > ${id%%_*}.hairpin.countsdone## bash count.hairpin.sh##cat >count.mature.shls *mature.sam  | while read iddosamtools view  -SF 4 $id |perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h }'  > ${id%%_*}.mature.countsdone## bash count.mature.sh

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

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

### step7: compare the results with paper'sGSM1470353: control-CM, experiment1; Homo sapiens; miRNA-Seq   SRR1542714GSM1470354: ET1-CM, experiment1; Homo sapiens; miRNA-Seq  SRR1542715GSM1470355: control-CM, experiment2; Homo sapiens; miRNA-SeqSRR1542716GSM1470356: ET1-CM, experiment2; Homo sapiens; miRNA-Seq SRR1542717GSM1470357: control-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542718GSM1470358: ET1-CM, experiment3; Homo sapiens; miRNA-Seq SRR1542719### 下面我用R语言来检验一下,我得到的分析结果跟文章发表的结果的区别。a=read.table("bowtie_bam/SRR1542714.mature.counts")b=read.table("paper_results/GSM1470353_iPS_010313_Unstim_known_miRNA_counts.txt")plot(log(tmp[,2]),log(tmp[,3]))cor(tmp[,2],tmp[,3])##[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列表

### step8: differential expression analysis by R package for miRNA expression patterns:## 文章里面提到的结果是:MicroRNA sequencing revealed over 250 known and 34 predicted novel miRNAs to be differentially expressed between ET-1 stimulated and unstimulated control hiPSC-CMs.## (FDR < 0.1 and 1.5 fold change)rm(list=ls())setwd('J:\\miRNA_test\\paper_results')  ##把从GEO里面下载的文献结果放在这里sampleIDs=c()groupList=c()allFiles=list.files(pattern = '.txt')i=allFiles[1]sampleID=strsplit(i,"_")[[1]][1]treat=strsplit(i,"_")[[1]][4]dat=read.table(i,stringsAsFactors = F)colnames(dat)=c('miRNA',sampleID)groupList=c(groupList,treat)for (i in allFiles[-1]){sampleID=strsplit(i,"_")[[1]][1]treat=strsplit(i,"_")[[1]][4]a=read.table(i,stringsAsFactors = F)colnames(a)=c('miRNA',sampleID)dat=merge(dat,a,by='miRNA')groupList=c(groupList,treat)}### 上面的代码只是为了把6个独立的表达文件给合并成一个表达矩阵## we need to filter the low expression level miRNAexprSet=dat[,-1]rownames(exprSet)=dat[,1]suppressMessages(library(DESeq2))exprSet=ceiling(exprSet)(colData <- data.frame(row.names=colnames(exprSet), groupList=groupList))## DESeq2就是这么简单的用dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ groupList)dds <- DESeq(dds)png("qc_dispersions.png", 1000, 1000, pointsize=20)plotDispEsts(dds, main="Dispersion plot")dev.off()res <- results(dds)## 画一些图,相当于做QC吧png("RAWvsNORM.png")rld <- rlogTransformation(dds)exprSet_new=assay(rld)par(cex = 0.7)n.sample=ncol(exprSet)if(n.sample>40) par(cex = 0.5)cols <- rainbow(n.sample*1.2)par(mfrow=c(2,2))boxplot(exprSet,  col = cols,main="expression value",las=2)boxplot(exprSet_new, col = cols,main="expression value",las=2)hist(exprSet[,1])hist(exprSet_new[,1])dev.off()library(RColorBrewer)(mycols <- brewer.pal(8, "Dark2")[1:length(unique(groupList))])# Sample distance heatmapsampleDists <- as.matrix(dist(t(exprSet_new)))#install.packages("gplots",repos = "http://cran.us.r-project.org")library(gplots)png("qc-heatmap-samples.png", w=1000, h=1000, pointsize=20)heatmap.2(as.matrix(sampleDists), key=F, trace="none",col=colorpanel(100, "black", "white"),ColSideColors=mycols[groupList], RowSideColors=mycols[groupList],margin=c(10, 10), main="Sample Distance Matrix")dev.off()png("MA.png")DESeq2::plotMA(res, main="DESeq2", ylim=c(-2,2))dev.off()## 重点就是这里啦,得到了差异分析的结果resOrdered <- res[order(res$padj),]resOrdered=as.data.frame(resOrdered)write.csv(resOrdered,"deseq2.results.csv",quote = F)##下面也是一些图,主要是看看样本之间的差异情况library(limma)plotMDS(log(counts(dds, normalized=TRUE) + 1))plotMDS(log(counts(dds, normalized=TRUE) + 1) - log(t( t(assays(dds)[["mu"]]) / sizeFactors(dds) ) + 1))plotMDS( assays(dds)[["counts"]] )  ## raw countplotMDS( 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包来下载作者分析好的表达矩阵,然后直接做差异分析。我这里选择的是后者,而且我跟作者分析方法有一点区别是,我先把探针都注释好了基因,然后只挑最大表达量的基因。而作者是直接对探针为单位的的表达矩阵进行差异分析,对分析结果里面的探针进行基因注释。我这里无法给出哪种方法好的绝对评价。代码如下

m(list=ls())library(GEOquery)library(limma)GSE60291 <- getGEO('GSE60291', destdir=".",getGPL = F)#下面是表达矩阵exprSet=exprs(GSE60291[[1]])library("annotate")GSE60291[[1]]## 下面是分组信息pdata=pData(GSE60291[[1]])treatment=factor(unlist(lapply(pdata$title,function(x) strsplit(as.character(x),"-")[[1]][1])))#treatment=relevel(treatment,'control')## 下面做基因注释platformDB='hgu133plus2.db'library(platformDB, character.only=TRUE)probeset <- featureNames(GSE60291[[1]])#EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))SYMBOL <-  lookUp(probeset, platformDB, "SYMBOL")## 下面对每个基因挑选最大表达量探针a=cbind(SYMBOL,exprSet)## remove the duplicated probesetrmDupID <-function(a=matrix(c(1,1:5,2,2:6,2,3:7),ncol=6)){exprSet=a[,-1]rowMeans=apply(exprSet,1,function(x) mean(as.numeric(x),na.rm=T))a=a[order(rowMeans,decreasing=T),]exprSet=a[!duplicated(a[,1]),]#exprSet=exprSet[!is.na(exprSet[,1]),]rownames(exprSet)=exprSet[,1]exprSet=exprSet[,-1]return(exprSet)}exprSet=rmDupID(a)rn=rownames(exprSet)exprSet=apply(exprSet,2,as.numeric)rownames(exprSet)=rnexprSet[1:4,1:4]#exprSet=log(exprSet) ## based on eboxplot(exprSet,las=2)## 下面用limma包来进行芯片数据差异分析design=model.matrix(~ treatment)fit=lmFit(exprSet,design)fit=eBayes(fit)#vennDiagram(decideTests(fit))DEG=topTable(fit,coef=2,n=Inf,adjust='BH')dim(DEG[abs(DEG[,1])>1.2 & DEG[,5]<0.05,])  ## 806 geneswrite.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网络分析

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

resOrdered=na.omit(resOrdered)DEmiRNA=resOrdered[abs(resOrdered$log2FoldChange)>log2(1.5) & resOrdered$padj <0.01 ,]write.csv(resOrdered,"deseq2.results.csv",quote = F)DEmiRNAexprSet=exprSet[rownames(DEmiRNA),]write.csv(DEmiRNAexprSet,'DEmiRNAexprSet.csv')DEmiRNAexprSet=read.csv('DEmiRNAexprSet.csv',stringsAsFactors = F)exprSet=as.matrix(DEmiRNAexprSet[,2:7])rownames(exprSet)=rownames(DEmiRNAexprSet)heatmap(exprSet)gplots::heatmap.2(exprSet)library(pheatmap)## 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预测的步骤,也是研究重点。

编辑校对:思考问题的熊

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏PPV课数据科学社区

学习R语言,一篇文章让你从懵圈到入门

在实际工作中,每个数据科学项目各不相同,但基本都遵循一定的通用流程。具体如下: 数据科学工作流程 数据导入 数据整理 反复理解数据 数据可视化 数据转换 ...

3164
来自专栏Bingo的深度学习杂货店

我分析了2837首歌曲,做了个信息检索与信息抽取系统

我把目标锁定在网易云音乐热门的华语男歌手、华语女歌手以及华语组合/乐队,每一类爬取20个热门歌手,这样我就有了60位歌手的信息。

681
来自专栏数据科学与人工智能

【Python环境】python的nltk中文使用和学习资料汇总帮你入门提高

nltk是一个python工具包, 用来处理和自然语言处理相关的东西. 包括分词(tokenize), 词性标注(POS), 文本分类, 等等现成的工具. 1....

3316
来自专栏知晓程序

厉害了!这个小程序,能让你说一口飘准的「普通发」

今天,知晓程序(微信号 zxcx0101)给大家推荐一款「普通话学习评分」小程序,它使用专业的普通话评分系统,你可以录音让它为自己的普通话打分。

964
来自专栏Renderbus云渲染农场

如何基于一张图片来创建3D模型?包含哪些步骤?

首先,基于一张图片(一个角度)来制作一个3D模型是不大现实的。因为三维物体是立体的、多维的,单从一个角度去观察很难判断物体其他视角的构造,制作出来的模型也就一个...

780
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 3

前面两篇文章(高通量数据分析必备|基因组浏览器使用介绍 - 1和高通量数据分析必备|基因组浏览器使用介绍 - 2)介绍了EPGG的基本使用、各部分特征、Trac...

895
来自专栏生信技能树

芯片数据分析,so easy?

我最早接触的高通量数据就是RNA-seq,后来接触的也基本是高通量测序结果而不是芯片数据,因此我从来没有分析过一次芯片数据,而最近有一个学员在看生信技能树在腾讯...

1193
来自专栏专知

【干货】如何为论文设计精致的插图

【导读】论文是经过数月或数年实验得到的数据的积淀。 但是,原始数据或描述本身并不能成为好的期刊文章。 数据可视化工具和免费绘图软件能帮助科学家阐述他们的工作。 ...

1333
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/6/3

1. Judea Pearl上次在NIPS有一张令人唏嘘的照片,不过现在他又回来了,发了新书也给了一个访谈,说深度学习就像是curve fitting(我觉得没...

734
来自专栏CDA数据分析师

学习R语言,一篇文章让你从懵圈到入门

在实际工作中,每个数据科学项目各不相同,但基本都遵循一定的通用流程。具体如下: ? 数据科学工作流程: 1.数据导入 2.数据整理 3.反复理解数据 数据可视...

1906

扫描关注云+社区