安裝

ShortStack依賴的軟件包括perl、samtools、bowtie、bowtie-build、gzip、RNAfold，檢查了一下，服務(wù)器中都符合要求，不需要額外安裝或更新了。

版本要求：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   #注意>>是向文件中追加內(nèi)容，>是覆蓋原有內(nèi)容
source ~/.bashrc

參數(shù)說明

ShortStack [參數(shù)] {readfile {bamfile cramfile}} genomefile
#bamfile和cramfile都是比對(duì)文件

• 一般參數(shù)
  --help
  --version
  --genomefile [string字符串] 給出基因組數(shù)據(jù)所在路徑
  --outdir [string] 給出輸出文件所在路徑
• 比對(duì)參數(shù)
  --readfile [string] 需要比對(duì)的reads文件路徑
  --adapter [string] 接頭文件，若不指定即默認(rèn)已經(jīng)去過接頭
  --bowtie_cores [integer] 傳遞給bowtie -p(線程數(shù))的參數(shù)。設(shè)定處理器核數(shù)，默認(rèn)為1。核與線程的關(guān)系：1個(gè)核對(duì)應(yīng)2個(gè)線程。
  --sort_mem [string] 傳遞給samtools的sort -m(內(nèi)存數(shù))的參數(shù)，即用sort對(duì)文件排序時(shí)使用的物理內(nèi)存空間。samtools中默認(rèn)是768M。
  --mismatches [integer] 比對(duì)中容忍的錯(cuò)配數(shù)，只能為0，1或2。默認(rèn)為1。如果有多個(gè)reads匹配上了，只保留錯(cuò)配數(shù)最低的reads。
  --cquals [string] 給出color-space類型測(cè)序數(shù)據(jù)的質(zhì)量值文件。color-space是一類特殊的測(cè)序數(shù)據(jù)，用顏色代表核苷酸，是Solid平臺(tái)產(chǎn)出數(shù)據(jù)的原始格式。
  --cram 默認(rèn)輸出sam格式，加上此參數(shù)可將比對(duì)結(jié)果輸出為cram格式。CRAM和BAM都是壓縮性質(zhì)的比對(duì)文件，CRAM在繼承了BAM優(yōu)點(diǎn)的基礎(chǔ)上，可對(duì)數(shù)據(jù)進(jìn)一步壓縮。
  --mmap [string] 處理multi-mapping reads的流程?？蛇xn (none), r (random), u (unique- seeded guide), or f(fractional-seeded guide). default: u
  --bowtie_m [string] max，指定最大比對(duì)到基因組的次數(shù)。超過50個(gè)alignments被標(biāo)記為unmapped。
  --ranmax [string] 當(dāng)--mmap指定u或f時(shí)，仍然可能出現(xiàn)多個(gè)權(quán)重最大的位置，這時(shí)候比對(duì)工具還是只能選擇隨機(jī)分配位置，但可能的位置數(shù)不能超過ranmax值。默認(rèn)為3。
  --align_only 出現(xiàn)此參數(shù)時(shí)，ShortStack只會(huì)執(zhí)行到比對(duì)這一步，不再進(jìn)行后續(xù)分析。
  --show_secondaries 出現(xiàn)此參數(shù)時(shí)，比對(duì)文件中除了有primary alignments，還有secondary alignments信息。
  --keep_quals 出現(xiàn)此參數(shù)時(shí)，會(huì)保留比對(duì)質(zhì)量值信息。
• 分析參數(shù)
  --bamfile sRNA的比對(duì)文件，和readfile、cramfile選一個(gè)作為輸入文件就行。
  --cramfile 同上。
  --dicermin Dicer酶切割后的sRNA的最小size。一般大于15，默認(rèn)為20.
  --dicermax Dicer酶切割后的sRNA的最大size。一般大于15，默認(rèn)為24.
  --foldsize 提取參考基因組序列進(jìn)行RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)的長度。范圍：200-1000.默認(rèn)為300.增加數(shù)值后運(yùn)行時(shí)間將大大增加，但也可能找到更多MIRNAs.這里指的應(yīng)該是miRNA的前體，一般認(rèn)為miRNA的前體長度不超過300nt。
  --locifile 包含特定的intervals信息的文本文件。和--locus互斥。
  
  loci文件示例
  
  --locus 同上
  --nohp 禁用miRNA search。
  --pad 如果sRNA clusters之間的距離小于等于pad值，sRNA clusters就會(huì)被合并。默認(rèn)值：75.
  --mincov sRNA clusters表達(dá)量的最低值。默認(rèn)值：0.5rpm
  --strand_cutoff strand的臨界值。默認(rèn)為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也就是與其互補(bǔ)的另一條鏈。注意與正義鏈、反義鏈的區(qū)分。

--total_primaries 直接告訴ShortStack bam文件中的primary alignments的數(shù)目，減少計(jì)算時(shí)間，包括沒有比對(duì)上的。

作者最后給出了一些運(yùn)行建議：比如bowtie -m和 --ranmax參數(shù)設(shè)置對(duì)文件size影響很大。

比對(duì)方法說明

ShortStack可以把比對(duì)過的文件作為輸入文件，也可以使用其內(nèi)置的比對(duì)方法：Improved Placement of Multi-mapping Small RNAs（doi:10.1534/g3.116.030452）
主要解決multi-mapping reads 如何放到合適位置的問題。

拓展01：multi-mapping reads在sRNA-seq里很常見，一是因?yàn)閟RNA reads很短，二是sRNA傾向于從多拷貝區(qū)域起源。以前的序列比對(duì)工具要么準(zhǔn)確性低（隨機(jī)選取一個(gè)匹配分?jǐn)?shù)高的位置），要么敏感性低（搞不清楚干脆就不要了，忽略掉所有multi-mapping reads），bowtie就是代表，它默認(rèn)隨機(jī)選擇，也可以設(shè)置主動(dòng)忽略。
拓展02：比對(duì)方法的基本原理——1)用bowtie找出每條read的best-scoring alignments，每條read最多允許有50個(gè)alignments；2）采用local-weighting為每個(gè)alignment計(jì)算概率；3）將概率作為權(quán)重，權(quán)重大的作為primary alignments，其余的作為secondary alignments；4)比對(duì)模式選擇，講了這么多就是告訴你選U最好。
拓展03：這里比對(duì)指的是比對(duì)到sRNA轉(zhuǎn)錄起源的地方，而非target sites。靶位點(diǎn)的預(yù)測(cè)一般會(huì)有一些species-specific parameters以及只考慮基因組的部分區(qū)域（如成熟轉(zhuǎn)錄區(qū)）。

ShortStack比對(duì)方法說明

基因組預(yù)處理

基因組文件需為fasta格式；
染色體名最好不要有空格，否則只保留空格前面的部分；
如果基因組擁有>50chromosomes/scaffolds/contigs，且N50<1Mb，ShortStack會(huì)默認(rèn)拼接短的染色體。但在3.8之后的版本中，報(bào)告的搜索結(jié)果仍然是基于原始基因組的。

Reads預(yù)處理

sRNA文件可以為壓縮格式；
同時(shí)將多個(gè)文件作為輸入文件，之間用逗號(hào)隔開；
不支持雙端文件；
No condensation?
所有reads必須有唯一名；
ShortStack也可以通過指定--adapter去接頭，不過可能比較簡陋，精細(xì)的去接頭要求可以考慮cutadapt或trimmomatic。

用戶指定的cluster

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

De novo cluster 鑒定

越寫越覺得翻譯只是過程，理解才是目的。很多術(shù)語翻譯過來有點(diǎn)奇怪的。（廢話.gif）
ShortStack中識(shí)別sRNA cluster的標(biāo)準(zhǔn)：1. 找到包含至少一個(gè)primary alignment的所有區(qū)域，alignment的兩端之間最大距離<=option --pad(default:75); 2. cluster中符合1的alignment的數(shù)目應(yīng)當(dāng)>=option --mincov(default:0.5rpm)

分析方法說明

不管是指定還是de novo，針對(duì)sRNA clusters的分析方法都是一樣的。

特定Read group的counts

測(cè)序時(shí)：
樣本建一個(gè)庫在同一條lane上完成測(cè)序產(chǎn)生reads sets可定義為一個(gè)Read group;
樣本建成分開獨(dú)立的庫測(cè)序得到的reads sets也就分屬于不同的Read groups;
引自：# GATK - Read groups

結(jié)果文件說明

1.Counts.txt

Locus:sRNA cluster的坐標(biāo)；
Name:sRNA cluster的名稱；
main:所有Read groups的reads總數(shù)；

2.Results.txt

1. Locus:同上；
2. Name:同上；
3. Length:cluster的長度，單位是nt；
4. Reads:primary alignments(權(quán)重大)的總數(shù)目；
5. RPM:上一列reads的標(biāo)準(zhǔn)化值；
6. UniqueReads:unique reads的數(shù)目；

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

7. FracTop:比對(duì)到前導(dǎo)鏈(top strand)的primary alignments的分?jǐn)?shù)(fraction)，對(duì)應(yīng)分析參數(shù)--strand_cutoff，高于0.8被識(shí)別為+鏈，低于0.2被認(rèn)為是-鏈；
8. Strand: +/-鏈識(shí)別的結(jié)果；
9. MajorRNA：cluster中reads數(shù)最多的RNA，如果第一不止一個(gè)，那就隨機(jī)選擇一個(gè)；
10. MajorRNAReads：MajorRNA的primary alignments數(shù)目；
11. Complexity: 計(jì)算公式為(n_distinct_read_sequences) / (abundance of all reads)，區(qū)間為(0,1]。數(shù)值越小說明該cluster中sRNA類型越少，反之亦同；
12. DicerCall:cluster中占主導(dǎo)地位的sRNA長度(20-24nt)，但如果超過80%的都不在此范圍內(nèi)，則被標(biāo)記為N，這類RNA很有可能是降解產(chǎn)物；
13. MIRNA: ShortStack 在排除假陽性這一點(diǎn)上太狠了，15條標(biāo)準(zhǔn)逐步篩選都通過了才能被認(rèn)定為Y(是miRNA)，不然哪一步失敗了就會(huì)標(biāo)記為N+失敗步驟；

搜索標(biāo)準(zhǔn)如下：

通過所有篩選后，得到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數(shù)目；
16.Long:primary alignments>--dicermax的reads數(shù)目；
17結(jié)束：不同RNA size對(duì)應(yīng)的primary alignments數(shù)目；

3.MIRNAs/

Original Location和Displayed Location的區(qū)別是什么？
里面有一些小寫字母表示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）歸一化后的表達(dá)量低于--mincov的sRNA; 2)multi-mapping中無法指定位置的sRNA.

8.ErrorLogs.txt

9. Log.txt

運(yùn)行記錄，和打印到屏幕上的信息一樣。

參考資料

https://github.com/MikeAxtell/ShortStack