source ./env.sh

###################################################################
cd $workdir

#链接基因组
ln -s 01.Repeat/08.RepeatMasker/genome.softmasked.fa .
ln -s 01.Repeat/08.RepeatMasker/genome.hardmasked.fa .
ln -s $datadir/genome.fasta  genome.fa


#拆分基因组不同染色体，方便后续软件使用多线程加快注释速度
##
seqkit split -p 2 genome.softmasked.fa --out-dir ref.split.softmasked
seqkit split -p 2 genome.hardmasked.fa --out-dir ref.split.hardmasked
##

mkdir 02.Braker
cd 02.Braker
#BRAKER2官方教程: https://github.com/Gaius-Augustus/BRAKER
#也就是GeneMakr, GenomeThreader, ProHint，我们还需要手动安装这些软件。但其中只有GeneMark是必须的，因为我们更多是利用RNA-seq数据进行模型训练，而GenomeThreader, ProHint是利用同源蛋白进行注释才用到，其中GenomeTrheads是处理近源蛋白，而ProHint是处理未知距离的蛋白。
#基因组命名应该简单，最好就是">contig1"或">tig000001"

#设置变量
threads=10                  # 线程数

#转录组数据
#fq1=$datadir/SRR7476717_1.fastq.gz
#fq2=$datadir/SRR7476717_2.fastq.gz
fq=$datadir/ath.rnaseq.fq.gz
est=$datadir/est.chr4.s.fa
protein=$datadir/protein.fasta
#genome=../genome.softmasked.fa
genome=../genome.fa

######################################################################
#Hist2_stringtie 比对与组装
###################################################################

# 为hisat2 构建参考序列
hisat2-build $genome $genome
# 使用 hisat2 进行比对，若果有多个样本的数据，就运行多次比对，注意修改输入的 fastq 文件和输出的 bam 文件
#双端read
for i in leaf root stem flower seed;do 
hisat2 --dta --new-summary -p $threads -x $genome -1 $datadir/${i}_1.fastq.gz -2 $datadir/${i}_2.fastq.gz  2>$i.hisat2.log| samtools sort -@ 10 > $i.hisat.bam
samtools index  $i.hisat.bam  #建立索引方便后续可视化
stringtie $i.hisat.bam   -p ${threads}   -o $i.rnaseq.gtf
done
#单端read
#hisat2 --dta --new-summary -p $threads -x $genome -U $fq   2>hisat2.log| samtools sort -@ 10 > rnaseq.bam
ls *.rnaseq.gtf >assembly_GTF_list.txt
stringtie  --merge -p $threads -o rnaseq.gtf   assembly_GTF_list.txt

gffread rnaseq.gtf -o rnaseq.gff
cat rnaseq.gff | sed 's/transcript/expressed_sequence_match/g' | \
sed 's/exon/match_part/g' | sed s/StringTie/est_gff\:est2genome/ >rnaseq_maker_input.gff
#TransDecoder提取orf
cufflinks_gtf_genome_to_cdna_fasta.pl rnaseq.gtf  ../genome.softmasked.fa > transcripts.fasta
cufflinks_gtf_to_alignment_gff3.pl rnaseq.gtf > transcripts.gff3
TransDecoder.LongOrfs -t transcripts.fasta
TransDecoder.Predict -t transcripts.fasta
cdna_alignment_orf_to_genome_orf.pl \
     transcripts.fasta.transdecoder.gff3 \
     transcripts.gff3 \
     transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3
#########################################
#运行braker  augustus 和geneMarkET 基因注释
#注意  genemark  需要key文件，
cp /share/work/biosoft/GeneMarkS/gm_key_64 ~/.gm_key
braker.pl --species=$species \
   --genome=../genome.softmasked.fa \
   --softmasking \
   --bam=stem.hisat.bam,flower.hisat.bam,root.hisat.bam,seed.hisat.bam,leaf.hisat.bam \
   --cores=$threads \
   --useexisting   \
   --gff3 --skip_fixing_broken_genes


#–fungus 基因组为真菌
#–softmasking 基因组是soft masked
#–useexisting 不从头训练，而是使用已有的物种参数
#--UTR=on 根据RNA-Seq 比对情况训练UTR参数，必须有--softmasking 和--bam