基因组结构预测

引言

基因组的结构注释目前已经有很多的流程，基本上都是结合转录组+蛋白+从头预测三方证据进行整合注释。因此，本次的结构注释选择Augustus（从头预测）+Miniprot（蛋白）+PASA&hisat2&Trinity（转录组）完成。

1. 软件下载

mamba create -n Gene_Anno

mamba activate 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::augustus

# 证据整合

mamba install bioconda::evidencemodeler

2. 软件使用

2.1 转录组部分

hisat2比对

hisat2-build genome.fasta ref

# 比对和排序

hisat2 --new-summary --summary-file rnaseq.hisat.summary --threads 10 -x ref -1 rnaseq_1.fastq.gz -2 rnaseq_2.fastq.gz | samtools sort -@ 5 -o rnaseq.bam -

samtools index rnaseq.bam

使用Trinity基于比对的bam文件进行转录组的组装

Trinity --max_memory 100G --genome_guided_bam ./rnaseq.bam --genome_guided_max_intron 20000 --CPU 10

运行PASA，第一步要修改PASA的配置文件

## templated variables to be replaced exist as <__var_name__>

# database settings
DATABASE=/tmp/Bra-pasa.sqlite # 数据库的路径，也是唯一需要硬性修改的地方，修改成你自己数据的路径

#######################################################
# 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

运行PASA pipline

### PASA pipline

Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g ./genome.fasta -t transcript.fasta --ALIGNERS minimap2 --CPU 10

#参数

# -c 配置文件

# -g 基因组文件

# -t 转录本文件

# --ALIGENERS 比对的软件

# --CPU 运行的线程数

# -C 从头创建数据库

# -R 运行组装和比对

# 输出结果

.pasa_assemblies.gff3：最重要的输出结果，最后EVM证据整合用到的输入文件

.assemblies.fasta：PASA的组装结果

# PASA本很还支持调用Transdecoder对结果进行成cds和pep的预测

/usr/local/src/PASApipeline/scripts/pasa_asmbls_to_training_set.dbi --pasa_transcripts_fasta Bra-pasa.sqlite.assemblies.fasta --pasa_transcripts_gff3 Bra-pasa.sqlite.pasa_assemblies.gff3

# --pasa_transcripts_fasta PASA的组装结果

# --pasa_transcripts_gff3 PASA的gff3格式的输出结果

# 输出结果

# .transdecoder.bed bed文件

# .transdecoder.cds 预测的cds文件

# .transdecoder.pep 预测的pep文件

# .transdecoder.gff3 带有cds的gff3文件

# .transdecoder.genome.bed 转换成基因组序列ID的bed文件

# .transdecoder.genome.gff3 转换成基因组序列ID的gff3文件 后续EVM的输入文件

2.2 蛋白部分

miniprot比对

# 将下载的多个近缘物种的蛋白序列进行合并

cat 1_pep.fasta 2_pep.fasta > total_pep.fasta

# 运行miniprot

miniprot -I -t 10 --gff --outc 0.5 --outn 100 .genome.fasta total_pep.fasta > miniprot.gff

# 参数

# -I 计算内含子最大长度的阈值

# -t 线程数

# --gff 生成gff3格式的结果

# --outc 比对的百分比阈值

# --outn 比对的输出结果阈值

# 输出结果

miniprot.gff： 比对的gff3格式的结果，后续输如到EVM

2.3 从头预测

训练模型（可选步骤，如果Augustus中有近缘物种的模型，可以直接应用）

# 训练模型

autoAugTrain.pl --genome genome.fasta --trainingset training.gff --species=Test2 --optrounds=5 --cpus=10
-v -v -v

# 参数

# --genome 基因组文件

# --trainingset 用作训练的训练集数据

# --species 生成的训练模型的名字，自己制定

# --optrounds 优化轮数

# --cpus 运行的线程数

# -v -v -v 日志的输出详细程度 3个 -v则是最详细

# 输出结果

# autoAugTrain：包含了训练过程的过程文件和最后的模型

运行Augustus（单线程）

# 直接运行Augustus

augustus --species=Test2 --genemodel=complete --softmasking=1 --codingseq=on --exonnames=on --outfile=augustus.gff --errfile=augustus.err genome.sm.fasta

# 参数

# --species 预测使用的模型

# --genemodel 只保留完整预测的基因

# --softmasking 是否使用软屏蔽后的基因序列

# --codingseq 结果中输出CDS序列

# --exonnames 结果中输出外显子名称

# --outfile 输出文件名

# --errfile 输出的报错信息文件

# 输出文件

augustus.gff：预测结果

# Augustus 还支持添加其他的证据信息进行预测（例如转录组）

# bam文件转化为Augustus的证据格式

bam2hints --in=../01.RNAseq/rnaseq.bam --out=hints.gff

# 运行Augustus

augustus --species=Test2 --hintsfile=hints.gff --extrinsicCfgFile=~/augustus_config/extrinsic/extrinsic.M.RM.E.W.cfg --genemodel=complete --softmasking=1 --codingseq=on --exonnames=on --outfile=augustus.gff --errfile=augustus.err genome.sm.fasta

# 得到gff文件后，还可以把预测的CDS和PEP序列直接提前出来

getAnnoFasta.pl --seqfile genome.sm.fasta augustus.gff

# 输出结果

# .aa PEP序列

# .codingseq CDS序列

# .cdsexons 外显子序列

多线程并行运行Augustus

# 切分基因组（按序列）

seqkit split -i -O split genome.sm.fasta

# 生成Augustus指定的序列文件列表

ls split/*.fasta | while read f ;do seqkit fx2tab -i -n -l $f |awk '{print "'$f'\t1\t" $2}'; done > chr.list

# 生成并行的Augustus运行命令

createAugustusJoblist.pl --sequences=chr.list --wrap="#" --overlap=100000 --chunksize=10000000 --outputdir=./augDir --joblist=jobs.lst --jobprefix=augSplit_ --command "augustus --species=Test1 --extrinsicCfgFile=./extrinsic.cfg --hintsfile=hints.gff --genemodel=complete --softmasking=1 --codingseq=on --exonnames=on "

# 参数

--sequences=chr.list 基因组序列文件信息

--wrap="#" 运行脚本的类型，因为我用的是服务器不是集群，所以是#

--overlap=100000 片段之间的覆盖区大小

--chunksize=10000000 拆分的片段大小

--outputdir=./augDir 输出文件夹

--joblist=jobs.lst 命令列表

--jobprefix=augSplit_ 输出的文件前缀

--command 具体的Augustus命令

# 输出结果

augSplit_X：每个并行任务的命令

jobs.lst：整合每个并行任务文件的列表文件

# 运行并行Augustus程序

parallel -j 8 "bash ./{}" < jobs.lst

#合并并行运行的结果

cat jobs.lst | xargs cat | sed -n 's/^.*--outfile=$\S*$\s.*$/\1/p' | xargs cat > augustus_concat.gtf

singularity exec ../../container/GENETools202309.sif join_aug_pred.pl < augustus_concat.gtf > augustus_combine.gtf

# 输出结果

augustus_combine.gtf：Augustus的运行结果，也是EVM的输入的从头预测的结果

2.4 EVM证据整合

预测结果文件格式调整（只需要调整Augustus和miniprot的，PASA不需要调整）

# augustus结果调整

perl /opt/EVidenceModeler/EvmUtils/misc/augustus_GTF_to_EVM_GFF3.pl augustus_combine.gtf > augustuts.gff3

# miniprot结果调整

python /opt/EVidenceModeler/EvmUtils/misc/miniprot_GFF_2_EVM_alignment_GFF3.py ../04.miniprot/miniprot.gff > miniprot.gff3

# 合并从头预测结果文件（PASA生成的结果文件然后通过Transdecoder预测的和Augustus从头预测的）

cat augustuts.gff3 pasa_transdecoder.gff3 > evm.gene_predictions.gff3

运行EVM

EVidenceModeler --CPU 10 --sample_id Bra --genome genome.fasta --weights evm_weights.txt --gene_predictions evm.gene_predictions.gff3 --protein_alignments evm.protein_alignments.gff3 --transcript_alignments evm.transcript_alignments.gff3 --repeats ../../P1.RepeatAnno/rmout_combine/genome.out.gff --segmentSize 500000 --overlapSize 50000

# 参数

--CPU 运行的线程数

--sample_id 输出文件前缀

--genome 基因组文件

--weights 权重文件

--gene_predictions 从头预测结果文件

--protein_alignments 蛋白预测结果文件

--transcript_alignments 转录组预测结果文件

--repeats 重复序列注释结果的gff文件

--segmentSize 片段大小

--overlapSize 两个连续片段的重复区域大小

# 输出结果

.EVM.pep：证据整合后预测的pep文件

.EVM.cds：证据整合后预测的cds文件

.EVM.gff3：证据整合后预测的gff3格式注释文件

.EVM.bed:证据整合后预测的bed文件

准备EVM的权重文件

# 第一列未证据名第二列是每个预测结果的第三列第三列的每个预测结果的权重

PROTEIN miniprot_protAln 2 # 蛋白预测结果
TRANSCRIPT assembler-Bra-pasa.sqlite 10 # 转录组预测结果
ABINITIO_PREDICTION Augustus 1 # 从头预测结果
OTHER_PREDICTION transdecoder 5 # 其他预测结果

2.5 注释文件中添加UTR和可变剪切

# 检查注释文件的格式

pasa_gff3_validator.pl Bra.EVM.gff3

# 加载gff3文件到数据库

Load_Current_Gene_Annotations.dbi -c alignAssembly.config -g genome.fasta -P Bra.EVM.gff3

# 更新gff3文件

Launch_PASA_pipeline.pl --CPU 10 -c annotCompare.config -A -g genome.fasta -t transcript.fasta

# 简化gff3文件

awk '$1 !~ /^#/' Bra-pasa.new.sqlite.gene_structures_post_PASA_updates.*.gff3 > Bra.EVM.update.gff3

# 提取PEP和CDS

gff3_file_to_proteins.pl Bra.EVM.update.gff3 genome.fasta prot > Bra.EVM.update.pep.fasta

gff3_file_to_proteins.pl Bra.EVM.update.gff3 genome.fasta CDS > Bra.EVM.update.cds.fasta

更改配置文件annotCompare.config

## templated variables to be replaced exist as <__var_name__>

# database settings
DATABASE=<__DATABASE__> # 更改为之前构建的数据库名称

#######################################################
# Parameters to specify to specific scripts in pipeline
# create a key = "script_name" + ":" + "parameter"
# assign a value as done above.

#script cDNA_annotation_comparer.dbi
cDNA_annotation_comparer.dbi:--MIN_PERCENT_OVERLAP=<__MIN_PERCENT_OVERLAP__>
cDNA_annotation_comparer.dbi:--MIN_PERCENT_PROT_CODING=<__MIN_PERCENT_PROT_CODING__>
cDNA_annotation_comparer.dbi:--MIN_PERID_PROT_COMPARE=<__MIN_PERID_PROT_COMPARE__>
cDNA_annotation_comparer.dbi:--MIN_PERCENT_LENGTH_FL_COMPARE=<__MIN_PERCENT_LENGTH_FL_COMPARE__>
cDNA_annotation_comparer.dbi:--MIN_PERCENT_LENGTH_NONFL_COMPARE=<__MIN_PERCENT_LENGTH_NONFL_COMPARE__>
cDNA_annotation_comparer.dbi:--MIN_FL_ORF_SIZE=<__MIN_FL_ORF_SIZE__>
cDNA_annotation_comparer.dbi:--MIN_PERCENT_ALIGN_LENGTH=<__MIN_PERCENT_ALIGN_LENGTH__>
cDNA_annotation_comparer.dbi:--MIN_PERCENT_OVERLAP_GENE_REPLACE=<__MIN_PERCENT_OVERLAP_GENE_REPLACE__>
cDNA_annotation_comparer.dbi:--STOMP_HIGH_PERCENTAGE_OVERLAPPING_GENE=<__STOMP_HIGH_PERCENTAGE_OVERLAPPING_GENE__>
cDNA_annotation_comparer.dbi:--TRUST_FL_STATUS=<__TRUST_FL_STATUS__>
cDNA_annotation_comparer.dbi:--MAX_UTR_EXONS=<__MAX_UTR_EXONS__>
cDNA_annotation_comparer.dbi:--GENETIC_CODE=<__GENETIC_CODE__>

2.6 提取最长转录组

### gff3转换成标准ensembl gtf格式

gffread -T genome.gff3 > tmp.genome.gtf

### 检查gtf文件是否有额外";",如有替换成"-"

# awk -F "\t" '$NF ~ /"\S+;\S*"/' tmp.genome.gtf |less

# mv tmp.genome.gtf tmp.gtf

# sed 's/$ "\S*$;$\S*"$/\1-\2/g' tmp.gtf > tmp.genome.gtf

# rm tmp.gtf

### 转换成标准ensembl格式

gtftk convert_ensembl -i tmp.genome.gtf > genome.gtf

## 统计

gtftk count -i genome.gtf > genome.gtf.count

## 提取编码基因gtf并统计

gtftk select_by_key -k feature -v CDS -i genome.gtf | gtftk tabulate -k transcript_id,gene_id --unique --no-header > coding_mapid.txt

perl ../script/gtf_extract.pl genome.gtf coding_mapid.txt > genes_coding.gtf

gtftk count -i genes_coding.gtf > genes_coding.gtf.count

## 筛选最长转录本

gtftk short_long -g -l -i genes_coding.gtf > genes_longest.gtf

gtftk count -i genes_longest.gtf > genes_longest.gtf.count

gffread -x genes_longest.cds.fasta -y genes_longest.pep.fasta -g genome.fasta genes_longest.gtf