该工具用于处理三代转录组(CCS)结果
具体请关注:
- https://github.com/Magdoll/cDNA_Cupcake/wiki
- https://github.com/Magdoll/cDNA_Cupcake/wiki/Cupcake:-supporting-scripts-for-Iso-Seq-after-clustering-step
如何获取高质量的转录本暂且跳过。。
1. 安装
- python 3.7
- BioPython
- bx-python
- Cupcake
## 创建一个环境
conda create -n cupcake
condo activate cupcake
## 安装所需工具
conda install -c bioconda Biopython bx-pyton cython -y
## 安装cupcake
git clone https://github.com/Magdoll/cDNA_Cupcake.git
cd cDNA_Cupcake
python setup.py build
python setup.py install
## 可以将cupcake中的脚本复制过来使用即可。
2. Collapse 冗余 isoforms (基于基因组)
在test_data中有测试数据,即
- 转录本数据hq_isoforms.fa
- ref.fa
=========================================================== - 比对(可以选用minimap2(v2.9)或者gamp都可以,这里以minimap2为例子)
## 比对
minimap2 -ax splice -t 20 -uf --secondary=no -C5 \
ref.fa hq_isoforms.fa >hq_isoforms.sam
## sort
sort -k3,3 -k4,4n hq_isoforms.sam >hq_isoforms.sorted.sam
- 利用cupcake中的脚本collapse_isoforms_by_sam.py 合并isoform
主要目的将冗余isoform以及未比对到ref上的isoform去掉。
基本用法
./collapse_isoforms_by_sam.py
usage: collapse_isoforms_by_sam.py [-h] [--input INPUT] [--fq] [-s SAM]
[-b BAM] -o PREFIX [-c MIN_ALN_COVERAGE]
[-I MIN_ALN_IDENTITY]
[--max_fuzzy_junction MAX_FUZZY_JUNCTION]
[--max_5_diff MAX_5_DIFF]
[--max_3_diff MAX_3_DIFF]
[--flnc_coverage FLNC_COVERAGE]
[--gen_mol_count] [--dun-merge-5-shorter]
[--cpus CPUS]
collapse_isoforms_by_sam.py: error: the following arguments are required: -o/--prefix
参数说明:
-i:默认0.95,相似度;
-c:默认0.99,coverage;
--max_5_diff: 5‘段外显子最大差异,默认100bp;
--max_3_diff: 3‘段外显子最大差异,默认100bp;
--dun-merge-5-shorter:不对5’端进行折叠,默认是进行折叠的
collapse_isoforms_by_sam.py --input hq_isoforms.fa \
-s hq_isoforms.sorted.sam -c 0.9 -i 0.9 \
--dun-merge-5-shorter --max_3_diff 500 \
-o hq_isoforms.sorted
结果主要得到4个文件:
(1)hq_isoforms.sorted.collapsed.gff
(2)hq_isoforms.sorted.collapsed.rep.fq
(3)hq_isoforms.sorted.collapsed.group.txt
PB.1.1 transcript/3049
PB.1.2 transcript/19325
PB.2.1 transcript/2250,transcript/2632
uniq的PB.2.1是合并了两个isoform;
PB.1.1,PB.1.2表示这个位点有两个isoforms
(4)hq_isoforms.sorted.collapsed.ids.txt (剔除的序列,由于没有比对到ref或者具有较低的coverage和identity,对此可以适当调整-c和-I参数)
3. 对唯一的isoform进行定量
获得唯一的isoform后,可以基于之前cluster步骤中得到的cluster_report.csv,可以检测比对到isoform中的FL数量。
get_abundance_post_collapse.py hq_isoforms.sorted.collapsed cluster_report.csv
共获得2个结果文件
- hq_isoforms.sorted.collapsed.read_stat.txt
id is_fl stat pbid
m64306_220512_094212/53741524/ccs Y unique PB.6356.1
m64306_220512_094212/101909358/ccs Y unique PB.6356.1
m64306_220512_094212/135922154/ccs Y unique PB.5419.1
m64306_220512_094212/26542671/ccs Y unique PB.5419.1
m64306_220512_094212/52037659/ccs Y unmapped NA
包括FL read(id)以及对应的isoform(pbid)。
- hq_isoforms.sorted.collapsed.abundance.txt
# Field explanation
# -----------------
# count_fl: Number of associated FL reads
# norm_fl: count_fl / total number of FL reads, mapped or unmapped
# Total Number of FL reads: 153557
#
pbid count_fl norm_fl
PB.1.1 2 1.3024e-05
PB.1.2 2 1.3024e-05
PB.2.1 4 2.6049e-05
该文件,第二列就是唯一的isoform所对应的FL数量。(默认唯一isoform最低有2个FL支持,当然也可以进行过滤,基于filter_by_count.py脚本,详见官方文档)。
4. 过滤掉降解的5’ isoform
基于脚本./filter_away_subset.py
在反转录获得cDNA时,容易出现5‘端降解的情况,这就导致存在多条isoform,其余序列均类似,只是5‘端降解,因此得以保留,该目的就是去除5’端降解的isoform。
例如下图中,isoforms PB10.1, PB10.2, PB10.3,PB10.4,PB10.5, PB10.6中,PB10.3,PB10.4,PB10.5, PB10.6均出现5‘端降解的情况,而经过过滤后,则只保留PB10.1, PB10.2。
filter_away_subset.py hq_isoforms.sorted.collapsed
主要获得3个结果
(1)hq_isoforms.sorted.collapsed.filtered.gff
(2)hq_isoforms.sorted.collapsed.filtered.rep.fq
(3)hq_isoforms.sorted.collapsed.filtered.abundance.txt