作者：Arno
審稿：童蒙
編輯：angelica
引言
Iso-seq測序經(jīng)初步分析獲得高質(zhì)量的轉(zhuǎn)錄本之后（全長轉(zhuǎn)錄本鑒定，全長轉(zhuǎn)錄本比對），便可以對轉(zhuǎn)錄本的結(jié)構(gòu)進行精確鑒定稿湿、注釋集峦。本次小編通過結(jié)合一個全長轉(zhuǎn)錄本后續(xù)分析的工具集Cupcake的使用方法，來給大家介紹一下Iso-seq測序得到高質(zhì)量轉(zhuǎn)錄本之后的一些分析评矩。
1. Cupcake安裝
在介紹Iso-seq后續(xù)分析之前叶堆，先來看一下Cupcake如何安裝使用。
Cupcake是一個Python和R腳本的合集斥杜，可直接通過python的形式安裝虱颗，其安裝方法可以參考下方代碼，目前均支持python2和python3版本蔗喂，以下示例均基于python2版本忘渔。
git clone https://github.com/Magdoll/cDNA_Cupcake.git
#export PATH=$PATH:<path_to_Cupcake>/sequence/
#export PATH=$PATH:<path_to_Cupcake>/rarefaction/
cd cDNA_Cupcake
#git checkout -b tofu2 tofu2_v21
python setup.py build
python setup.py install
2. 轉(zhuǎn)錄本去冗余
Iso-seq測序后經(jīng)過smrtlink和IsoSeq3軟件進行全長轉(zhuǎn)錄本鑒定后（cluster&polish），就獲得了高質(zhì)量的轉(zhuǎn)錄本（hq.isoform.fq）缰儿。這些高質(zhì)量轉(zhuǎn)錄本均為全長非嵌合的高質(zhì)量轉(zhuǎn)錄本（包含polyA畦粮，準確率>=0.99，至少2個全長序列支持）乖阵，但是因為IsoSeq3的聚類算法的敏感性和特異性宣赔，以及天然RNA5‘端的易降解的特性，得到的高質(zhì)量轉(zhuǎn)錄本中仍然存在冗余的轉(zhuǎn)錄本瞪浸，所以需要進一步去除儒将。
可使用Cupcake的“collapse_isoforms_by_sam.py”腳本，具體代碼以及說明可參考如下：
python ~/software/Python-2.7.8/bin/collapse_isoforms_by_sam.py --input sample.hq.fasta \
--dun-merge-5-shorter --sam sample.sort.sam --prefix sample_name --min-coverage 0.85 \
--min-identity 0.95 2>sample.collapse_isoforms.log
# --min-coverage --min-identity 為去冗余時的覆蓋率和一致性对蒲，默認為0.99和0.85可根據(jù)實際情況調(diào)整
# --dun-merge-5-shorter
# 得到的結(jié)果中sample.collapsed.group.txt為記錄合并冗余后的轉(zhuǎn)錄本信息,轉(zhuǎn)錄本格式為：PB.<loci_index>.<isoform_index>
# sample.ignored_ids.txt為去除的轉(zhuǎn)錄本信息
# sample.collapsed.rep.fq和sample.collapsed.gff分別為非冗余的轉(zhuǎn)錄本序列及其gff文件
此種去冗余方式是針對的有參考基因組序列的樣本椅棺，需要用到跟參考基因組比對的Sam文件，如果沒有參考基因組齐蔽，可以使用CD-HIT對序列進行聚類去冗余两疚，具體方式可參考：https://github.com/Magdoll/cDNA_Cupcake/wiki/Tutorial:-Collapse-redundant-isoforms-without-genome
3. 轉(zhuǎn)錄本定量
得到unique的轉(zhuǎn)錄本之后，再結(jié)合前邊聚類分析得到report文件cluster_report.csv是可以計算出來每個unique轉(zhuǎn)錄本的count數(shù)目的含滴，Cupcake提供了計算的腳本诱渤。但是目前PacBio測序得到的CCS對于做定量來說數(shù)據(jù)還是不太夠的，建議用ONT平臺測序的數(shù)據(jù)去做定量谈况，ONT的數(shù)據(jù)reads足夠長勺美，數(shù)據(jù)量足夠多的。
python ~/software/Python-2.7.8/bin/get_abundance_post_collapse.py sample.collapsed sample.cluster_report.csv
# sample.collapsed為樣本去冗余后的文件前綴名稱
4. 過濾5'端降解的轉(zhuǎn)錄本
去完冗余之后的轉(zhuǎn)錄本仍然存在一部分轉(zhuǎn)錄本比對到參考基因組的位置一致碑韵，但5'端長度不一致的轉(zhuǎn)錄本赡茸，這種情況是因為建庫過程中，使用的cDNA試劑盒并不會對5'端進行加帽處理祝闻，所以再整個過程中很可能會發(fā)生5'端的降解占卧，而發(fā)生降解的這些轉(zhuǎn)錄本是沒有任何生物學意義的，所以可以將5'端降解的轉(zhuǎn)錄本過濾掉，再用于后續(xù)的分析华蜒。過濾也可以使用Cupcake工具包提供的腳本辙纬。
python ~/software/Python-2.7.8/bin/filter_away_subset.py sample.collapsed
# sample.collapsed為樣本去冗余后的文件前綴名稱
# 得到輸出結(jié)果文件sample.collapsed.filtered.gff, sample.collapsed.filtered.abundance.txt,
# sample.collapsed.filtered.rep.fq
5. 融合基因分析
基因融合在基因組層面上可能由于基因組變異（染色體易位、中間缺失叭喜、染色體倒位）使得兩個不同基因的部分序列或全部序列融合到一起贺拣，形成一個新的基因，可能表達也可能不表達捂蕴；轉(zhuǎn)錄組層面上可能由于兩個基因轉(zhuǎn)錄產(chǎn)生的RNA譬涡，由于某種原因融合在一起，形成新的融合RNA啥辨，當然該RNA可能編碼蛋白也可能不編碼蛋白涡匀。
對于Iso-seq測序得到的轉(zhuǎn)錄本數(shù)據(jù)，尋找融合基因委可，可以采用Cupcake 中的“fusion_finder.py” 這個腳本進行渊跋，鑒定的默認標準有如下4點：
- 比對到2個或更多位置；
- 比對到的每一個位置至少覆蓋5%的轉(zhuǎn)錄本着倾；
- 融合轉(zhuǎn)錄本（各個位置的相加）比對率至少99%以上拾酝；
- 每一個比對位置的距離至少10kb以上。
## Best practice for fusion transcript finding
## https://github.com/Magdoll/cDNA_Cupcake/wiki/Best-practice-for-fusion-transcript-finding
gmap -D [dir] -d hg38 -f samse -n 0 input.fasta > input.fasta.gmap.sam
minimap2 -ax splice -uf --secondary=no hg38.fa input.fasta > input.fasta.minimap2.sam
sort -k 3,3 -k 4,4n input.fasta.minimap2.sam > input.fasta.minimap2.sorted.sam
fusion_finder.py --input input.fasta -s input.fasta.minimap2.sorted.sam \
--cluster_report cluster_report.csv \
-o output.fusion \
--min_locus_coverage_bp 500 -d 1000000
6. 結(jié)語
除了我們介紹的卡者，Cupcake有很多強大的cDNA序列分析功能蒿囤，關于其詳細的介紹可以查閱其githup倉庫：https://github.com/Magdoll/cDNA_Cupcake 。
