使用Mikado挑选最好的转录本进行注释

本文主要是介绍使用Mikado挑选最好的转录本进行注释，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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 configure创建配置文件, 包括如下内容

• 配置文件名: -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  stringtie -od Dmelanogaster --name Dmelanogaster \-o daijin.yaml --prot-db Reference/Aedes_aegypti.fasta;

注: 在-as参数中不推荐使用CLASS2，在我的运行经验里，CLASS2的组装得到GTF文件为空，导致后续步骤会报错。也不推荐用Trinity，有一个bug导致其无法顺利运行。

这里面的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 -nd -J 20 --C 20 -t 50 -nd daijin.yaml
# -nd 一般服务器都没有DRMAA
# -J/--jobs: 每次投递多少任务
# -C/--cores: 运行多少个任务并行
# -t/--threads: 每个任务最多多少个线程，一般在daijin configure 用--threads指明

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

• 对于某些服务器而言，即便在参数将任务投递系统设置为空，程序依旧会投递，解决方案就是加上-nd
• 运行过程中会用到gnuplot进行绘图，如果报错，就找到对应行将其注释, 即下面的plot-bamstats部分。

...
rule bam_stats:input:bam=rules.bam_sort.output,idx=rules.bam_index.outputoutput: ALIGN_DIR+"/output/{align_run}.sorted.bam.stats"params:load=loadPre(config, "samtools"),#plot_out=ALIGN_DIR+"/output/plots/{align_run}/{align_run}"threads: 1message: "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 -J 10 -C 20 -nd Dmelanogaster/mikado.yaml

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

参考资料

• 官方文档

这篇关于使用Mikado挑选最好的转录本进行注释的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---