友情提示：本文共有 10744 个字，阅读大概需要 22 分钟。
microRNAs靶基因数据库哪家强
使用miRNAtap数据源提取miRNA的预测靶基因结果
对miRNA进行go和kegg等功能数据库数据库注释
很多粉丝留言想听miRNA-seq数据分析流程，主要是上游分析流程，因为下游其实就是表达矩阵分析。跟RNA-seq拿到的counts矩阵是类似的分析策略，只不过是miRNA-seq热度已经过去了，我也仅仅是五年前接触过一次。其实miRNA-seq数据上游分析有两个方案，一个是仅仅是针对已知的miRNA进行定量，这样的话无需比对到物种参考基因组，仅仅是比对到miRNA序列合集即可。另外一个挖掘新的miRNA，主要是考虑到人类对miRNA的认知的不停的进步。当然，如果你想 系统性学习，可以参考我五年前在生信菜鸟团自学miRNA-seq分析系列笔记：自学miRNA-seq分析第八讲~miRNA-mRNA表达相关下游分析 | 生信菜鸟团
自学miRNA-seq分析第七讲~miRNA样本配对mRNA表达量获取 | 生信菜鸟团
自学miRNA-seq分析第六讲~miRNA表达量差异分析 | 生信菜鸟团
自学miRNA-seq分析第五讲~miRNA表达量获取 | 生信菜鸟团
自学miRNA-seq分析第四讲~测序数据比对 | 生信菜鸟团
自学miRNA-seq分析第三讲~公共测序数据下载 | 生信菜鸟团
自学miRNA-seq分析第二讲~学习资料的搜集 | 生信菜鸟团
自学miRNA-seq分析第一讲~文献选择与解读 | 生信菜鸟团
当然了，如果你是微信阅读，也可以点击：一篇文章学会miRNA-seq分析首先下载miRBase的miRNA参考序列文件
miRBase是miRNA研究领域最权威的数据库，提供miRNAs序列以及注释查询、定位、发卡序列查询，以及提供命名服务。截止到现在（2020-04-28 ），miRBase 22.1 共收录了271个物种的总共38589 条miRNA信息。其中人类的共收录了1917条pre-miRNA(前体)，以及2656条成熟的miRNAs。见：http://www.mirbase.org/前体miRNA和成熟的miRNA的区别，前体miRNA长度一般是70–120 碱基，一般是茎环结果，也就是发夹结构，所以叫做hairpin。成熟之后，一般22 个碱基。在 http://www.mirbase.org/ 对应的ftp网站下载如下所示文件：mkdir-p~/reference/miRNA
cd~/reference/miRNA
wgetftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.gz##38589reads
wgetftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.zip##48885reads
wgetftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.zip
wgetftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3
unziphairpin.fa.zip
unzipmature.fa.zip
grepsapiensmature.fa|wc-l#2656
grepsapienshairpin.fa|wc#1917
##Homosapiens
perl-alne'{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};nextif$tmp!=1;s/U/T/gif!/>/;print}'hairpin.fa>hairpin.human.fa
perl-alne'{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};nextif$tmp!=1;s/U/T/gif!/>/;print}'mature.fa>mature.human.fa
#这里值得一提的是miRBase数据库下载的序列，居然都是用U表示的，也就是说就是miRNA序列，而不是转录成该miRNA的基因序列，而我们测序的都是基因序列。
最后得到的文件如下：5.9MMar122018hairpin.fa
1.5MApr2915:09hairpin.fa.gz
1.5MApr2915:11hairpin.fa.zip
263KApr2915:13hairpin.human.fa
523KApr2915:11hsa.gff3
3.7MMar122018mature.fa
786KApr2915:09mature.fa.zip
196KApr2915:13mature.human.fa
构建miRNA-seq数据分析环境只需要使用上面的的文件即可。然后使用conda配置好软件环境
仍然是参考我五年前整理的流程，使用bowtie软件和SHRiMP软件进行比对，然后fastx_toolkit 和fastqc进行质量控制。condacreate-nmiRNA
condaactivatemiRNA
condainstall-y-cbiocondahisat2bowtiesamtoolsfastpbowtie2fastx_toolkitfastqc
#sra-toolkits,subreads也可以一并下载
#condasearchfastx_toolkit
#耗费约1.2G的磁盘空间
因为SHRiMP软件太多年没有更新，所以放弃它，反正比对这个过程，有bowtie就ok了。针对miRBase数据库下载的序列构建bowtie索引
需要注意的是bowtie和bowtie2不一样哦：http://bowtie-bio.sourceforge.net/index.shtml
http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
bowtie-buildmature.human.fahsa-mature-bowtie-index
bowtie-buildhairpin.human.fahsa-hairpin-bowtie-index
得到的文件如下：4.2MApr2915:42hsa-hairpin-bowtie-index.1.ebwt
20KApr2915:42hsa-hairpin-bowtie-index.2.ebwt
17KApr2915:42hsa-hairpin-bowtie-index.3.ebwt
39KApr2915:42hsa-hairpin-bowtie-index.4.ebwt
4.2MApr2915:42hsa-hairpin-bowtie-index.rev.1.ebwt
20KApr2915:42hsa-hairpin-bowtie-index.rev.2.ebwt
4.2MApr2915:41hsa-mature-bowtie-index.1.ebwt
7.1KApr2915:41hsa-mature-bowtie-index.2.ebwt
24KApr2915:41hsa-mature-bowtie-index.3.ebwt
15KApr2915:41hsa-mature-bowtie-index.4.ebwt
4.2MApr2915:41hsa-mature-bowtie-index.rev.1.ebwt
7.1KApr2915:41hsa-mature-bowtie-index.rev.2.ebwt
下载测试数据
这里我们仍然是使用 RNA expression profiling of human iPSC-derived cardiomyocytes in a cardiac hypertrophy model. PLoS One 2014;9(9):e108051. PMID: 25255322 文章里面的 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60292 数据集，如下：GSM1470353control-CM,experiment1
GSM1470354ET1-CM,experiment1
GSM1470355control-CM,experiment2
GSM1470356ET1-CM,experiment2
GSM1470357control-CM,experiment3
GSM1470358ET1-CM,experiment3
对应的SRR1542714
SRR1542715
SRR1542716
SRR1542717
SRR1542718
SRR1542719
仍然是参考：一篇文章学会miRNA-seq分析，简单看看规律ftp://ftp.sra.ebi.ac.uk/vol1/srr/SRR154/004/SRR1542714；ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR154/004/SRR1542714/SRR1542714.fastq.gz
#从14到19
#可以对这6个样本，分开prefetch
~/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/prefetchSRR1542714
脚本如下：mkdir-p~/reference/miRNA
cd~/reference/miRNA
#step1:downloadrawdata
mkdirmiRNA_test&&cdmiRNA_test
echo{14..19}|sed's//n/g'|whilereadid;
do(~/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/prefetchSRR15427${id});
done
#下面是另外一个课题，可以参考比较
echo{44..59}|sed's//n/g'|whilereadid;
do(~/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/prefetchSRR77077${id});
done
#另外，这样的下载方式，在中国大陆会非常慢！！！
#建议换成aspera哈
得到的sra文件如下：33MApr2916:07SRR1542714.sra
31MApr2916:07SRR1542715.sra
50MApr2916:07SRR1542716.sra
24MApr2916:07SRR1542717.sra
52MApr2916:07SRR1542718.sra
34MApr2916:07SRR1542719.sra
然后批量转为fq文件， 代码如下：dump=/home/yb77613/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump
echo{14..19}|sed's//n/g'|whilereadid;
do($dump-O./--gzip--split-3SRR15427${id});
done
文件如下：50MApr2916:14SRR1542714.fastq.gz
47MApr2916:15SRR1542715.fastq.gz
75MApr2916:15SRR1542716.fastq.gz
35MApr2916:15SRR1542717.fastq.gz
81MApr2916:16SRR1542718.fastq.gz
50MApr2916:16SRR1542719.fastq.gz
质控和清洗测序数据
清洗前后，都走一下fastqc图表，清洗主要是fastx_toolkit修剪，代码如下：ls*gz|whilereadid
do
echo$id
#fastqc$id
zcat$id|fastq_quality_filter-v-q20-p80-Q33-i--otmp;
fastx_trimmer-v-f1-l27-m15-itmp-Q33-z-o${id%%.*}_clean.fq.gz;
#fastqc${id%%.*}_clean.fq.gz
done
fastq_quality_filter和fastx_trimmer两个命令，都是来自于fastx_toolkit软件包：fastq_quality_filter 代表根据质量过滤序列
fastx_trimmer 代表缩短FASTQ或FASTQ文件中的读数，根据质量修剪（剪切）序列。
其参数详解如下：#运行一次查看是否可以成功运行
#fastq_quality_filter代表根据质量过滤序列
#-v,即verbose,报告序列的数目
#-q,即quality,保留碱基所要求的最小质量评分，低于此值将会去除
#-p,即percent,即序列中超过最小质量评分的碱基数目的最小百分率，低于这个比率将删除
#-Q33,即phred33,默认是按照phred64作为参考的
#-i,即infile,代表输入文件，注意不能是压缩文件，可以是FASTA/FASTQ都行。默认是STDIN，即标准输入
#-o,即outfile,代表输出文件（注意不是outputdir即输出目录，此处输出是一个文件而不是文件夹），默认是STDOUT即标准输出，指定则输出到指定文件
fastq_quality_filter-v-q20-p80-Q33-iSRR1542714.1.fastq-otmp#生成第一步处理的临时文件
#中间文件是 tmp ，到时候可以删除掉。
#fastx_trimmer 代表缩短FASTQ或FASTQ文件中的读数，根据质量修剪（剪切）序列。
#-v,即verbose,报告序列的数目
#-f,即first,代表保留第几个碱基，默认是保留第一个
#-l,即last,代表保留最后的碱基，默认是整个reads
#-i,即infile,代表输入文件，注意不能是压缩文件，可以是FASTA/FASTQ都行。默认是STDIN，即标准输入
#-z,即compress,代表的是使用Gzip压缩输出
#-o,即outfile,代表输出文件（注意不是outputdir即输出目录，此处输出是一个文件而不是文件夹），默认是STDOUT即标准输出，指定则输出到指定文件
fastx_trimmer-v-f1-l27-itmp-Q33-z-oSRR1542714.1_clean.fq.gz#生成最后的文件
所以 fastq_quality_filter 和 fastx_trimmer命令是有替代品的，就是需要去自行搜索，如果你是要建立好用的miRNA-seq数据分析环境，就需要自己耗费大量时间在里面哦。最后得到的干净的测序数据文件如下：39MApr2916:28SRR1542714_clean.fq.gz
37MApr2916:28SRR1542715_clean.fq.gz
49MApr2916:29SRR1542716_clean.fq.gz
18MApr2916:29SRR1542717_clean.fq.gz
68MApr2916:30SRR1542718_clean.fq.gz
30MApr2916:30SRR1542719_clean.fq.gz
前面的步骤，需要自行研究里面的细节哦。比对 miRBase数据库下载的序列 +定量
mature=/home/yb77613/reference/miRNA/hsa-mature-bowtie-index
hairpin=/home/yb77613/reference/miRNA/hsa-hairpin-bowtie-index
ls*_clean.fq.gz|whilereadid
do
echo$id
bowtie-n0-m1--best--strata$mature$id-S${id}_matrue.sam
bowtie-n0-m1--best--strata$hairpin$id-S${id}_hairpin.sam
done
ls*.sam|whilereadid;do(samtoolssort-Obam-@5-o$(basename${id}'.sam').bam${id});done
rm*.sam
使用miRNA序列比对的推荐参数：-n 0 -m1 --best --strata ，理由参考：对其中一个样本的比对的log日志如下：SRR1542714_clean.fq.gz
#readsprocessed:1520320
#readswithatleastonereportedalignment:331645(21.81%)
#readsthatfailedtoalign:1145386(75.34%)
#readswithalignmentssuppresseddueto-m:43289(2.85%)
Reported331645alignments
#readsprocessed:1520320
#readswithatleastonereportedalignment:331482(21.80%)
#readsthatfailedtoalign:1008271(66.32%)
#readswithalignmentssuppresseddueto-m:180567(11.88%)
Reported331482alignments
可以看到，比对率还是蛮低的，但是可以调整参数（允许错配碱基）的数量。得到的bam文件如下：29MApr2920:15SRR1542714_clean.fq.gz_hairpin.bam
28MApr2920:15SRR1542714_clean.fq.gz_matrue.bam
27MApr2920:15SRR1542715_clean.fq.gz_hairpin.bam
26MApr2920:15SRR1542715_clean.fq.gz_matrue.bam
36MApr2920:15SRR1542716_clean.fq.gz_hairpin.bam
35MApr2920:15SRR1542716_clean.fq.gz_matrue.bam
13MApr2920:15SRR1542717_clean.fq.gz_hairpin.bam
13MApr2920:15SRR1542717_clean.fq.gz_matrue.bam
49MApr2920:15SRR1542718_clean.fq.gz_hairpin.bam
48MApr2920:15SRR1542718_clean.fq.gz_matrue.bam
22MApr2920:15SRR1542719_clean.fq.gz_hairpin.bam
22MApr2920:15SRR1542719_clean.fq.gz_matrue.bam
批量走 samtools idxstats 得到表达量矩阵：ls*.bam|xargs-isamtoolsindex{}
ls*.bam|whilereadid;do(samtoolsidxstats${id}>${id}.txt);done
#samtoolsviewmatrue.bam|cut-f3|sort|uniq-c
比对参考基因组+定量
index=/home/yb77613/reference/index/bowtie/hg38
ls*_clean.fq.gz|whilereadid
do
echo$id
bowtie2-p4-x$index-U$id|samtoolssort-@4-o${id}_genome.bam-
done
使用默认参数，对其中一个样本的比对的log日志如下：SRR1542715_clean.fq.gz
1520320reads;ofthese:
1520320(100.00%)wereunpaired;ofthese:
123178(8.10%)aligned0times
333700(21.95%)alignedexactly1time
1063442(69.95%)aligned>1times
91.90%overallalignmentrate
可以看到，这个时候的比对率不得了，得到的bam文件如下：34MApr2920:30SRR1542714_clean.fq.gz_genome.bam
31MApr2920:30SRR1542715_clean.fq.gz_genome.bam
42MApr2920:31SRR1542716_clean.fq.gz_genome.bam
15MApr2920:31SRR1542717_clean.fq.gz_genome.bam
57MApr2920:32SRR1542718_clean.fq.gz_genome.bam
25MApr2920:32SRR1542719_clean.fq.gz_genome.bam
然后定量，需要使用下面的命令和参数，都是需要看软件文档摸索的。gtf=/home/yb77613/reference/miRNA/hsa.gff3
bin_featureCounts='/home/yb77613/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts';
$bin_featureCounts-T4-Fgff-M-tmiRNA-gName-a$gtf-oall.counts.mature.txt*genome*1>counts.mature.log2>&1
$bin_featureCounts-T4-Fgff-M-tmiRNA_primary_transcript-gName-a$gtf-oall.counts.hairpin.txt*genome*1>counts.hairpin.log2>&1
这样我们就有不同的定量流程拿到的不同的表达矩阵啦。比较两个定量方式的表达矩阵差异
在featureCounts软件的结果文件 all.counts.id.txt 里面的信息如下：hsa-miR-12136chr1632668632685-1826215067819
hsa-miR-34a-5pchr191517359151756-22300994943467413652997097
hsa-miR-34a-3pchr191516939151714-228815310086126132
查询其中一个，发现featureCounts软件定量少了80%，比如hsa-miR-12136在featureCounts流程，是 18 26 21 50 6 78 19，但是在miRBase数据库流程，都是五倍以上的表达量。SRR1542714_clean.fq.gz_matrue.bam.txt:hsa-miR-12136181040
SRR1542715_clean.fq.gz_matrue.bam.txt:hsa-miR-12136181820
SRR1542716_clean.fq.gz_matrue.bam.txt:hsa-miR-12136181090
SRR1542717_clean.fq.gz_matrue.bam.txt:hsa-miR-1213618440
SRR1542718_clean.fq.gz_matrue.bam.txt:hsa-miR-12136182350
SRR1542719_clean.fq.gz_matrue.bam.txt:hsa-miR-1213618920
查询另外2个，有意思的是其中一个发现featureCounts软件定量多了1倍，另外一个多了2倍。SRR1542714_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2217450
SRR1542715_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2258250
SRR1542716_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2220150
SRR1542717_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2226050
SRR1542718_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2230810
SRR1542719_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-5p2244130
SRR1542714_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22250
SRR1542715_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22690
SRR1542716_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22210
SRR1542717_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22270
SRR1542718_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22370
SRR1542719_clean.fq.gz_matrue.bam.txt:hsa-miR-34a-3p22570
比如，hsa-miR-34a-5p在featureCounts流程，是3009 9494 3467 4136 5299 7097，但是在miRBase数据库流程，都只有一半的表达量。对hsa-miR-34a-3p来说，在featureCounts流程，是88 153 100 86 126 132，但是在miRBase数据库流程，都是两倍的表达量。这样的定量差异，理论上是不会改变该miRNA的差异分析结果，因为是同样的膨胀系数。但是，这样的不确定性，让我们的miRNA-seq数据分析，蒙上了一层阴影。
收集不易，本文《构建miRNA-seq数据分析环境》知识如果对你有帮助，请点赞收藏并留下你的评论。
