作者：Arno
審稿：童蒙
編輯：angelica
前面我們介紹了PacBio三代全長轉(zhuǎn)錄組測序相關(guān)的全長轉(zhuǎn)錄本鑒定、全長轉(zhuǎn)錄本比對、全長轉(zhuǎn)錄本結(jié)構(gòu)分析上篇。
今天我們繼續(xù)介紹包括新轉(zhuǎn)錄本鑒定、可變剪切以及可變多聚腺苷酸化APA等全長轉(zhuǎn)錄本結(jié)構(gòu)分析。
新轉(zhuǎn)錄本鑒定
通過將去除冗余后的unique轉(zhuǎn)錄本與參考基因組進(jìn)行比較，可以對轉(zhuǎn)錄本進(jìn)行結(jié)構(gòu)注釋，從而可以發(fā)現(xiàn)新的未知的轉(zhuǎn)錄本。
MatchAnnot軟件是一款可以將比對結(jié)果跟注釋文件或者注釋文件和注釋文件進(jìn)行比較的Python軟件，可以鑒定已知和新的全長轉(zhuǎn)錄本，同時(shí)基于其輸出結(jié)果還可以進(jìn)行基因的可視化。下面我們看看具體怎么使用。
1. 比對
MatchAnnot需要比對并排序后的sam文件作為輸入，所以再運(yùn)行之前，需要先進(jìn)行全長轉(zhuǎn)錄本的比對，可以參考我們之前介紹的比對方面，下面給大家再提供一個(gè)示例。
## step1. 比對（可以使用gmap或minimap2均可）
gmap -D [dir] -d hg38 -f samse -n 0 sample.hq.fasta > sample.raw.sam
minimap2 -ax splice -uf --secondary=no hg38.fa sample.hq.fasta > sample.raw.sam
samtools view -bS -o sample.raw.bam sample.raw.sam
samtools sort -@ 8 sample.raw.bam ./sample
/samtools index sample.bam
samtools view -h sample.bam sample.sam
2. 運(yùn)行MatchAnnot
得到排序后的sam文件后，進(jìn)行MatchAnnot分析。MatchAnnot會(huì)對三代測序到轉(zhuǎn)錄本，逐個(gè)exon進(jìn)行分析，找出注釋文件中對應(yīng)的所有轉(zhuǎn)錄本，以及新發(fā)現(xiàn)的轉(zhuǎn)錄本。具體的輸出文件解讀或格式，可以參考：https://github.com/TomSkelly/MatchAnnot/wiki/How-to-Interpret-matchAnnot-Output
## step2. MatchAnnot
python ~/MatchAnnot-master/matchAnnot.py --gtf hg38.gtf --format alt sample.sam --outpickle sample.pick >sample.matchAnnot.xls
# --format gtf的格式默認(rèn)為GENCODE標(biāo)準(zhǔn)的GTF格式，同時(shí)也支持alternate format(包含更多信息的gtf格式)或者python的pickle格式
# --outpickle 如果需要進(jìn)行后續(xù)的可視化，需要輸出pickle格式
# 輸入的Sam文件需要是按照染色體排序的
3. 結(jié)果可視化
為了更好的幫助理解MatchAnnot的輸出結(jié)果，該軟件同時(shí)提供了可視化的方法，可以針對性的可視化感興趣的基因。運(yùn)行示例，以及結(jié)果示例如下：
## step3 clusterView
python ~/MatchAnnot-master/clusterView.py --gtf hg38.gtf --format alt \
--gene=target_gene --matches sample.pickle --output target_gene.png \
--title "target gene"
可變剪切鑒定
真核生物中，基因轉(zhuǎn)錄產(chǎn)生的mRNA前體可以通過外顯子跳躍、內(nèi)含子保留等不同的剪切形式產(chǎn)生多種轉(zhuǎn)錄本異構(gòu)體（isoforms），大大增加了轉(zhuǎn)錄本多樣性。全長轉(zhuǎn)錄組測序憑借讀長優(yōu)勢，能夠直接獲得由5’端至3’端poly(A)尾的完整mRNA序列，從而可以準(zhǔn)確鑒定基因的不同剪切形式的轉(zhuǎn)錄本。
可變剪切（Alternative splicing），通?？梢詣澐忠韵聨追N常規(guī)類型：
- 外顯子跳躍(Skipped exon, SE)
- 5'端可變外顯子(Alternative 5’ splice site, A5SS)
- 3'端可變外顯子(Alternative 3’ splice site, A3SS)
- 內(nèi)含子保留(Intron Retention, IR)
- 互斥外顯子(Mutually Exclusive Exons, MEE)
除此常規(guī)類型之外，還有可變轉(zhuǎn)錄起始位點(diǎn)、可變轉(zhuǎn)錄終止位點(diǎn)等類型?？傊虻目勺兗羟行问蕉喾N多樣，從而導(dǎo)致轉(zhuǎn)錄本和蛋白結(jié)構(gòu)與功能的多態(tài)性，是一種重要的調(diào)控機(jī)制。
目前常見的全長轉(zhuǎn)錄本的可變剪切分析工具有：
- SpliceMap-LSC-IDP pipeline
- SUPPA2
- AStalavista等。
AStalavista相較SpliceMap-LSC-IDP pipeline工具檢出效率要高很多，建議大家使用AStalavista。
AStalavista（alternative splicing and transcriptional landscape visualization）工具的分析，基于鑒定的轉(zhuǎn)錄本的gtf注釋文件，可以在線（http://astalavista.sammeth.net/ ）或者本地命令行兩種使用形式，使用簡單方便。其中，本地命令行執(zhí)行方式可參考如下：
~/AStalavista/astalavista-4.0/bin/astalavista -t asta -i sample.transcript.gtf --threads 2 -o result.gtf 1>./astalavista.o 2>&1
## -t 指定astalavista使用的分析工具，默認(rèn)為asta，進(jìn)行可變剪切事件鑒定，另外還可以使用sortBED、sortGTF、subsetter等工具
## -i 輸入文件
## --threads使用線程數(shù)
## 結(jié)果輸出形式為gtf形式，其中structure部分，分別用數(shù)字、"^","-"符號(hào)表示可變剪切發(fā)生的相對位置、供體和受體位點(diǎn)。如'0,1-2^'代表SE類型；'1-,2-'代表A3SS類型；'1^,2^'代表A5SS類型；'1-2^,3-4^'代表MXE類型；'0,1^2-'代表IR類型等。
SUPPA2軟件也可基于全長轉(zhuǎn)錄本gtf文件進(jìn)行可變剪切事件鑒定，同時(shí)可以基于二代數(shù)據(jù)，進(jìn)行可變剪切定量以及差異可變剪切分析，功能比較強(qiáng)大。
## step1. 基于二代數(shù)據(jù)進(jìn)行轉(zhuǎn)錄本定量
~/salmon/bin/salmon index -t unigene.fa -i unigene_index
~/salmon/bin/salmon quant -i unigene_index -l ISF --gcBias -1 R1.fq -2 R2.fq -p 4 -o sample
python ~/SUPPA/multipleFieldSelection.py -i ~/Salmon_output/*/quant.sf -k 1 -f 4 -o sample_tpm.txt ## 整理salmon結(jié)果
## step2. 使用generateEvents命令根據(jù)全長轉(zhuǎn)錄本gtf文件生成所有的可變剪切事件，結(jié)果為ioe格式
python ~/SUPPA/suppa.py generateEvents -i unigene.transcript.gtf -o sample.events -e SE SS MX RI FL -f ioe
### -i 輸入的gtf文件
### -o 輸出的文件前綴
### -e 輸出可變剪切的類型
### -f 設(shè)置輸出格式,將不同的可變剪切事件合并成一個(gè)結(jié)果
## step3. 計(jì)算PSI值,結(jié)果為.psi格式
python ~/SUPPA/suppa.py psiPerEvent -i sample.events.ioe -e sample_tpm.txt -o sample_events
## step4. 兩個(gè)樣本間的差異分析，結(jié)果為.dpsi格式
python ~/SUPPA/suppa.py diffSplice -m empirical -gc -i unigene.transcript.gtf -p sample1_events.psi sample2_events.psi -e sample1_tpm.txt sample2_tpm.txt -o sample_diffSplice
## step5. 對可變剪切事件進(jìn)行聚類分析
python ~/SUPPA/suppa.py clusterEvents --dpsi sample_diffSplice.dpsi --psivec sample.psivec --sig-threshold 0.05 --eps 0.2 --separation 0.11 -dt 0.2 --min-pts 10 --groups 1-2,4-6 -c OPTICS -o ./
可變多聚腺苷酸化APA分析
可變多聚腺苷酸化（alternative polyadenylation，APA）是指一個(gè)基因上有多個(gè)多聚腺苷酸化位點(diǎn)，使得一個(gè)基因可以產(chǎn)生多條帶有不同長度3’UTR的mRNA，即polyA尾長度不一致，可能產(chǎn)生不同編碼序列的轉(zhuǎn)錄本，從而使得轉(zhuǎn)錄本存在多樣性。
APA是一種非常常見的轉(zhuǎn)錄后修飾和調(diào)控方式，polyA的長度對于mRNA的穩(wěn)定性以及蛋白翻譯影響很大，太短的mRNA穩(wěn)定性較差，不同物種間的polyA尾長度差異也比較大。
一般APA可以分為四種類型：
1 3’UTR APA：發(fā)生在末端外顯子內(nèi)，產(chǎn)生具有不同長度3’UTR的轉(zhuǎn)錄本，不影響蛋白編碼功能，是最常見的APA形式；
2 可變末端外顯子APA：這種APA產(chǎn)生了末端外顯子不同的轉(zhuǎn)錄本，影響蛋白編碼功能；
3 內(nèi)含子APA：發(fā)生于在內(nèi)含子區(qū)，延長了某個(gè)內(nèi)部外顯子并使之成為末端外顯子；
4 內(nèi)部外顯子APA：在編碼區(qū)域內(nèi)部發(fā)生剪切和多聚腺苷酸化。
TAPIS(Transriptome AnalysisPipeline from Isoform Sequencing)可以用來做全長轉(zhuǎn)錄組可變剪切以及APA分析，研究者使用此軟件研究了高粱的可變polyA。依賴于Python2.7，其可用于三代測序數(shù)據(jù)的糾錯(cuò)、比對、鑒定可變剪切以及識(shí)別APA位點(diǎn)，原理流程如下圖所示。
其用于鑒定APA位點(diǎn)的使用方法可參考如下。
## step1 alignPacBio.py 將全長轉(zhuǎn)錄本回帖到基因組，依賴于Gmap軟件比對
usage: alignPacBio.py [-h] [-v] [-i ITERATIONS] [-e EDR] [-o OUTDIR]
[-p PROCS] [-K MAXINTRON]
indexesDir indexName reference fasta
python alignPacBio.py -p 4 -o ./sample ~/Genome_Index ref ref.genome.fa sample.collapsed.hq.fa
# ~/Genome_Index 為GMAP軟件比對建立的參考基因組索引所在文件夾
# ref 索引名稱
# ref.genome.fa 索引序列
# sample.collapsed.hq.fa 去冗余的轉(zhuǎn)錄本
## step2 run_tapis.py 分析可變剪切和APA
usage: run_tapis.py [-h] [-v] [-p] [-o OUTDIR] [-t TRIMMAX] [-w W]
[-m MINDIST] [-s MINSUPPORT]
geneModel bamfile
python run_tapis.py -o ./sample ref.gtf ~/Sampe/aligned.bam
# 兩個(gè)輸入文件：參考基因組的注釋文件以及前一步的比對結(jié)果
# 輸出文件為 assembled.gtf 和 polyA_summary.csv 記錄了每個(gè)基因polyA的個(gè)數(shù)以及位置
參考文獻(xiàn)
- Foissac S, Sammeth M (2007) ASTALAVISTA: dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic Acids Research 35 (Web Server issue): W297-299
- https://github.com/comprna/SUPPA/wiki/SUPPA2-tutorial#differential-splicing-with-local-events
- Abdel-Ghany,S. E., Hamilton, M., Jacobi, J. L., Ngam, P., Devitt, N., Schilkey, F.,Ben-Hur, A., and Reddy, A. S. N. A survey of the sorghum transcriptome usingsingle-molecule long reads. Nature Communications, 2016 7:11706.
- https://github.com/TomSkelly/MatchAnnot/wiki/How-to-Interpret-matchAnnot-Output
