CPAT——lncRNA编码能力预测（一）

最新推荐文章于 2024-09-12 11:05:36 发布

pingu的生信备忘录

最新推荐文章于 2024-09-12 11:05:36 发布

阅读量1k

点赞数 19

分类专栏：工具篇数据分析文章标签：数据分析

本文链接：https://blog.csdn.net/import__/article/details/141647321

版权

工具篇同时被 2 个专栏收录

5 篇文章 0 订阅

订阅专栏

数据分析

5 篇文章 0 订阅

订阅专栏

文章目录

前言
CPAT——lncRNA编码能力预测（一）

前言

对于转录组测序的数据而言，组装得到转录本之后，首先要做的就是区分蛋白编码和非蛋白编码的RNA。

目前针对这一问题，有多种解决方案，基本可以分为以下两类：

alignment-based
alignment-free

第一种算法基于序列比对，可以较好的识别保守性较好的蛋白编码基因，包括CPC ，PhyloCSF 等软件；第二种算法不需要比对，而是通过 coding 和 non-coding 转录本的序列特征来进行区分，包括 CNCI ，CPAT，PLEK等。

lncRNA 为长于 200 核苷酸的非编码蛋白质的转录物，具有链特异性，在物种间的保守性较差，另外部分 lncRNA 的染色体位置和蛋白编码基因存在重叠，通过序列比对的方式来区分容易造成误判。除此之外，基于序列比对的软件，其运行速度相对较慢，所以采用第二种算法的软件综合效果更好。

CPAT——lncRNA编码能力预测（一）

CPAT 使用 Logistic 回归模型，根据 RNA 序列四个特征预测 RNA 的编码概率：

开放阅读框大小（ORF size）
开放阅读框覆盖率（ORF coverage）
Fickett TESTCODE 统计量
六聚体使用偏差（Hexamer usage bias）

其中，前两个因素都是针对开放阅读框定义的，第一个因素是开放阅读框的大小，第二个因素是开放阅读框占转录本总长度的比例，第三个因素基于序列的碱基组成和密码子分布进行定义，第四个因素基于序列中六聚体的频率进行定义。

相较于其他基于比对的软件，如 CPC 和PhyloCSF，CPAT 具有高准确性，且其速度快了约四个数量级，使其用户能够在几秒钟内处理数千个转录本。该软件接受输入序列为 FASTA 或 bed 格式数据文件，且开发了 web 界面，允许用户提交序列并立即接收预测结果。

官网：https://cpat.readthedocs.io/en/latest/

网页版：https://wlcb.oit.uci.edu/cpat

安装

使用 pip3，安装之前需要保证 python >= 3.5，同时需要 numpy、pysam、R。

pip3 install CPAT

对于人类 human、小鼠 mouse、斑马鱼 zebrafish 和蝇类 fly，用户需要下载预建的 logit 模型和 hexamer 表；对于其他物种，可自行构建。

下载地址：https://sourceforge.net/projects/rna-cpat/files/prebuilt_models/

使用

1. 输入

用户需要提供一个基因文件（-g）、一个 logit 模型文件（-d）、一个 hexamer 频率表文件（-x），并指定输出文件名（-o）。人类 human、小鼠 mouse、斑马鱼 zebrafish 和蝇类 fly 的 logit 模型和 hexamer 频率表文件可以直接下载使用，其他物种可以参考进阶版对应内容完成构建。

基因文件可以是 BED 或 FASTA 格式，其中 BED 格式需要指定参考基因组序列文件（-r）。

FASTA 格式示例如下：

>hg19_ct_UserTrack_3545_NM_001014980 range=chr1:1177826-11821025'pad=0 3'pad=0 strand=- repeatMasking=none
CTCGCCGCGCTGAGCCGCCTCGGGACGGAGCCATGCGGCGCTGGGCCTGG
GCCGCGGTCGTGGTCCTCCTCGGGCCGCAGCTCGTGCTCCTCGGGGGCGT
CGGGGCCCGGCGGGAGGCACAGAGGACGCAGCAGCCTGGCCAGCGCGCAG
ATCCCCCCAACGCCACCGCCAGCGCGTCCTCCCGCGAGGGGCTGCCCGAG
GCCCCCAAGCCATCCCAGGCCTCAGGACCTGAGTTCTCCGACGCCCACAT
GACATGGCTGAACTTTGTCCGGCGGCCGGACGACGGCGCCTTAAGGAAGC
GGTGCGGAAGCAGGGACAAGAAGCCGCGGGATCTCTTCGGTCCCCCAGGA
CCTCCAGGTGCAGAAGTGACCGCGGAGACTCTGCTTCACGAGTTTCAGGA
GCTGCTGAAAGAGGCCACGGAGCGCCGGTTCTCAGGGCTTCTGGACCCGC
TGCTGCCCCAGGGGGCGGGCCTGCGGCTGGTGGGCGAGGCCTTTCACTGC
CGGCTGCAGGGTCCCCGCCGGGTGGAC

BED 格式示例如下：

chr1	13709021378262	NM_199121	0	+	1371128137282303299,163,3802,	0,1799,3558,
chr1	14475221470067	NM_001170535	0	+	14476481469452016331,77,102,60,70,166,70,156,57,126,125,52,71,168,109,762,	0,3869,5168,5573,6778,7998,8405,10601,11368,11696,12130,13097,14318,15552,17080,21783,
chr1	14475221470067	NM_018188	0	+	14476481469452016331,77,246,60,70,166,70,156,57,126,125,52,71,168,109,762,	0,3869,5024,5573,6778,7998,8405,10601,11368,11696,12130,13097,14318,15552,17080,21783,

2. 基础版

此处 -d 和 -x 参数对应的文件为 logit 模型和 hexamer 频率表文件。

# 输入 fasta 文件
cpat.py -g transcript.fa \
  -d dat/Human_logitModel.RData \
  -x dat/Human_Hexamer.tsv \
  -o output.txt

# 输入 bed 文件
cpat.py -r /database/hg19.fa \
  -g mRNA_hg19.bed \
  -d dat/Human_logitModel.RData \
  -x dat/Human_Hexamer.tsv \
  -o output.txt

其他可选参数：

--min-orf=MIN_ORF_LEN：最小 ORF 长度（以核苷酸为单位），默认为 75
--top-orf=N_TOP_ORF：报告的候选 ORF 数量，RNA 可能有几十个推定 ORF，在大多数情况下，真正的 ORF（按大小）被排在前几位，默认为 5
--width=LINE_WIDTH：FASTA 格式输出 ORF 的长度，默认为 100
--best-orf=MODE：选择最佳 ORF 的标准，“l” 根据 "ORF 长度 "选择，“p” 根据 "编码概率 "选择，默认为 “p”
--antisense：确定是否从反义链搜索 ORF 的逻辑值，有义链（或编码链）是指在 5′至 3′方向上携带可翻译代码的 DNA 链，默认为 False（即只搜索有义链上的 ORF）

3. 进阶版

3.1 构建 hexamer 表

make_hexamer_tab 根据 fasta 格式的 CDS 序列计算 hexamer（6mer）频率， CDS 是不含 3’ UTR 和 5’ UTR 区域的 mRNA 序列，cpat 需要此表来计算六聚体使用率得分。

make_hexamer_tab -c Human_coding_transcripts_CDS.fa.gz \
  -n Human_noncoding_transcripts_RNA.fa.gz \
  > Human_Hexamer.tsv

输出 Human_Hexamer.tsv 文件：

hexamer coding  noncoding
AAAAAA  0.00064711067360927860.001606589931772997
AAAAAC  0.000420923732220075660.0005113004850646316
AAAAAG  0.00081336231124085570.0006870944872085282
AAAAAT  0.00059172875865302710.0009504638599970318
AAAACA  0.00049346027475359820.0007256901384894673
AAAACC  0.00040038053623247950.0003686803641407804
AAAACG  9.064420497619743e-050.00010448394168197091
AAAACT  0.00040683999476466180.0004784022870680216
AAAAGA  0.00042865390390612990.000774026596998453

**补充：**此处 CDS 序列可以通过数据库下载，也可以根据基因组注释文件和序列文件使用 TBtools 工具提取。如果目标物种的基因组缺乏足够的 "编码 "和 "非编码 "基因注释来构建训练数据集，则可以考虑利用进化相关物种的数据来构建模型。

3.2 构建 logit 模型

# 输入 fasta 文件
make_logitModel -x Human_Hexamer.tsv \
  -c Human_coding_transcripts_mRNA.fa.gz \
  -n Human_noncoding_transcripts_RNA.fa.gz \
  -o Human
  
# 输入 bed 文件
make_logitModel -x Human_Hexamer.tsv \
	-c Human_coding_transcripts_hg19.bed \
	-n Human_noncoding_transcripts_hg19.bed \ 
	-r /database/hg19.fa \ 
	-o Human

3.3 检测 ORF

将蛋白编码基因序列存为 fa 文件作为输入，使用 -g 参数即可检测该蛋白编码基因的ORF，此处需要指定--antisense，否则将只搜索编码链。

cpat -x Human_Hexamer.tsv \
	-d  Human_logitModel.RData \ 
	--top-orf=100 \ 
	--antisense \
	-g test.fa \
	-o output

输出 output.ORF_prob.tsv 文件：

ID      mRNA    ORF_strand      ORF_frame       ORF_start       ORF_end ORF     Fickett Hexamer Coding_prob
NM_013387.4_ORF_1       916     -       2       327     1       327     1.103   0.28998918917275        0.792763525921043
NM_013387.4_ORF_2       916     +       2       209     430     222     1.1605  0.0674464550896935      0.271842476390681
NM_013387.4_ORF_3       916     -       1       889     695     195     0.9192  -0.32000518247443       0.0113140534730678
NM_013387.4_ORF_4       916     +       1       31      222     192     1.2952  0.600469985268255       0.915129459422605
NM_013387.4_ORF_5       916     -       1       337     197     141     1.1626  0.133867810597757       0.185245402415541
NM_013387.4_ORF_6       916     -       3       119     3       117     1.2673  0.442351820001225       0.618496534888714
NM_013387.4_ORF_7       916     -       3       842     735     108     0.5832  -0.19401829042094       0.00290794398512764
NM_013387.4_ORF_8       916     +       3       684     761     78      0.7415  -0.154613060436537      0.00454929869181486

**注意：**虽然对于大多数 mRNA 来说，最大的 ORF 通常是最有可能的，但情况并非总是如此。在 NM_013387.4 的例子中，尽管 ORF_4 不是最大的 ORF，但从其最高的编码概率来看，它最有可能编码蛋白质。

3.4 cutoff 设置

最佳 cutoff 是根据 TG-ROC 确定的。例如，在人类中，编码概率（CP）截止值 >= 0.364 表示编码序列，CP < 0.364 表示非编码序列。

species	CP threshold	Sensitivity & Specificity
Human (Homo sapiens)	0.364	0.966
Mouse (Mus musculus)	0.44	0.955
Fly (Drosophila melanogaster)	0.39	0.963
Zebrafish (Danio rerio)	0.38	0.984

最佳训练数据集显示出平衡性，即编码序列的数量与非编码序列的数量大致相等。

输入的训练 data 的格式：

Human_train.dat (20,000 rows includes 10,000 coding and 10,000 noncoding genes): 
	Column-1: gene ID
	Column-2: mRNA size
	Column-3: ORF size
	Column-4: Fickett score
	Column-5: hexamer usage
	Column-6: Label (1: coding, 0: noncoding)

R 脚本与数据可以从此处下载：https://sourceforge.net/projects/rna-cpat/files/Figure3_data/，注意训练数据应该与脚本在同一路径下

Rscript 10Fold_CrossValidation.r

性能评估采用 10 倍交叉验证，考虑的数据集包括 10,000 个编码基因和 10,000 个非编码基因。蓝色虚线表示单个 10 倍交叉验证的结果，红色实线表示 10 次验证运行的平均曲线。评估指标包括 (A) ROC 曲线，(B) 精确度-召回率 (PR) 曲线，© 精确度与截断值的关系图，以及 (D) 两条旨在确定最佳截断值的 ROC 曲线图。

_images/CPAT_performance.png

使用此处的 R 脚本可以获得 cutoff，如上图所示，最佳 cutoff 为 0.364，即 Sensitivity 与 Specificity 相交处。

4. 输出

共输出6个文件，文件名前缀为 -o 设置的部分。在这里插入图片描述

output.ORF_seqs.fa：FASTA 格式的 top ORF 序列（至少 75 个核苷酸长）
output.ORF_prob.tsv：ORF 信息（strand, frame, start, end, size, Fickett TESTCODE score, Hexamer score and coding probability）

在这里插入图片描述

output.ORF_prob.best.tsv：最佳 ORF 的信息，该文件是 output.ORF_prob.tsv` 的子集
output.no_ORF.txt：未找到 ORF 的序列 ID 或 BED 条目，应视为非编码
output1.r：Rscript 文件
CPAT_run_info.log：日志文件

参考

官方文档：https://cpat.readthedocs.io/en/latest/pdf
Wang, L., Park, H. J., Dasari, S., Wang, S., Kocher, J. P., & Li, W. (2013). CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic acids research, 41(6), e74. https://doi.org/10.1093/nar/gkt006
【生信修炼手册】CPAT:转录本蛋白编码能力预测软件：https://mp.weixin.qq.com/s/6lBLb_SzJ_eM7qUYDJ_rQg
【恒峰基因】SEQ 4. 转录本蛋白编码能力预测软件(CPAT)：https://mp.weixin.qq.com/s/pR_WijyjYy6XLil9NiSrPg
【生信森林】CPAT | 序列编码能力预测软件 | lncRNA：https://mp.weixin.qq.com/s/QTVE7zWAz6FAGzLgqQtolQ