美文网首页
2023-05-23 augustus训练

2023-05-23 augustus训练

作者: Athena404 | 来源:发表于2023-05-23 09:57 被阅读0次

    正在更新。英文是真的看不习惯。
    参考:使用MAKER进行基因注释(高级篇之AUGUSTUS模型训练) - 简书 (jianshu.com)
    Training AUGUSTUS (wisc.edu)
    Augustus 进行基因注释 - 斩毛毛 - 博客园 (cnblogs.com)

    首先是我曾经理解上的错误:
    1.“若存在已经被训练的物种(augustus --species=help查看,或者augustus的config/species路径下有文件夹列表),则直接使用一下代码进行预测基因,”
    说实话,我一开始理解错了。我以为是使用其他近缘物种(比如我的近缘物种是rice)也可以对我自己新测的物种参考。但是他们想说的是,如果你测序的是水稻,那你可以直接用水稻的augustus训练好的文件。
    2.“b. 这些基因的基因结构一定要足够的准确。”
    我一开始就更不理解了。我都找你来训练了,我怎么能确认我的基因结构是一定准确的?大概意思需要转录组RNA数据就可以?

    总结训练集要求:基因之间不重复(只要一个基因的一个转录本),至少随机100-200个基因,要准确。

    理论的事我就不多说了。我直接运行算了。

    我是已经用genewise进行过了同源注释。用了四个单子叶的做近缘物种,genewise.gff的BUSCO还都挺高(水稻的低一些),如果没有RNA数据用同源注释结果应该也是可以的。


    水稻genewise的BUSCO结果 同源注释的gffBUSCO结果
    同源注释的gffBUSCO结果 同源注释的gffBUSCO结果

    准备好gff之后,运行一个prepare脚本后,可以开始训练了。

    #!bin/bash
    genome=genome.fa #基因组文件
    species=Hb21 #我自己设置的物种名
    unset PERL5LIB; export PATH=/useful/perl-5.30.2/bin:$PATH
    export PATH=/share/app/blat-319/blat:$PATH
    source activate augustus
    export AUGUSTUS_CONFIG_PATH=/01.bin/Augustus/augustus_config #公共软件无法写入新文件夹,可以设置自己的路径
    
    #perl /bin/perfect_gene/perfect_gene.pl --sco 99 --start 0 --stop 1 $genome ../Oryza_sativa/Oryza_sativa.IRGSP-1.0.pep.all.fa.genewise.gff ../Brachypodium_distachyon/Brachypodium_distachyon.Brachypodium_distachyon_v3.0.pep.all.fa.genewise.gff ../Musa_acuminata/Musa_acuminata_v2.pep.all.fa.genewise.gff ../Sorghum_bicolor/Sorghum_bicolor.Sorghum_bicolor_NCBIv3.pep.all.fa.genewise.gff
    #上面这个命令应该是用同源注释做训练集。
    perl /autoAugTrain.pl --genome=$genome --species=$species --trainingset=genome.fa.gff.nr.gff --cpus 30
    

    上面这个perfect_gene.pl脚本是别人写的,我不好分享。内容大概是:(待补充)
    /autoAugTrain.pl这个好像在augustus的script文件夹下
    最后会在你设置的AUGUSTUS_CONFIG_PATH生成这些文件:

    #在AUGUSTUS_CONFIG_PATH/species/Hb21下
     Hb21_weightmatrix.txt
     Hb21_metapars.cfg
     Hb21_metapars.utr.cfg
     Hb21_metapars.cgp.cfg
     Hb21_parameters.cfg.orig1
     Hb21_parameters.cfg
     Hb21_intron_probs.pbl
     Hb21_exon_probs.pbl
     Hb21_igenic_probs.pbl
     Hb21_exon_probs.pbl.withoutCRF
     Hb21_igenic_probs.pbl.withoutCRF
     Hb21_intron_probs.pbl.withoutCRF
    
    #运行的脚本路径autoAugTrain/training下
     utr
     training.gff
     training.gb
     training.gb.train
     training.gb.test
     training.gb.onlytrain
     training.gb.train.test
     train.err
     train.out
     optimize.out
     tmp_opt_Hb21
     train.withoutCRF.err
     train.withoutCRF.out
     test
    #
    

    就是你的augustus比标准数据库里多了你训练的物种。
    然后我发现我只有训练集,没有测试集对这个进行测试。而且根据参考(avrilomics: Training the Augustus gene-finding software),他说默认会进行5次参数的优化,我这只有一次Hb21_parameters.cfg.orig1。所以我想再运行几次,再生成一个测试集(使用训练集的测试结果会虚高)。

    如果要用RNA结果训练augustus...RNA注释是有两种策略,,一种是使用HISAT2 + StringTie先比对再组装, 一种是从头组装,然后使用PASA将转录本比对到基因组上。(基因结构注释(3):转录组预测 - 简书 (jianshu.com)
    RNA下机数据经过对基因组index,

    samtools faidx genome.fasta
    bwa index genome.fasta
    java -jar picard/2.23.8/picard.jar CreateSequenceDictionary R=genome.fasta O=genome.fasta.dict
    

    去接头

    conda activate rna
    java -jar /01.Software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 16 -phred33 RNA/*85_1.fq.gz RNA/*85_2.fq.gz /cleandata/*85_1.paired.fq.gz /cleandata/*85_1.unpaired.fq.gz /cleandata/*85_2.paired.fq.gz /cleandata/*85_2.unpaired.fq.gz ILLUMINACLIP:/01.Software/Miniconda/envs/rna/adapter/adapter.fa:2:35:4:12:true  LEADING:3 TRAILING:3 SLIDINGWINDOW:5:15 MINLEN:50 2>trimming.log 
    #我的conda下环境rna安了很多处理RNA的软件。
    #路径是简化的。trimmomatic运行命令可去查找别的教程。
    #ILLUMINACLIP的fa文件和测序方联系获得。就是你接头的序列。trimmomatic软件的目录下也自带了一些接头。trimmomatic参数参考其他教程吧,不赘述。
    

    之后运行

    #!usr/bin/bash
    hisat2-build -p 4 genome.fasta genome.fasta
    hisat2 -p 4 --max-intronlen 500000 --sensitive --dta  --dta-cufflinks --phred33  --no-discordant --no-mixed  -x  genome.fasta -1 cleandata/*85_1.paired.fq.gz  -2 cleandata/*85_2.paired.fq.gz -S genome_rna.sam
    samtools view -bF 4 -S genome_rna.sam -b -o genome_rna.bam
    samtools sort genome_rna.bam  -o genome_rna.bam.sort
    stringtie genome_rna.bam.sort -p 1 -o genome_rna.transcript.gtf
    

    xzg是用的PASA,虽然我用的stringtie,但是估计也要回到PASA,因为最后要用EVM整合。我问了chatGPT有关这两种方法的区别,回答如下(警告:AI自动生成结果,不一定对,注意甄别)

    “PASA和StringTie是两个常用的生物信息学分析软件,用于转录组数据的分析和注释。它们有一些区别,包括其用途、功能和特点。
    PASA(Program to Assemble Spliced Alignments)是一个用于转录组注释和组装的软件。它的主要功能是利用已知的基因组和转录组测序数据,通过比对和组装转录本来提供基因结构注释。PASA能够将已知的基因组注释与转录组测序数据进行比对,通过多轮迭代的方法,生成准确的转录本注释,并提供详细的转录本结构信息。
    相比之下,StringTie是一个转录组装和表达定量分析工具。它可以将转录组测序数据进行拼接并组装成转录本,同时对每个转录本进行定量表达分析。StringTie具有高效的组装算法,能够有效地处理复杂的基因结构和低表达水平的转录本。它还能够生成基因和转录本的表达水平估计,提供丰富的定量信息。
    以下是PASA和StringTie之间的一些区别:
    用途:PASA主要用于转录组注释和组装,旨在提供准确的基因结构注释。StringTie主要用于转录组装和表达定量分析,旨在提供转录本的组装和表达量估计。
    功能:PASA具有注释转录本结构的功能,能够比对和组装转录本,并生成准确的注释结果。StringTie具有高效的组装算法,能够处理复杂的基因结构,并提供表达量估计和可视化分析。
    灵活性:PASA相对较为复杂,需要多轮迭代和大量的参数调整来生成准确的注释结果。StringTie具有简单的工作流程和默认参数,适用于一般的转录组装和表达量分析。
    可扩展性:PASA可以与其他转录组装软件和注释工具进行集成,如BLAT、Exonerate和BLAST等。StringTie可以与其他转录组装和表达量分析软件进行集成,如Cufflinks和Ballgown等。
    综上所述,PASA和StringTie在转录组分析中有不同的用途和功能。选择哪个软件取决于具体的研究目标和需求。如果需要准确的转录本注释和结构信息,PASA可能是一个更好的选择。如果主要关注转录组装和表达量分析,StringTie可能更适合。”

    相关文章

      网友评论

          本文标题:2023-05-23 augustus训练

          本文链接:https://www.haomeiwen.com/subject/stnzsdtx.html