使用Mikado进行基因结构注释-CSDN博客

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文件有一个为空那么就会报错，把这个组装软件在参数中删掉

参考资料

官方文档