1. 使用ShortStack鉴定,注释产生小RNA的基因
1.1 下载安装
$ git clone https://github.com/MikeAxtell/ShortStack.git
$ cd ShortStack/
$ ls
LICENSE README ShortStack
1.2 下载依赖的软件(包)
perl
, samtools
, bowtie
& bowtie-build
这几个都是平时常用的,不用安装了。
另外的RNALfold
& RNAeval
(ViennaRNA
package)以及einverted
(EMBOSS
package)以前没接触过,需要安装。见这里。
那这几个软件有什么作用呢?
RNALfold & RNAeval:Predicts locally stable RNA secondary structures
einverted:Finds inverted repeats in DNA sequences
bowtie & bowtie-build:建立索引,无gap的短reads比对
samtools:处理sam格式文件
1.3 简单使用
ShortStack \
--readfile /srna_project/data/out.SRR1266859.fastq.gz /srna_project/data/out.SRR1266860.fastq.gz \
--genomefile /srna_project/ref/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--bowtie_cores 6 --sort_mem 10G
会自动生成一个以ShortStack_开头的目录,最终生成如下文件:
$ ls -t | cat
Log.txt
ShortStack_All.gff3
ShortStack_D.gff3
ShortStack_N.gff3
Counts.txt
Unplaced.txt
Results.txt
MIRNAs/
merged_alignments.bam
ErrorLogs.txt
中间会生成out.SRR1266860_readsorted.sam.gz(以queryname排序),out.SRR1266860.bam(以坐标排序)等文件,最后会删除;
参考基因组所在的文件夹下还会生成samtools faidx和bowtie-build产生的索引文件。
2. 结果解读
2.1 ErrorLogs.txt
报错信息(如果有的话);还包含一些运行信息,如bowtie-build, bowtie, samtools sort, and samtools merge的详细运行信息。比如下面的就比较重要:
# reads processed: 20394922
# reads with at least one reported alignment: 16290733 (79.88%)
# reads that failed to align: 1705936 (8.36%)
# reads with alignments suppressed due to -m: 2398253 (11.76%)
Reported 60107848 alignments
# reads processed: 27228269
# reads with at least one reported alignment: 22673414 (83.27%)
# reads that failed to align: 2373241 (8.72%)
# reads with alignments suppressed due to -m: 2181614 (8.01%)
Reported 75813835 alignments
2.2 Results.txt
包含主要结果信息,每一列的具体含义见说明文档。
Locus Name Length Reads RPM UniqueReads FracTop Strand MajorRNA MajorRNAReads 'Complexity DicerCall MIRNA PhaseScore' Short Long 20 21 22 23 24
2.3 Unplaced file
显示了每一条满足--mincov的无法定位的小RNA reads。其中第五列表示能比对到参考基因组上的位置数,>=0的任意数;第一列第二列表示序列和长度,如果很小则可能是小RNA的片段(一部分)。
2.4 Counts file
计数文件,可以用来做差异表达分析。
2.5 MIRNAs directory
包含鉴定出来的miRNA基因簇信息,依据的是“Results.txt”中MIRNA这一列为“Y”(yes)的记录会被确定为一个"new miRNA family"。
每一个文件具体信息包括:这段位置的序列,以点和括号表示的二级结构,miRNA and miRNA-star的位置;下方的记录表示其它的小RNA,当比对到相反链时,用"<"表示;l: length of RNA, a: number of alignments.
2.6 gff3文件
根据Results.txt中DicerCalls这一列是否为'N',得到ShortStack_D.gff3,ShortStack_N.gff3,ShortStack_ALL.gff3这几个粗略注释文件
3. ShortStack详细使用指南
Usage: ShortStack [options] {--readfile <r> | {--bamfile <b> | --cramfile <c>}} --genomefile <g>
<r> : *.fasta or *.fa; colorspace-fasta (.csfasta); *.fastq or *.fq; their gzip-compressed versions(多文件时,空格分隔)
<b> : *.bam
<c> : *.cram
<g> : FASTA formatted (.fa or .fasta) genome file
Alignment Options:
--adapter:至少由8个ATGC碱基组成的序列
--bowtie_cores:bowtie比对时所用线程数
--sort_mem:samtools sort所用的内存
--mismatches:比对时能容忍的错配数,默认是1
--mmap:处理多比对reads的方式,n (none), r (random), u (unique-seeded guide), or f (fractional-seeded guide). default: u
'n,忽略多比对reads;r,随机给多比对reads分配一个位置;
u,假设一条READS比对到了基因组的A, B, C三个位置,在这三个位置上划定50nt的窗口,计算窗口内特异性比对reads的数量。READS更有可能落在数量大的那个窗口对应的基因组位置上;
f,在上面比较,不同的是,除了计算窗口内特异性比对reads的数量,还要计算多比对reads的数量;'
—— 配图1
--bowtie_m:多比对最多比对多少次,默认50;超过此值就定为“未比对上”
--ranmax:当--mmap使用u, f时,也可能出现不能准确定位的情况,比如,同时存在多个权重最大的位置。这时,则会随机分配可能位置,但是也有条件,就是可能的位置数不能大于--ranmax,一旦超过就会定为“未比对上”,默认是3
--align_only:跑完比对+定相(处理多比对reads的方法)之后就停止
--show_secondaries:显示primary alignments以及secondary alignments,如果已经比对定相了,这个选项还有用吗?
Analysis Options:
--dicermin:Dicer酶切割得到的小RNA的最小长度,默认20
--dicermax:Dicer酶切割得到的小RNA的最大长度,默认24 ——配图2
--foldsize:在基因组上搜寻MIRNA基因时,考虑到空间结构的折叠(反向重复,茎环结构),片段大小不应该大于一个阈值,一旦大于则不会被看作是MIRNA,默认是300
--locifile:制表符分隔的普通文本文件。第一列Chr:start-stop,第二列为名称(可选)
--locus:与--locifile作用相似,且不能同时使用。Chr:start-stop[,Chr:start-stop,Chr:start-stop......]
--nohp:不执行MIRNA搜索
--pad:如果小RNA的基因簇之间的距离小于等于pad值则会被合并,默认值75 —— 背景知识见下文引用
--mincov:小RNA的基因簇最少的比对reads数,也可以是标准化值,如rpm, rpmm(这时要带小数点),默认值是0.5rpm
--strand_cutoff:对于基因组上一个区段,其实有两条链,需要根据--strand_cutoff这个阈值来定义正负(+-)链。该区段上超过80%的reads比对到参考基因组链,则称为+链;小于等于20%reads比对到参考基因组链定义为-链
关于基因簇:大多数miRNA 基因以单拷贝、多拷贝或基因簇的形式存在于基因组中。基因簇(gene cluster)指基因家族中的各成员紧密成簇排列成大串的重复单位,位于染色体的特殊区域。基因簇少则可以是由重复产生的两个相邻相关基因所组成,多则可以是几百个相同基因串联排列而成。他们属于同一个祖先的基因扩增产物。也有一些基因家族的成员在染色体上排列并不紧密,中间还含有一些无关序列。但总体是分布在染色体上相对集中的区域。——百度词条