PacBio測(cè)序表現(xiàn)出了非常強(qiáng)大的優(yōu)勢(shì)，比如測(cè)序read長(zhǎng)、無(wú)GC偏好、直接檢測(cè)變異、直接檢測(cè)堿基修飾等。尤其在轉(zhuǎn)錄中的應(yīng)用已經(jīng)相對(duì)成熟，下面重點(diǎn)介紹如何去進(jìn)行序列的比對(duì)和去冗余。
關(guān)鍵詞：PacBio、minimap2、GMAP、collapse、full length isoforms
因?yàn)槲覀冎繧so-seq測(cè)序后經(jīng)過(guò)smrtlink和IsoSeq3軟件的處理會(huì)得到高質(zhì)量的全長(zhǎng)轉(zhuǎn)錄本序列，針對(duì)有參考基因組的物種，首先要進(jìn)行序列的回帖（比對(duì)）。
回帖軟件有minimap2, GMAP, STAR, BLAT等，重點(diǎn)介紹minimap2和GMAP。
一、比對(duì)
1. minimap2
軟件鏈接：https://github.com/lh3/minimap2
該軟件支持剪切比對(duì)和非剪切比對(duì)，所以非常適合轉(zhuǎn)錄本的比對(duì)，推薦使用2.9及以上版本，且支持不建索引模式的比對(duì)。
使用示例如下：
minimap2 -t 30 -ax splice -uf --secondary=no -C5 hg38.fasta hq_isoforms.fasta? > hq_isoforms.fasta.sam? 2> hq_isoforms.fasta.sam.log
-ax spliced比對(duì)模式
--secondary=no 只輸出最好的比對(duì)結(jié)果
2.GMAP
軟件鏈接：http://research-pub.gene.com/gmap/
索引文件必需提前建好（gmap_build命令），推薦使用2018-03-20或更高版本（版本以日期命名）。
使用示例：
gmap -D /gmap_db/ -d hg38 -f samse -n 0 -t 16? --cross-species --max-intronlength-ends 200000 -z sense_force hq_isoforms.fasta > hq_isoforms.fasta.sam? 2> hq_isoforms.fasta.sam.log
1和2步驟中，獲得了比對(duì)文件后（SAM格式）可進(jìn)一步轉(zhuǎn)換為BAM格式：
samtools view -bS hq_isoforms.fasta.sam > hq_isoforms.fasta.bamsamtools sort hq_isoforms.fasta.bam > hq_isoforms.fasta.sorted.bamsamtools index hq_isoforms.fasta.sorted.bam
二、去冗余
比對(duì)文件中記錄了低質(zhì)量比對(duì)和同于的基因、轉(zhuǎn)錄本異構(gòu)體，需要進(jìn)一步進(jìn)行過(guò)濾。
是的，沒錯(cuò)，CupCake可以處理。
軟件鏈接：https://github.com/Magdoll/cDNA_Cupcake
直接使用里面的collapse_isoforms_by_sam.py進(jìn)行處理：
usage: collapse_isoforms_by_sam.py [-h]
[--input INPUT] [--fq] -s SAM -o PREFIX [-c MIN_ALN_COVERAGE] [-i MIN_ALN_IDENTITY] [--max_fuzzy_junction MAX_FUZZY_JUNCTION] [--flnc_coverage FLNC_COVERAGE] [--dun-merge-5-shorter]
比對(duì)處理和去冗余：
gmap -D <gmap_db_location> -d hg38 -f samse -n 0 -t 12? -z sense_force hq_isoforms.fastq > hq_isoforms.fastq.sam
sort -k 3,3 -k 4,4n hq_isoforms.fastq.sam > hq_isoforms.fastq.sorted.sam
collapse_isoforms_by_sam.py --input hq_isoforms.fastq --fq -s hq_isoforms.fastq.sorted.sam --dun-merge-5-shorter -o test -c 0.95 -i 0.85
后續(xù)可以通過(guò)再IGV中觀察hq_isoforms.fastq.bam和生成的去冗余后的test.collapsed.gff文件。
另外，
test.group.txt為記錄合并冗余后的對(duì)應(yīng)文件，可以知道一共保留了多少個(gè)非冗余轉(zhuǎn)錄本。test.ignored_ids.txt為沒有比對(duì)上而被丟棄的轉(zhuǎn)錄本編號(hào)。
參考：
https://github.com/Magdoll/cDNA_Cupcake/wiki/Best-practice-for-aligning-Iso-Seq-to-reference-genome:-minimap2,-GMAP,-STAR,-BLAT#refgmap
https://github.com/Magdoll/cDNA_Cupcake/wiki/Cupcake-ToFU:-supporting-scripts-for-Iso-Seq-after-clustering-step#collapse
