11、比對軟件STAR（https://github.com/alexdobin/STAR）

轉載：https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484731&idx=1&sn=b15fbee5910b36341bf366860ee5df53&scene=21#wechat_redirect

這次給大家帶來的是ENCODE project的御用比對軟件STAR，ENCODE項目是一個由美國國家人類基因組研究所(NHGRI)在2003年9月發起的一項公共聯合研究項目，旨在找出人類基因組中所有功能組件[。這是既完成人類基因組計划后國家人類基因組研究所開始的最重要的項目之一。所有在該項目中產生的數據都會被迅速的在公共數據庫中公開。
在我之前的那篇RNA-seq數據分析—-方法學文章的實戰練習文章里關於比對軟件的比較中STAR也展現了不俗的表現。所以在處理比對時我也考慮了將HISAT2與STAR共同使用，查看它們的表現情況，選取適合的比對工具。

---

STAR的安裝

cd biosoft && mkdir STAR && cd STAR
wget https://github.com/alexdobin/STAR/archive/2.5.3a.tar.gz
tar -xzf 2.5.3a.tar.gz
cd STAR-2.5.3a

# for easy use, add bin/ to your PATH

下載需要參考基因組並進行index構建

# downloading dna index fasta file
nohup wget -r -np -nH -nd -R index.html -L ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna_index/ &

# download gft annotation file
nohup wget ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/Homo_sapiens.GRCh38.90.chr_patch_hapl_scaff.gtf.gz &

mkdir STAR_index && cd STAR_index
STAR --runMode genomeGenerate --genomeDir ~/reference/STAR_index/ --genomeFastaFiles ~/reference/genome/hg38/Homo_sapiens.GRCh38.dna.toplevel.fa --sjdbGTFfile ~/reference/genome/hg38/Homo_sapiens.GRCh38.90.chr_patch_hapl_scaff.gtf --sjdbOverhang 199

# --sjdbOverhang 數值為reads長度-1
# Mode 為generate
# --genomeFastaFiles --sjdbGTFfile 分別對應fasta文件和GTF文件

---

STAR的使用

# STAR的manual里面給了最基本的比對參數示例
STAR
--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2 ]

# 基本示例，
針對fastq.gz文件增加--readFilesCommand gunzip -c 參數/--readFilesCommand zcat參數，針對bzip2文件使用--readFilesCommand bunzip2 -c參數
STAR --runThreadN 20 --genomeDir ~/reference/STAR_index/ --readFilesCommand zcat --readFilesIn ~/RNA-seq/LiuPing_data/RNA-seq/SC_w2q20m35_N_1.fq.gz ~/RNA-seq/LiuPing_data/RNA-seq/SC_w2q20m35_N_2.fq.gz

# 輸出unsorted or sorted bam file
--outSAMtype BAM Unsorted 實際上就是-name 的sort，下游可以直接接HTSeq
--outSAMtype BAM SortedByCoordinate
--outSAMtype BAM Unsorted SortedByCoordinate 兩者都輸出

額外參數說明

# 單獨指定注釋文件，而不用在構建的時候使用
--sjdbGTFfile /path/to/ann.gtf
--sjdbFileChrStartEnd /path/to/sj.tab

# ENCODE參數

# 減少偽junction的幾率
--outFilterType BySJout

# 最多允許一個reads被匹配到多少個地方
--outFilterMultimapNmax 20

# 在未有注釋的junction區域，最低允許突出多少個bp的單鏈序列
--alignSJoverhangMin 8

# 在有注釋的junction區域，最低允許突出多少個bp的單鏈序列
--alignSJDBoverhangMin 1

# 過濾掉每個paired read mismatch數目超過N的數據，999代表着忽略這個過濾
--outFilterMismatchNmax 999

# 相對paired read長度可以允許的mismatch數目，如果read長度為100，數值設定為0.04，則會過濾掉100*2*0.04=8個以上的數據
--outFilterMismatchNoverReadLmax 0.04

# 最小的intro長度
--alignIntronMin 20

# 最大的intro長度
--alignIntronMax 1000000

# maximum genomic distance between mates，翻譯不出來，自行理解
--alignMatesGapMax 1000000

---

STAR的輸出

STAR可以根據你的參數設定輸出多個結果文件，包含各種信息，下面對默認參數情況下的輸出文件做了一個詳細的展示，有些不好翻譯的地方我選擇使用原汁原味的manual text

• Aligned.out.sam
  Aligned.out.sam當然就是我們的比對結果啦！

E00516:168:H37WKCCXY:8:1101:6400:59130    99    1    92836373    255    20M1063N129M    =    92837548    4244    GGCTTGTCTATCCCTCACAGTACCAAACGATTCCCTGGTTATGATTCTGAAAGCAAGGAATTTAATGCAGAAGTACATCGGAAGCACATCATGGGCCAGAATGTTGCAGATTACATGCGCTACTTAATGGAAGAAGATGAAGATGCTTA    AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ    NH:i:1    HI:i:1    AS:i:289    nM:i:0

# 我截取了一條比對信息
我們來看一下最后面的 NH:i:1  HI:i:1  AS:i:289    nM:i:0
NH:i:后面的數值代表着此條read比對到幾個loci，1代表着unique map，數值大於1代表着multi-mappers
HI:i:后面的數值attrbiutes enumerates multiple
alignments of a read starting with 1，下游分析接cufflinks or stringtie的時候需要使用參數--outSAMattrIHstart 0
AS:i:的數值代表着local alignment score (paired for paired-edn reads)
nM:i:的數值代表着the number of mismatches per (paired) alignment, not to be confused with NM, which is the number of mismatches in each mate
關於下游處理工具的兼容性還需要使用者自己仔細參考manual

• Log.out文件
  Log.out文件記錄了程序運行時的信息，可以用來回溯錯誤信息。

tail Log.out
Joined thread # 12
Completed: thread #13
Joined thread # 13
Joined thread # 14
Joined thread # 15
Joined thread # 16
Joined thread # 17
Joined thread # 18
Joined thread # 19
ALL DONE!

• Log.progress.out文件
  Log.progress.out報告比對進程情況，1分鍾記錄一次

tail Log.progress.out
Sep 08 17:57:52     33.1    23115987      285    94.1%    284.0     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 17:58:53     34.0    24349711      285    94.1%    284.0     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:00:23     33.5    24789186      285    94.1%    284.1     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:01:51     33.3    25493588      285    94.1%    284.0     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:02:58     33.5    26284824      285    94.1%    284.1     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:04:23     33.7    27163519      285    94.1%    284.1     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:05:36     33.1    27428080      285    94.1%    284.1     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:06:54     33.8    28659661      285    94.1%    284.1     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:08:00     34.3    29741743      285    94.1%    283.9     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
ALL DONE!

• Log.final.out文件
  Log.final.out包含了比對結束后比對統計的信息

head Log.progress.out 
           Time    Speed        Read     Read   Mapped   Mapped   Mapped   Mapped Unmapped Unmapped Unmapped Unmapped
                    M/hr      number   length   unique   length   MMrate    multi   multi+       MM    short    other
Sep 08 17:17:47      2.9       88583      288    94.2%    287.4     0.1%     4.0%     0.1%     0.0%     1.7%     0.0%
Sep 08 17:18:53     14.5      711158      282    94.1%    281.9     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
Sep 08 18:08:00     34.3    29741743      285    94.1%    283.9     0.2%     4.0%     0.1%     0.0%     1.8%     0.0%
ALL DONE!

• SJ.out.tab文件
  SJ.out.tab包含了剪切信息，其實目前我還沒怎么看懂，等以后再來補坑。

head SJ.out.tab 
1    14830    14969    2    2    0    1    9    69
1    14844    14969    2    2    0    0    2    30
1    15039    15795    2    2    1    2    7    53
1    15948    16606    2    2    1    1    1    41
1    16028    16606    2    2    0    0    1    57
1    16311    16606    2    2    0    2    0    67
1    16766    16853    2    2    0    2    0    43
1    16766    16857    2    2    1    17    108    73
1    16766    16875    2    2    0    0    1    61
1    16789    16875    2    2    0    0    1    53

# 參數釋義
column 1: chromosome
column 2: first base of the intron (1-based)
column 3: last base of the intron (1-based)
column 4: strand (0: undened, 1: +, 2: -)
column 5: intron motif: 0: non-canonical; 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5:AT/AC, 6: GT/AT
column 6: 0: unannotated, 1: annotated (only if splice junctions database is used)
column 7: number of uniquely mapping reads crossing the junction
column 8: number of multi-mapping reads crossing the junction
column 9: maximum spliced alignment overhang

---

寫在最后

其實我探究STAR的最終目的實現利用STAR的Chimeric and circular alignments. 我自己處理的數據里面存在着fusion-protein，而其余的比對軟件暫時還沒發現有這個功能的

當使用—chimSegmentMin參數的時候，STAR可以把read拆分為兩部分，分別進行比對

STAR-Fusion是一個package，可以承接STAR的chimeric output，點我看代碼

當然STAR還可以做2-pass mapping，可以detect more splicesreads mapping to novel junctions

使用—quantMode GeneCounts參數還可以達到HTSeq的效果哦，可以幫你生成count matrix，省去你HTSeq的功夫， 有空回來做一個比對，看HTSeq和GeneCounts的效率。