基因组组装和注释结果的整理 |
您所在的位置：网站首页 › 基因序列注释怎么写的啊 › 基因组组装和注释结果的整理 |
基因组组装和注释结果的整理
|
1. 基因组整理【组装之后，注释之前】1.1. 基因组碱基整理 基因组碱基-u修改成大写字母形式，同时建议使用-w 0把序列输出限定为一行序列。 seqkit seq -u -w 0 genome.old.fa>genome.new.fa 1.2. 基因组碱基小写带来的问题： edta运行TIR预测时，报错并中断运行。 1234EDTA/bin/TIR-Learner2.5/Module2/RunGRF.py", line 79, in if (len(str(records[0].seq))>int(length)+500):IndexError: list index out of rangecp: cannot stat 'TIR-Learner/*-p': No such file or directory geta生成基因注释结果减少。前期运行不受影响，会在geta生成的/out.tmp/6.combineGeneModels/中合并前面步骤预测结果时把许多基因预测结果识别为partial而筛除掉本为完整的基因，可看到/out.tmp/6.combineGeneModels/genome.gff3的gene数量是正常的，而out.tmp/6.combineGeneModels/genome.completed.gff3明显降低，/out.tmp/6.combineGeneModels/genome.partial.gff3明显升高。out.tmp/6.combineGeneModels/genome.filter.gff3是最后的基因注释结果，即out.gff3也相应减少。 maker的某些运行步骤也会出错。 1.3. 基因组序列ID修改1.3.1. 排序【optional】–根据需要选择按长度对基因组contigs进行从长到短的排序 seqkit sort -l -r -w 0 genome.old.fa >genome.sort.fa 1.3.2. 重命名【推荐】最好重命名成长度一致的数字，避免后期分析遇到contig1和contig11同时被匹配的问题。 solution A【速度最快，推荐】序列IDs按顺序命名为1-n的数字，–nr-width可以设置数字的长度，长度为5时第一条序列名称为00001。seqkit replace -p '.*' -r {nr} --nr-width 5 genome.sort.fa solution B把序列id重命名成1-n的数字awk 'BEGIN{i=0 ; FS="," ; OFS=","}{ if(/>/){gsub($1,">"++i,$1);print $0}else{print $0}}' genome.sort.fa>genome.fa 基因组的序列id改为scaf001的形式 sed -i "s/>/>scaf/g" genome.fa solution C把基因组序列id的数字1-100改成0001-0100这种数字位数相同的格式。这个没有solution A快，只推荐在不想做大规模更改的情况下使用。用下面脚本实现： 12345678#!/bin/bashj=1for i in `seq -w 1 0100` #把相同数字位数的0001-0100依次赋值给变量ido sed -i -E -e "s/scaffold$j$/scaffold_$i/g" -e "s/scaffold$j([^0-9])/scaffold_$i\1/g" species.fa #两次替换，第一次替换scaffold$j为行尾的字符串（比如在基因组序列文件中），第二次替换scaffold$j不为行尾的字符串，[^0-9]代表不为数字的任意一个字符，\1代表替换前括号([^0-9])中的内容。 sed -i -E -e "s/scaffold$j$/scaffold_$i/g" -e "s/scaffold$j([^0-9])/scaffold_$i\1/g" species.gff #同上，替换其他文件，比如gff文件。 j=$(($j+1)) #把1-100依次赋值给变量jdone 2. 基因组注释整理【注释结果整理】GFF3toolkit【推荐】推荐使用GFF3toolkit来进行gff3注释文件的整理，具体介绍参考博客GFF3toolkit blog。 GFF3toolkit包含许多模块： gff3_QC：检测gff3格式错误 gff3_fix：修正gff3格式错误 gff3_merge：合并两个gff3文件 gff3_sort：根据scaffold，coordinates坐标来排序gff3文件 gff3_to_fasta：根据基因组fasta和注释gff生成gene/cds/protein/exon等序列 2.1. 合并多个注释结果基因组的基因结构注释，如果使用了多款软件进行，想要合并多套注释结果，并让注释的基因根据染色体和位置信息排序，可参考这个办法。 合并两个基因注释文件（若是多个就依次合并） 输出A.gff中符合要求的行，即A.gff中与B.gff中位置有overlap的行。bedtools intersect -a A.gff -b B.gff -wa >A.dup.gff 输出只在A.gff存在的行，即把A.dup.gff从A.gff中删除，需要两个文件都已排序，若未排序先sort排序再处理。comm -1 A.gff A.dup.gff >A.filter.gff 合并B.gff和A.gff中与B.gff不重复的部分。cat B.gff A.filter.gff >sample.gff 2.2. 重命名注释结果【可选】根据需要选择【非必要的】，根据染色体位置顺序对注释的基因结果进行排序和重命名整理。除了基因注释外，不同分析源都可能发现新的基因，增加注释的基因数量，从而打乱顺序。 2.2.1. 排序 把sample.gff3中gene替换成1gene，mRNA替换成2mRNA，exon替换成3exon，CDS替换成4CDS；目的是确保排序后单个基因内部的顺序是1-2-3-4。sed -i -e "s/gene/1gene/g" -e "s/mRNA/2mRNA/g" -e "s/exon/3exon/g" -e "s/CDS/4CDS/g" sample.gff3 根据染色体位置排序。按照第一列第四个字符，第四列数值，第五列数值逆序依次排序。cat sample.gff3 |sort -k 1.4n -k 4n -k 5nr >sample.sort.gff3 获得的sample.sort.gff3文件再1gene替换成gene，234类似替换。sed -i -e "s/1gene/gene/g" -e "s/2mRNA/mRNA/g" -e "s/3exon/exon/g" -e "s/4CDS/CDS/g" sample.sort.gff3 2.2.2. 重命名 复制一份注释文件用于修改和替换cp sample.sort.gff3 sample.new.gff3 获取旧ID的list。获取sample.sort.gff3第九列中的ID值。cat sample.new.gff3 | awk '$3=="gene" {print $9}'|sed -e "s/;.*//g" -e "s/ID=//g" >old.name 获取新ID的list。for i in $(seq -w 1 cat old.name|wc -l); do echo scaf$i; done >new.name 合并旧的和新的ID。paste -d "/" old.name new.name > old2new.name 生成替换脚本old2new.sh。sed -e "s/^/sed -i \"s\//g" -e "s/$/\/g sample.new.gff3/g\"" old2new.name >old2new.sh 运行替换脚本，会直接替换sample.new.gff3文件的旧ID为新ID，可能会运行较长时间。sh old2new.shps：这部分代码可实现，但是是非常冗余的代码量（水平有限），且未必适用所有注释文件，谨慎参考。 2.3. 从注释文件提取序列用gffread根据注释文件从基因组提取基因序列: 从基因组genome.fa和注释文件sample.gff3获取cds序列gffread -x sample.cds.fa -g genome.fa sample.gff3 从基因组genome.fa和注释文件sample.gff3获取exon序列gffread -w sample.exon.fa -g genome.fa sample.gff3 从基因组genome.fa和注释文件sample.gff3获取protein序列gffread -y sample.protein.fa -g genome.fa sample.gff3 2.4. 从注释文件提取intron信息脚本extract_intron_info.pl 欢迎关注微信公众号：生信技工 公众号主要分享生信分析、生信软件、基因组学、转录组学、植物进化、生物学概念等相关内容，包括生物信息学工具的基本原理、操作步骤和学习心得。
|
【本文地址】
今日新闻 |
推荐新闻 |