使用Mikado进行基因结构注释

简介: Mikado是基于Python3写的基因组结构注释工具,它主要做的是从多个转录组组装工具得到的转录本里挑选出最好的结果作为基因组的结构注释。此外,它还会基于同源蛋白比对结果对转录本打分。

Mikado是基于Python3写的基因组结构注释工具,它主要做的是从多个转录组组装工具得到的转录本里挑选出最好的结果作为基因组的结构注释。此外,它还会基于同源蛋白比对结果对转录本打分。换句话说这个软件主要是根据转录组数据进行注释,没有 ab inito 预测。

软件安装比较方法,我们可以使用bioconda进行安装:

conda create -n mikado mikado
# 打开Python进行测试, 注意大小写
# import Mikado
# Mikado.test()

使用Daijin准备Mikado所需文件

第一步: 准备输入

如下是下载参考序列和对应的GTF注释文件。

mkdir -p Reference
cd Reference
wget ftp://ftp.ensembl.org/pub/release-89/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.89.gtf.gz
wget ftp://ftp.ensembl.org/pub/release-89/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.dna.toplevel.fa.gz
wget "http://www.uniprot.org/uniprot/?sort=score&desc=&compress=yes&query=taxonomy:diptera%20NOT%20taxonomy:%22Drosophila%20(fruit%20flies)%20[7215]%22%20AND%20taxonomy:%22Aedes%20aegypti%22&fil=&format=fasta&force=yes" -O Aedes_aegypti.fasta.gz
gunzip *gz
cd ../

如下代码下载转录组数据

mkdir -p Reads
cd Reads
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/003/ERR1662533/ERR1662533_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/003/ERR1662533/ERR1662533_2.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/004/ERR1662534/ERR1662534_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/004/ERR1662534/ERR1662534_2.fastq.gz
cd ../

第二步:创建配置文件

使用daijin comfigure创建配置文件, 包括如下内容

  • 配置文件名: -o OUT
  • 每个任务的线程数: --threads N
  • 物种名和参考序列:--name, --genome, --transcriptome
  • 二代测序数据: --sample-sheet
  • 比对软件: -al [{gsnap,star,hisat,tophat} [{gsnap,star,hisat,tophat} ...]]
  • 组装工具: -as [{class,cufflinks,stringtie,trinity,scallop} [{class,cufflinks,stringtie,trinity,scallop} ...]]
  • 输出文件夹: -od OUT_DIR
  • 打分文件用于Mikado: --scoring
  • 蛋白数据库: --prot-db
  • 转录本的间距: --flank
  • 集群任务投递工具: --scheduler
  • 集群任务投递配置文件:-c CLUSTER_CONFIG
daijin configure --scheduler "" \
     --scoring dmelanogaster_scoring.yaml \
     --copy-scoring dmelanogaster_scoring.yaml \
     -m permissive --sample-sheet sample_sheet.tsv \
     --flank 500 -i 50 26000 --threads 2 \
     --genome Reference/Drosophila_melanogaster.BDGP6.dna.toplevel.fa \
     -al hisat -as class stringtie -od Dmelanogaster -name Dmelanogaster \
     -o daijin.yaml --prot-db Reference/Aedes_aegypti.fasta;

这里面的samples_sheet.tsv内容如下. 第一列和第二列是双端测序的read, 第三列是样本名, 第四列表示是否为链特异性建库, 包括非链特异性(fr-unstranded), 链特异性数据且第一个reads是正向链第二个reads是反向链(fr-firststrand ), 链特异性数据且第二个reads是正向链第一个reads是反向链(fr-secondstrand), 仅正向链(f)和仅反向链(r), 最后一列表示是否非为二代测序结果(False表示为二代测序)

Reads/ERR1662533_1.fastq.gz Reads/ERR1662533_2.fastq.gz     ERR1662533      fr-unstranded   False
Reads/ERR1662534_1.fastq.gz Reads/ERR1662534_2.fastq.gz     ERR1662534      fr-unstranded   False

第三步:运行

执行组装步骤。

daijin assemble --cores 20 -nd

运行时出现的问题和解决方案:

  • 对于某些服务器而言,即便在参数将任务投递系统设置为空,程序依旧会投递,解决方案就是加上-nd
  • 运行过程中会用到gnuplot进行绘图,如果报错,就找到对应行将其注释, 即下面的plot-bamstats部分。
...
rule bam_stats:
        input:
                bam=rules.bam_sort.output,
                idx=rules.bam_index.output
        output: ALIGN_DIR+"/output/{align_run}.sorted.bam.stats"
        params:
                load=loadPre(config, "samtools"),
                #plot_out=ALIGN_DIR+"/output/plots/{align_run}/{align_run}"
        threads: 1
        message: "Using samtools to collected stats for: {input}"
        shell: "{params.load} samtools stats {input.bam} > {output}"
               #" && plot-bamstats -p {params.plot_out} {output}"
...

运行结束之后得到如下文件

Dmelanogaster/3-assemblies/output/class-0-hisat-ERR1662533-0.gtf
Dmelanogaster/3-assemblies/output/class-0-hisat-ERR1662534-0.gtf
Dmelanogaster/3-assemblies/output/stringtie-0-hisat-ERR1662533-0.gtf
Dmelanogaster/3-assemblies/output/stringtie-0-hisat-ERR1662534-0.gtf

同时将总的统计结果存放在了"Dmelanogaster/3-assemblies/assembly.stats"下

第四步:运行Mikado

上一步提供了组装得到的GTF文件就可以作为Mikado的输入进行结构注释, 其中mikado要求的输入文件在dmelanogaster_scoring.yaml, 里面的内容

daijin mikado -nd Dmelanogaster/mikado.yaml

最后的结果在Dmelanogaster/5-mikado/pick/permissive/mikado-permissive.loci.gff3

如果组装结果的GTF文件有一个为空那么就会报错,把这个组装软件在参数中删掉

参考资料

目录
相关文章
|
8月前
|
JSON 数据可视化 API
Python 中调用 DeepSeek-R1 API的方法介绍,图文教程
本教程详细介绍了如何使用 Python 调用 DeepSeek 的 R1 大模型 API,适合编程新手。首先登录 DeepSeek 控制台获取 API Key,安装 Python 和 requests 库后,编写基础调用代码并运行。文末包含常见问题解答和更简单的可视化调用方法,建议收藏备用。 原文链接:[如何使用 Python 调用 DeepSeek-R1 API?](https://apifox.com/apiskills/how-to-call-the-deepseek-r1-api-using-python/)
|
8月前
|
Shell Linux
【linux】Shell脚本中basename和dirname的详细用法教程
本文详细介绍了Linux Shell脚本中 `basename`和 `dirname`命令的用法,包括去除路径信息、去除后缀、批量处理文件名和路径等。同时,通过文件备份和日志文件分离的实践应用,展示了这两个命令在实际脚本中的应用场景。希望本文能帮助您更好地理解和应用 `basename`和 `dirname`命令,提高Shell脚本编写的效率和灵活性。
624 32
|
11月前
|
数据可视化 Python
基因组之全局互作热图可视化
基因组之全局互作热图可视化
基因组之全局互作热图可视化
|
数据采集 存储 数据管理
cdga|数据治理策略:击破壁垒,迈向纵向一体化的新纪元
企业将逐步击破数据壁垒,实现数据的纵向一体化。这意味着企业能够更高效地整合内外部数据资源,形成全面、准确、及时的数据视图,为管理层提供有力的决策支持。同时,数据的一体化也将促进业务流程的优化和创新,推动企业向智能化、数字化转型迈进。
cdga|数据治理策略:击破壁垒,迈向纵向一体化的新纪元
|
Shell Linux C语言
【Shell 命令集合 文本处理工具】Linux 字段连接 join 命令使用指南
【Shell 命令集合 文本处理工具】Linux 字段连接 join 命令使用指南
231 1
|
机器学习/深度学习 存储 数据可视化
手把手教你绘制和解读实用R列线图(Nomogram):从入门到精通
手把手教你绘制和解读实用R列线图(Nomogram):从入门到精通
2966 1
|
存储 数据可视化 数据挖掘
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
知识点丨重测序数据进行kinship亲缘关系分析、构建IBS矩阵的方法与介绍
ggbiplot | 带箭头的主成分分析(PCA)图绘制
ggbiplot | 带箭头的主成分分析(PCA)图绘制
674 0
|
并行计算 算法 数据挖掘
GATK4丨如何设置线程数和内存大小?
GATK4丨如何设置线程数和内存大小?
|
数据采集 弹性计算 数据可视化
GWAS全基因组关联分析工具GAPIT最新版!详细安装教程与报错解决方案笔记(Linux版)
GWAS全基因组关联分析工具GAPIT最新版!详细安装教程与报错解决方案笔记(Linux版)