全长转录组 | Iso-Seq 三代测序数据分析流程 (Pa

通过isoseq collapse和pigeon，我们最终能得到带有基因注释的表达矩阵，后续基因的差异分析，差异基因的KEGG注释等都和二代RNA-seq类似了。

Isoseq 数据分析第一部分我们最后使用了isoseq cluster 获得了聚类后高质量的转录本，但是我们仍然不知道这些经过聚类的转录本在基因组的位置以及属于哪些基因？这些转录本是已经注释的还是新的isoform？每个聚类是否能够进一步合并？每个isoform的表达量情况？下面我们通过使用isoseq collapse和 pigeon对转录本（isoforms）进行在参考基因组指导下的进一步合并（collapse），注释，分类和定量。

一、Iso-Seq Collapse

在isoseq cluster完成以后，我们首先需要将高质量全长isoforms回贴到参考基因组上，然后进行isoseq collapse。

图1. Map and Collapse

加上--do-not-allow-extra-5exons和默认条件下Isoforms合并原则 （图2）。

图2. 5'端外显子合并原则

下面三种参数下的合并（collapse）原则（图3）：
--max-fuzzy-junction INT Ignore mismatches or indels shorter than or equal to N. [5]
--max-5p-diff INT Maximum allowed 5' difference if on same exon. [50]
--max-3p-diff INT Maximum allowed 3' difference if on same exon. [100]

图3. Collapse Maximum 5p/3p Differences

• 下载参考基因组
  这里我们使用gencode版的人类参考基因组作为示例，任何版本的参考基因组原理类似。按照图4所示下载相应的人类参考基因组以及注释文件，并重新命名为GRCh38.fa和GRCh38.gtf。

图4. 人类参考基因组和注释文件的下载

单个样本数据示例

• pbmm2回贴转录本，map reads using pbmm2

# 建立参考基因组文件夹，并将参考基因组和注释文件放入其中
$ mkdir human_ref
#参考序列建立索引
$ pbmm2 index GRCh38.fa GRCh38.mmi

# 单个样本数据示例

$ pbmm2 align --preset ISOSEQ --sort <input.bam> <ref.fa> <mapped.bam>  #使用说明

$ pbmm2 align --preset ISOSEQ --sort human_80K.transcripts.bam human_ref/GRCh38.fa human_80K.mapped.bam

• 对回贴转录本进行collapse， collapse mapped reads into unique isoforms

# 单个样本数据示例
$ isoseq collapse --do-not-collapse-extra-5exons <mapped.bam> <flnc.bam> <collapsed.gff>  #使用说明
$ isoseq collapse --do-not-collapse-extra-5exons human_80K.mapped.bam human_80K.flnc.bam human_80K.collapsed.gff

多样本或者混样数据示例

一个样本不同组织间的混样用以注释物种Isoforms时，可以合并成一个样本，按照单个样本数据来注释分析。

如果是多样本，后期需要做定量或差异表达分析的话，如果没有样本名称的话（bam文件中 SM tag为样本名称)，需要在lima拆分或者去除primers的时候将每个样本重新命名，这样*collapsed.flnc_count.txt才会将每个样本的每个转录本数量给统计出来。*collapsed.flnc_count.txt文件里的矩阵，可以用于差异表达分析，如DeSeq2，的输入文件（转录本还未注释）。下面代码展示了样本重新命名的步骤：

# lima拆barcode和样本重命名
$ lima --isoseq m64307e_230628_025302.hifi_reads.bam IsoSeq_v2_primers_12.fasta UHRR.fl.bam --biosample-csv samplenames.csv --overwrite-biosample-names 
# --biosample-csv samplenames.csv --overwrite-biosample-names 要一起使用才能修改样本名称。
#samtools view -H  x.bam   检查SM tag

# Combine inputs
$ ls UHRR.fl.IsoSeqX*bam > all.fofn

# Remove poly(A) tails and concatemer
$ isoseq refine all.fofn IsoSeq_v2_primers_12.fasta UHRR.flnc.bam --require-polya

#cluster，转录本聚类
$ isoseq cluster2 UHRR.flnc.bam UHRR.transcripts.bam

#pbmm2回贴转录本，map reads using pbmm2
$ pbmm2 align --preset ISOSEQ --sort UHRR.transcripts.bam human_ref/GRCh38.fa UHRR.mapped.bam

#对回贴转录本进行collapse， collapse mapped reads into unique isoforms
$ isoseq collapse --do-not-collapse-extra-5exons UHRR.mapped.bam UHRR.flnc.bam UHRR.collapsed.gff

二、Pigeon使用方法

Pigeon是一个PacBio转录工具包，包含了用于将全长转录本isoforms按照参考基因组注释进行分类和过滤的工具。Pigeon基于SQANTI3开发，其输出与单细胞Seurat软件下游分析兼容。

官方网站：https://isoseq.how/classification/pigeon.html

• 软件安装:

$ conda install -c bioconda pbpigeon

• 准备参考基因组文件:

#使用说明
$ pigeon prepare <gencode.annotation.gtf> <reference.fa> <cage.bed> <intropolis.tsv>

#实际运行
$ pigeon prepare human_ref/GRCh38.gtf human_ref/GRCh38.fa

• 准备转录本注释（transcript GFF）文件 ，isoseq collapse 输出文件：

#使用说明
$ pigeon prepare <collapsed.gff>

#实际运行
$ pigeon prepare UHRR.collapsed.gff

Isoforms 归类 （Classify Isoforms）

• 转录本分类（Transcript classification）

#使用说明
$ pigeon classify <collapsed.sorted.gff> <annotations.gtf> <reference.fa> --cage-peak refTSS.bed --poly-a polyA.list

#实际运行
$ pigeon classify UHRR.collapsed.sorted.gff human_ref/GRCh38.sorted.gtf human_ref/GRCh38.fa

测序聚类所得的转录本根据参考基因组注释文件进行以下分类（图5），分类标准参照SQANTI3：

1. FSM (Full Splice Match)
2. ISM (Incomplete Splice Match)
3. NIC (Novel In Catalog)
4. NNC (Novel Not in Catalog)
5. Antisense
6. Genic Intron
7. Genic Genomic
8. Intergenic

图5. 转录本分类

• 将FLNC定量结果加入到分类文件中（Adding FLNC counts to classification output）

#使用说明
$ pigeon classify <collapsed.sorted.gff> <annotations.gtf> <reference.fa> --fl abundance.txt

#实际运行
$ pigeon classify UHRR.collapsed.sorted.gff human_ref/GRCh38.sorted.gtf human_ref/GRCh38.fa --fl UHRR.collapsed.flnc_count.txt

有了基因注释的表达矩阵(*_classification.txt)就可以愉快的去做差异表达分析了！

Isoforms 过滤 （Filter isoforms）

#使用说明
$ pigeon filter <classification.txt>

#实际运行
$ pigeon filter UHRR_classification.txt

如果想生成过滤后的.gff文件，需要输入 collapsed.sorted.gff。

#使用说明
$ pigeon filter <classification.txt> --isoforms <collapsed.sorted.gff>

#实际运行
$ pigeon filter UHRR_classification.txt --isoforms UHRR.collapsed.sorted.gff

参考文献

1. Transcript Classification