1. 前言
去年年初，首次完成了第一个组装的基因组的结构注释，并且把流程记录下来写成了博客，美中不足是基于二代的转录组完成的，所以今天来记录一次基于二代和三代的RNA-seq数据进行结构注释的过程。另外，本次的denovo预测从Agustus换成了ANNEVO。这个博客更多的是一个流程，一些细节不会很详细，如果看不懂，可以先看之前关于基因结构注释的博客。
2. 软件下载
git clone https://github.com/xjtu-omics/ANNEVO.git
cd ANNEVO
mmaba env create -f ANNEVO.yml -n Anno
mamba activate Anno
mamba create -n Gene_Anno
mamba install bioconda::hisat2
mamba install bioconda::trinity
mamba install bioconda::stringtie
mamba install bioconda::pasa
mamba install bioconda::samtools
mamba install bioconda::miniprot
mamba install bioconda::evidencemodeler
mamba install bioconda::minimap2
3. 结构注释
3.1 Denovo
python annotation.py \
--threads 20 \
--genome Final_genome_complete_sort.fasta \
--lineage Embryophyta \
--output Final_genome_complete_sort_denovo.gff
# --threads：CPU 数量
# --genome：基因组，不需要屏蔽
# --lineage：BUSCO/Annotation 参考数据库
# --output：输出 GFF
3.2 RNA-seq
全长转录组比对到基因组
minimap2 -ax splice \
-uf \
-t 20 \
ref.fa \
ont.filtered.fq.gz \
| samtools sort \
-@ 8 \
-o ont.sorted.bam
# -ax 比对的模式，这里是长reads比对
# -uf 设置 GT-AG 的格式，f 为标准的转录组格式
# -t 20 比对的线程数
# -@ 排序的线程数
组装三代转录本
stringtie -p 40 \
-o stringtie.gtf \
-L \
DH2-YL1.bam DH2-OL1.bam
# -p 线程数
# -o 转录本gff文件
# -L 长读长模式
# 提取转录本
gtf_genome_to_cdna_fasta.pl \
stringtie.gtf \
../Final_genome_complete_sort.fasta > stringtie.fasta
# -x 参考基因组
# -1 -2 双端测序文件
# -o 比对结果
组装二代转录本
hisat2 --new-summary /
--summary-file rnaseq.hisat.summary/
--threads 20 /
-x ../Final_genome_complete_sort.fasta /
-1 Conbine_1.fq.gz /
-2 Conbine_2.fq.gz /
| samtools sort -@ 20 /
-o rnaseq.bam -
# -x 参考基因组
# -1 -2 双端测序文件
# -o 比对结果
# 组装转录本
Trinity --max_memory 100G /
--genome_guided_bam ./rnaseq.bam /
--genome_guided_max_intron 20000 /
--CPU 10 /
--output Trinity_result
# --max_memory 最大内存
# --genome_guided_bam 输入的bam文件
# --genome_guided_max_intron 内含子最大的长度
# --CPU 运行的线程数
# --output 输出文件夹
合并二三代转录本，并且使用PASA进行转录本的组装
cat stringtie.fasta ./Trinity_result/Trinity-GG.fasta > Merge_transcript.fasta
Launch_PASA_pipeline.pl -c pasa.alignAssembly.Template.txt /
-C /
-R /
-g Final_genome_complete_sort.fasta /
-t Merge_transcript.fasta /
--ALIGNERS minimap2 /
--CPU 20
# -g 基因组文件
# -t 转录本，这里应该是2+3的
# --ALIGNERS 转录本和基因组的比对软件
# --CPU 运行的线程数
# -c 配置文件
pasa的配置文件
## templated variables to be replaced exist as <__var_name__>
# database settings
DATABASE=/home/majunpeng/sda2/Xianyu/Xigua/Xigua/Genome_asmmbly/Anno/Gene/Final_genome_complete_sort.fasta.sqlit
#######################################################
# Parameters to specify to specific scripts in pipeline
# create a key = "script_name" + ":" + "parameter"
# assign a value as done above.
#script validate_alignments_in_db.dbi
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=80
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=95
#script subcluster_builder.dbi
subcluster_builder.dbi:-m=50
2.3 EVM证据整合
到这一步，就是要整合来自denovo，同源和RNA-seq的三方注释结果。其中同源，也就是miniprot的结果需要进行格式转换。
python3 miniprot_GFF_2_EVM_alignment_GFF3.py miniprot.gff > EVM_miniprot.gff
此外还要整理一个配置文件，设定每种证据的权重
PROTEIN miniprot_protAln 3
TRANSCRIPT assembler-Bra-pasa.sqlite 10
ABINITIO_PREDICTION ANNEVO 5
证据整合
EVidenceModeler --CPU 40 \
--sample EVM \
--genome Final_genome_complete_sort.fasta \
--weights EVM_Weight.txt \
--gene_predictions Final_genome_complete_sort_denovo.gff \
--protein_alignments EVM_miniprot.gff \
--transcript_alignments Final_genome_complete_sort.fasta.sqlit.pasa_assemblies.gff3 \
--segmentSize 100000 \
--overlapSize 10000 \
--repeats ../Final_genome_complete_sort.fasta.out.gff
# --sample 输出结果前缀
# --genome 基因组
# --weights 权重文件
# --gene_predictions Dedovo预测gff文件
# --protein_alignments 同源预测gff文件
# --transcript_alignments RNA-seq预测gff文件
# --repeats 重复序列注释gff文件
# --segmentSize 基因组划分成多个片段，每个片段的大小
# --overlapSize 相邻分块间保留的重叠区域大小