ShortStack的安装与运行

2020-04-23 本文已影响0人 guguaihezi

安装

ShortStack依赖的软件包括perl、samtools、bowtie、bowtie-build、gzip、RNAfold，检查了一下，服务器中都符合要求，不需要额外安装或更新了。

版本要求：This release of ShortStack (3.8.4) tested on Mac OSX (10.10.5), perl
5.18.2, samtools 1.3.1, bowtie 1.2.0, RNAfold 2.3.2. It is known that
samtools 1.x and higher is critical (no old 0.x samtools allowed).

pwd
echo 'export PATH=/mnt/data/hezi/biosoft/ShortStack-3.8.5:$PATH' >>~/.bashrc   #注意>>是向文件中追加内容，>是覆盖原有内容
source ~/.bashrc

参数说明

ShortStack [参数] {readfile {bamfile cramfile}} genomefile
#bamfile和cramfile都是比对文件

一般参数
--help
--version
--genomefile [string字符串] 给出基因组数据所在路径
--outdir [string] 给出输出文件所在路径
比对参数
--readfile [string] 需要比对的reads文件路径
--adapter [string] 接头文件，若不指定即默认已经去过接头
--bowtie_cores [integer] 传递给bowtie -p(线程数)的参数。设定处理器核数，默认为1。核与线程的关系：1个核对应2个线程。
--sort_mem [string] 传递给samtools的sort -m(内存数)的参数，即用sort对文件排序时使用的物理内存空间。samtools中默认是768M。
--mismatches [integer] 比对中容忍的错配数，只能为0，1或2。默认为1。如果有多个reads匹配上了，只保留错配数最低的reads。
--cquals [string] 给出color-space类型测序数据的质量值文件。color-space是一类特殊的测序数据，用颜色代表核苷酸，是Solid平台产出数据的原始格式。
--cram 默认输出sam格式，加上此参数可将比对结果输出为cram格式。CRAM和BAM都是压缩性质的比对文件，CRAM在继承了BAM优点的基础上，可对数据进一步压缩。
--mmap [string] 处理multi-mapping reads的流程。可选n (none), r (random), u (unique- seeded guide), or f(fractional-seeded guide). default: u
--bowtie_m [string] max，指定最大比对到基因组的次数。超过50个alignments被标记为unmapped。
--ranmax [string] 当--mmap指定u或f时，仍然可能出现多个权重最大的位置，这时候比对工具还是只能选择随机分配位置，但可能的位置数不能超过ranmax值。默认为3。
--align_only 出现此参数时，ShortStack只会执行到比对这一步，不再进行后续分析。
--show_secondaries 出现此参数时，比对文件中除了有primary alignments，还有secondary alignments信息。
--keep_quals 出现此参数时，会保留比对质量值信息。
分析参数
--bamfile sRNA的比对文件，和readfile、cramfile选一个作为输入文件就行。
--cramfile 同上。
--dicermin Dicer酶切割后的sRNA的最小size。一般大于15，默认为20.
--dicermax Dicer酶切割后的sRNA的最大size。一般大于15，默认为24.
--foldsize 提取参考基因组序列进行RNA二级结构预测的长度。范围：200-1000.默认为300.增加数值后运行时间将大大增加，但也可能找到更多MIRNAs.这里指的应该是miRNA的前体，一般认为miRNA的前体长度不超过300nt。
--locifile 包含特定的intervals信息的文本文件。和--locus互斥。
loci文件示例
--locus 同上
--nohp 禁用miRNA search。
--pad 如果sRNA clusters之间的距离小于等于pad值，sRNA clusters就会被合并。默认值：75.
--mincov sRNA clusters表达量的最低值。默认值：0.5rpm
--strand_cutoff strand的临界值。默认为0.8 At default of 0.8, a locus must have 80% of more of its reads on the top strand to be called a + strand locus, or 20% or less on the top strand to be a - strand locus. All others receive no strand call(e.g. '.').+strand指的是参考基因组中的那条链，-strand也就是与其互补的另一条链。注意与正义链、反义链的区分。

--total_primaries 直接告诉ShortStack bam文件中的primary alignments的数目，减少计算时间，包括没有比对上的。

作者最后给出了一些运行建议：比如bowtie -m和 --ranmax参数设置对文件size影响很大。

比对方法说明

ShortStack可以把比对过的文件作为输入文件，也可以使用其内置的比对方法：Improved Placement of Multi-mapping Small RNAs（doi:10.1534/g3.116.030452）
主要解决multi-mapping reads 如何放到合适位置的问题。

拓展01：multi-mapping reads在sRNA-seq里很常见，一是因为sRNA reads很短，二是sRNA倾向于从多拷贝区域起源。以前的序列比对工具要么准确性低（随机选取一个匹配分数高的位置），要么敏感性低（搞不清楚干脆就不要了，忽略掉所有multi-mapping reads），bowtie就是代表，它默认随机选择，也可以设置主动忽略。
拓展02：比对方法的基本原理——1)用bowtie找出每条read的best-scoring alignments，每条read最多允许有50个alignments；2）采用local-weighting为每个alignment计算概率；3）将概率作为权重，权重大的作为primary alignments，其余的作为secondary alignments；4)比对模式选择，讲了这么多就是告诉你选U最好。
拓展03：这里比对指的是比对到sRNA转录起源的地方，而非target sites。靶位点的预测一般会有一些species-specific parameters以及只考虑基因组的部分区域（如成熟转录区）。

ShortStack比对方法说明

基因组预处理

基因组文件需为fasta格式；
染色体名最好不要有空格，否则只保留空格前面的部分；
如果基因组拥有>50chromosomes/scaffolds/contigs，且N50<1Mb，ShortStack会默认拼接短的染色体。但在3.8之后的版本中，报告的搜索结果仍然是基于原始基因组的。

Reads预处理

sRNA文件可以为压缩格式；
同时将多个文件作为输入文件，之间用逗号隔开；
不支持双端文件；
No condensation?
所有reads必须有唯一名；
ShortStack也可以通过指定--adapter去接头，不过可能比较简陋，精细的去接头要求可以考虑cutadapt或trimmomatic。

用户指定的cluster

通过指定--locifile或--locus可以提供cluster的信息。

De novo cluster 鉴定

越写越觉得翻译只是过程，理解才是目的。很多术语翻译过来有点奇怪的。（废话.gif）
ShortStack中识别sRNA cluster的标准：1. 找到包含至少一个primary alignment的所有区域，alignment的两端之间最大距离<=option --pad(default:75); 2. cluster中符合1的alignment的数目应当>=option --mincov(default:0.5rpm)

分析方法说明

不管是指定还是de novo，针对sRNA clusters的分析方法都是一样的。

特定Read group的counts

测序时：
样本建一个库在同一条lane上完成测序产生reads sets可定义为一个Read group;
样本建成分开独立的库测序得到的reads sets也就分属于不同的Read groups;
引自：# GATK - Read groups

结果文件说明

1.Counts.txt

Locus:sRNA cluster的坐标；
Name:sRNA cluster的名称；
main:所有Read groups的reads总数；

2.Results.txt

Locus:同上；
Name:同上；
Length:cluster的长度，单位是nt；
Reads:primary alignments(权重大)的总数目；
RPM:上一列reads的标准化值；
UniqueReads:unique reads的数目；

Uniquely mapped reads map to exactly one location within the reference genome.
引自：ENCODE

FracTop:比对到前导链(top strand)的primary alignments的分数(fraction)，对应分析参数--strand_cutoff，高于0.8被识别为+链，低于0.2被认为是-链；
Strand: +/-链识别的结果；
MajorRNA：cluster中reads数最多的RNA，如果第一不止一个，那就随机选择一个；
MajorRNAReads：MajorRNA的primary alignments数目；
Complexity: 计算公式为(n_distinct_read_sequences) / (abundance of all reads)，区间为(0,1]。数值越小说明该cluster中sRNA类型越少，反之亦同；
DicerCall:cluster中占主导地位的sRNA长度(20-24nt)，但如果超过80%的都不在此范围内，则被标记为N，这类RNA很有可能是降解产物；
MIRNA: ShortStack 在排除假阳性这一点上太狠了，15条标准逐步筛选都通过了才能被认定为Y(是miRNA)，不然哪一步失败了就会标记为N+失败步骤；

搜索标准如下：

标记	拒绝标准
N0	添加了--nohp参数就不会执行miRNA search
N1	比对上的reads为0
N2	在DicerCall范围内的reads不超过80%
N3	MajorRNAReads数目<2
N4	MajorRNA的长度不在DicerCall范围内
N5	locus size的长度>`--foldsize`
N6	locus比对不上任何一条链，即0.2<`--strand_cutoff`<0.8
N7	RNA folding attempt failed at locus
N8	可能的成熟miRNA所在链与locus相反，即为miRNA*
N9	Retrieval of possible mature miRNA position failed
N10	无法计算miRNA* 的一般性错误
N11	可能的成熟miRNA在预测的前体二级结构中具有> 5个未配对碱基
N12	单个预测的hairpin结构中未包含可能的成熟miRNA
N13	可能的miRNA/miRNA-star包含超过2个budges，且/或buldge>3nts
N14	Reads for possible miRNA, miRNA-star, and their 3p variants added up to less than 50% of the total reads at the locus
N15	如果miRNA* 没有被找到，一样是拒绝

通过所有筛选后，得到de novo annotation的new miRNA family.

14.PhaseScore：

A score of ~30 or more indicates a well-phased locus. Not all loci are subject to phasing analysis. Loci with no reads at all aligned, a DicerCall of anything except 21 or 24, a Locus Size of < 3*DicerCall, and stranded loci (>= 80% of reads on top strand OR <= 20% of reads on top strand) are not analyzed. These are assigned a PhaseScore of -1.

15.Short:primary alignments<--dicermin的reads数目；
16.Long:primary alignments>--dicermax的reads数目；
17结束：不同RNA size对应的primary alignments数目；

3.MIRNAs/

Original Location和Displayed Location的区别是什么？
里面有一些小写字母表示sRNA序列不同于参考序列的碱基。
l: length of RNA, a: number of alignments.

Below this top line, all other small RNAs aligned to the locus are shown. Those aligned to the opposite strand have '<' as delimiters instead of '.'.

4.ShortStack_All.gff3

所有的locus信息：ShortStack_D.gff3+ShortStack_N.gff3

5.ShortStack_D.gff3

Dicercall不为N的locus信息；

6.ShortStack_N.gff3

Dicercall为N的locus信息；

7.Unplaced.txt

未确定位置的sRNA有2种情况：1）归一化后的表达量低于--mincov的sRNA; 2)multi-mapping中无法指定位置的sRNA.