(9)Augustus模型训练

琉璃月月呀

已于 2024-03-25 10:08:57 修改

阅读量216

点赞数 1

文章标签： python 深度学习 linux

于 2023-05-26 11:32:15 首次发布

本文链接：https://blog.csdn.net/weixin_44562189/article/details/130882912

版权

一、激活环境

二、格式转换、创建物种名、尝试训练、捕捉错误、过滤错误

三、自身比对，去除高相似度蛋白、转换格式

四、随机分组、修正密码子频率、第一次training

一、激活环境

$ conda activate funannotate 
#先激活该环境，augustus在该环境下使用 
~/anaconda3/envs/funannotate/bin 这个里面有超多脚本

二、格式转换、创建物种名、尝试训练、捕捉错误、过滤错误

$ gff2gbSmallDNA.pl chum_rename.gff3 Umbilicaria_muehlenbergii_genomic.fa 1000 Umgenes.raw.gb 
#回到我的目录下，进行格式转换，将GFF文件转换成genebank格式，左右各加1000bp序列
$ new_species.pl --species=UMfor_bad_genes_removing 
#创建一个新的物种名，由于按课本上的会提示XX文件已经存在，所以从此步骤起，将以UM作为前缀
$ etraining --species=UMfor_bad_genes_removing --stopCodonExcludedFromCDS=UMfalse genes.raw.gb 2> UMtrain.err 
#尝试训练
$ cat UMtrain.err | perl -pe 's/.*in sequence (\S+): .*/$1/' > UMbadgenes.lst 
#捕捉错误
$ filterGenes.pl UMbadgenes.lst UMgenes.raw.gb > UMgenes.gb 
#过滤掉可能错误的基因结构

三、自身比对，去除高相似度蛋白、转换格式

$ grep '/gene' UMgenes.gb |sort |uniq  |sed 's/\/gene=//g' |sed 's/\"//g' |awk '{print $1}' >UMgeneSet.lst 
#提取上一步过滤掉的UMgenes.gb中的蛋白
$ seqkit grep -f UMgeneSet_new.lst chum_use.protein.fa >UMproSet.lst.fa 
#利用seqkit提取，注意列表里的名字要和蛋白序列文件相对应，这样才能成功

$ makeblastdb -in UMproSet.lst.fa -dbtype prot -parse_seqids -out UMproSet.lst.fa 
#将得到的蛋白序列进行建库
$ blastp -db UMproSet.lst.fa -query UMproSet.lst.fa -out UMgeneSet.lst.fa.blastp -evalue 1e-5 -outfmt 6 -num_threads 8 
#自身blastp比对

#将比对结果导到本地进行查看，identity100%的自己和自己比上了，而我们需要删除其他identity >= 70%的序列
#本地Excel排序后，query序列和比对上的序列两列保留一列就行，（取前面那列）然后借助Excel去除重复项，大概有209条蛋白
#利用本地notepad++，保存要去除蛋白的名称，就像上面的那个lst一样，然后上传服务器

$ seqkit grep -f UMidentity70.lst UMproSet.lst.fa -v >UMproSet.new.lst.fa 
#利用seqkit按照要去除的identity >= 70%的蛋白list去从蛋白序列文件中去除，并取反，identity 高的 proteins 序列仅保留一条
#先在本地上准备好要剔除的序列的文件 
$ python3 um_delete.py  
#运行这个脚本，利用已经弄好的list文件，从原始的gff3里面去除identity >= 70%的序列

#该Python脚本内容！！！好好学

  1 #!/user/bin/env python    #开头必须写个这个
  2 # -*- coding:utf-8 -*-  #这个也是开头
  3 
  4 idlist = []  #建立一个列表
  5 for line in open('/ifs1/User/wuqi/Qxy/knowngenome/Umbilicaria_muehlenbergii/augustus_UM/identity70_2.lst','r'):  #对于文件里的每一行，建立一个list
    #这里用绝对路径就可以保证该脚本在哪里都能调用该文件，r 表示只读read
  6     idlist.append(line.strip())  #这个就是个指令，要写
  7 #print(idlist)  #我也不知道为啥要打印这个
  8 output = open('/ifs1/User/wuqi/Qxy/knowngenome/Umbilicaria_muehlenbergii/augustus_UM/chum_rename.new2.gff3','w') #新建一个输出文件，w表示可以写入 write
  9       
 10 for line in open('/ifs1/User/wuqi/Qxy/knowngenome/Umbilicaria_muehlenbergii/augustus_UM/chum_rename.gff3','r'): #对于打开的XX文件中的每一行
 11     conserve = True  #默认是都保存  开始做下面的事情
 12     for i in idlist:  #对于idlist里的每一个变量i
 13         if i not in line:  #如果这个i不在这个行里
 14             pass  #那么就通过
 15         else:  #否则
 16             conserve = False  #保存失败，就是不保存
 17     if conserve:  #如果保存这个选项是true（默认true嘛）
 18         output.write(line) #那么就输出这一行
 19     else:  #否则
 20         pass  #就pass掉，相当于就剔除了含有查找的该行，达到目的
 21 output.close()  #关闭

$ gff2gbSmallDNA.pl chum_rename.new2.gff3 Umbilicaria_muehlenbergii_genomic.fa 1000 Umgenes.new.gb 
#将之前费劲千辛万苦做了前边步骤，剔除一致性高达70%的序列后的gff3格式文件转为genebank格式

四、随机分组、修正密码子频率、第一次training

$ randomSplit.pl Umgenes.new.gb 100 
#将Umgenes.new.gb中的基因模型随机分成2个部分，分别放入genes.gb.test和genes.gb.train 
前者是随机抽出的100个基因，用于检测HMM参数的准确性，后者是剩下基因，进行training
$ new_species.pl --species=UM_species 
#创建初始化的HMM参数文件
$ etraining --species=UM_species Umgenes.new.gb.train > UMtrain.out 
#使用前面随机分的大块XX.gb.train文件进行training  注意名称的更改

##根据输出结果 train.out 来修正参数文件 species_parameters.cfg 中终止密码子的频率 
$ tag=$(perl -ne 'print "$1" if /tag:\s+\d+\s+\((\S+)\)/' UMtrain.out)
$ perl -p -i -e "s#/Constant/amberprob.*#/Constant/amberprob                   $tag#" /ifs1/User/wuqi/anaconda3/envs/funannotate/config/species/UM_species/UM_species_parameters.cfg
$ taa=$(perl -ne 'print "$1" if /taa:\s+\d+\s+\((\S+)\)/' UMtrain.out)
$ perl -p -i -e "s#/Constant/ochreprob.*#/Constant/ochreprob  $taa#" /ifs1/User/wuqi/anaconda3/envs/funannotate/config/species/UM_species/UM_species_parameters.cfg                 
$ tga=$(perl -ne 'print "$1" if /tga:\s+\d+\s+\((\S+)\)/' UMtrain.out)
$ perl -p -i -e "s#/Constant/opalprob.*#/Constant/opalprob  $tga#" /ifs1/User/wuqi/anaconda3/envs/funannotate/config/species/UM_species/UM_species_parameters.cfg
#这里总共有三组代码，每一组第二行都是连在一起的 注意更改train.out文件名称，XXparameters.cfg文件的地址也要找对

五、第一次检验

$ augustus --species=UM_species Umgenes.new.gb.test | tee firsttest.out 
#对初次的training，进行第一遍精确性检测，结果文件是firsttest.out
$ grep -A 22 Evaluation firsttest.out 
#该命令展示了评价部分，输出一个表格
#表格内容解读：nucleotide level，sensitivity（预测到的百分率），specificity（其中正确的百分率）
#exon level，pred total/unique（预测得到unique外显子总数），anno total/unique（实际unique外显子总数），TP（正确的预测），FP（假阳性），FN（假阴性）

###用XX.gb.train去进行training，然后用XX.gb.test去进行检测，看training的准确性（这也是最开始要随机分为两个集的目的）
#分别在核酸、外显子、基因水平检测其敏感性（sensitivity）和特异性（specificity）

六、优化参数

$ randomSplit.pl Umgenes.new.gb.train 800
#再将 training.gb.train 分成两份： 第一份 800 个基因，剩下基因为另外一份。这两份基因都用于 Optimizing meta parameters of AUGUSTUS
$ mv Umgenes.new.gb.train.train training.gb.onlytrain
使用上面的 2 份基因模型，进行 Augustus training 的优化

$ optimize_augustus.pl --species=UM_species --rounds=5 --cpus=16 --kfold=16 Umgenes.new.gb.train.test --onlytrain=training.gb.onlytrain --metapars=/ifs1/User/wuqi/anaconda3/envs/funannotate/config/species/UM_species/UM_species_metapars.cfg > optimize.out

#按如上参数，则程序会对 my_species_metapars.cfg 文件中的 28 个参数进行优化，总共优化 5 轮或有一轮找不到可优化的参数为止。每进行一个参数的优化时：将 training.gb.train.test 文件中 800 个基因随机分成 16 等份，取其中 15 等份和 training.gb.onlytrain 中的基因模型一起进行 training，剩下的 1 等份用于精确行评估要对每个等份都进行一次精确性评估；
#使用 16 个 CPU 对 16 个等份并行进行 training 和精确性评估，得到平均的精确值；优化的每个参数会分别 3~4 个值，每个值得到一个 training 的精确值，对参数的多个设定值进行比较，找到最佳的值。
#此外， training 的精确值的算法： accuracy value = (3*nucleotide_sensitivity + 2*nucleotide_specificity + 4*exon_sensitivity + 3*exon_specificity + 2*gene_sensitivity + 1*gene_specificity) / 15 。
##具体每个参数的介绍，可以看教材！

七、第二次training和第二次检验

#优化参数完毕后，需要再此使用 training.gb.train 进行一遍 training！！即进行第二次training
$ etraining --species=UM_species Umgenes.new.gb.train

#再进行第二遍的精确性评估，一般进行优化后，精确性会有较大幅度的提高
$ augustus --species=UM_species Umgenes.new.gb.test | tee secondtest.out 
#用此代码进行精确性评估
$ grep -A 22 Evaluation secondtest.out 
#用表格展示评估内容
###meta parameters调整好后必须重新Trainning模型，否则没有任何意义

八、进行CRF 训练

#进行CRF Training
#在进行 Training with CRF 之前，先备份非 CRF 的参数文件，把3个XX.pbl格式的文件进行备份
$ cd ~/anaconda3/envs/funannotate/config/species/UM_species
$ cp UM_species_exon_probs.pbl UM_species_exon_probs.pbl.withoutCRF
$ cp UM_species_igenic_probs.pbl UM_species_igenic_probs.pbl.withoutCRF
$ cp UM_species_intron_probs.pbl UM_species_intron_probs.pbl.withoutCRF
$ cd -  #转到上一次访问的目录，如果经常在两个目录之间切换，可以用这样，方便一些

#在 training 的时候，加入 --CRF=1 的参数，这样进行 training 比较耗时
$ etraining --species=UM_species Umgenes.new.gb.train --CRF=1 1>train.CRF.out 2>train.CRF.err
$ etraining --species=UM_species --CRF=1 Umgenes.new.gb.train
#再次进行精确行评估
$ augustus --species=UM_species Umgenes.new.gb.test | tee secondtest.out.withCRF
#比较 CRF 和 非 CRF 两种情况下的精确性。一般情况下，CRF training 的精确性要高些。若 CRF training 的精确行低些，则将备份的参数文件还原回去即可。