2023-05-23 augustus训练

正在更新。英文是真的看不习惯。
参考：使用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可能更适合。”