AUGUSTUSを用いた非モデル生物ゲノムへの遺伝子アノテーション付

busco --augustus --long -m genome \
    -i genome.fa  \
    -o [output] -l [分類群] -c [CPU] 

BuildDatabase -name [database名] genome.fa
RepeatModeler -database [database名] -threads [CPU]
RepeatMasker -xsmall -dir [output] -pa [CPU] -lib RM_*/consensi.fa.classified genome.fa 

seqkit split2 genome.fa.masked -s 1 -O genome_splits

cat genome.fa.out | tail -n +3 | \
perl -ne 'chomp; s/^\s+//; @t = split(/\s+/);print $t[4]."\t"."repmask\tnonexonpart\t".$t[5]."\t".$t[6]."\t0\t.\t.\tsrc=RM\n";' | \
sort -n -k 1,1 > repeat.hints

GENOME="dir/genome.fa.masked"
GENOMESPLIT="genome_splits/"
RNAdir="RNAdir/"
READS=(DRR1 DRR2 DRR3)
SPECIES="BUSCO_Sp"
CPU=10

mkdir RNA_hints

hisat2-build -p 12 ${GENOME} genome_index
mkdir ${RNAdir}/bam

for read in ${READS[@]}; do
    hisat2 -p ${CPU} -x genome_index \
    -1 ${RNAdir}/${READS}_1.fastq.gz -2 ${RNAdir}/${READS}_2.fastq.gz \
    -S ${RNAdir}/bam/${READS}.sam
    samtools sort -@ ${CPU} ${RNAdir}/bam/${READS}.sam > ${RNAdir}/bam/${READS}.bam
    rm ${RNAdir}/bam/${READS}.sam
done

samtools merge merged.bam ${RNAdir}/bam/*bam &

samtools sort -@ ${CPU} -n merged.bam -o merged.s.bam 
filterBam --uniq --in merged.s.bam --out merged.s.f.bam 
samtools view -H merged.s.f.bam > header.txt
samtools sort -@ ${CPU} merged.s.f.bam -o both.merged.s.f.bam
bam2hints --intronsonly --in=both.merged.s.f.bam --out=hints.gff

nohup ls -1 "${GENOMESPLIT}" | xargs -n 1 -P ${CPU} -I {} \
  sh -c 'augustus "$1/$2" --species="$3" --extrinsicCfgFile=extrinsic.E.cfg \
  --alternatives-from-evidence=true --hintsfile=hints.gff \
  --allow_hinted_splicesites=atac --introns=on --genemodel=complete \
  > aug1.$2.out' _ "${GENOMESPLIT}" "{}" "${SPECIES}" &

cat aug1.{}.out > aug1.out

cat aug1.out | tee aug.prelim.gff | grep -P "\tintron\t" > aug1.introns.gff
cat hints.gff aug1.introns.gff | perl -ne '@array = split(/\t/, $_);print "$array[0]:$array[3]-$array[4]\n";'| sort -u > introns.lst
intron2exex.pl --introns=introns.lst --flank 150 --seq=${GENOME} --exex=exex.fa --map=map.psl

bowtie2-build exex.fa exex1
for read in ${READS[@]}; do
    bowtie2 -p ${CPU} -x exex1 \
    -1 ${RNAdir}/${READS}_1.fastq.gz \
    -2 ${RNAdir}/${READS}_2.fastq.gz -S ${read}.sam
done

samtools view -@ ${CPU} -S -F 4 <(cat *.sam ) > bowtie.F.sam
samMap.pl bowtie.F.sam map.psl 150 > bowtie.global.sam

bamtools filter -in merged.bam -out merged.noN.bam \
-script operation_N_filter.txt

cat header.txt bowtie.global.sam > bowtie.global.h.sam
samtools view -@ ${CPU} -bS -o bowtie.global.h.bam bowtie.global.h.sam
samtools merge -@ ${CPU} both.bam     bowtie.global.h.bam   merged.noN.bam
samtools sort -@ ${CPU} -n both.bam -o both.s.bam
filterBam --uniq --in both.s.bam --out both.sf.bam
samtools sort -@ ${CPU} both.sf.bam -o both.ssf.bam
bam2hints --intronsonly --in=both.ssf.bam --out=RNA.hints

GENOME="dir/genome.fa.masked"
CPU=10

seqkit split2 protein.faa -p 100 -O for_exonerate_faa
nohup ls -1 for_exonerate_faa/ | xargs -n 1 -P ${CPU} -I {} \
sh -c 'exonerate -q for_exonerate_faa/"$1" \
  -t "$2" --model protein2genome --querytype protein --targettype dna \
  --showvulgar no --showalignment no --showtargetgff yes --showcigar no \
  --verbose 0 --softmasktarget TRUE > "$1".gff' _ {} ${GENOME} &
  
cat *.gff > exonerate.gff
exonerate2hints.pl --in=exonerate.gff--source=P --out=exonerate.hints

GENOMESPLIT="genome_splits/"
SPECIES="BUSCO_Sp"
CPU=10

cat repeat.hints RNA.hints exonerate.hints > all.hints

nohup ls -1 ${GENOMESPLIT} | xargs -n 1 -P ${CPU}  -I {} \
sh -c 'augustus "$1/$2" --species="$3" --softmasking=1 \
--extrinsicCfgFile=extrinsic.M.RM.E.P.cfg --alternatives-from-evidence=true \
--hintsfile=all.hints --allow_hinted_splicesites=atac \
> {}.augfin.out' _ ${GENOMESPLIT} "{}" ${SPECIES} &

cat *.augfin.out > augfin.out

gffread -O augfin.out -o augfin.gff3

#ID を一意にする
awk -F'\t' '{
    OFS="\t";
    if ($3 == "transcript") {
        gsub(/ID=/, "ID="$1, $9);                     # IDに一列目 (scaffold名) を追加
    } else if ($3 == "CDS") {
        sub(/Parent=/, "Parent="$1, $9);             # Parentに一列目 (scaffold名) を追加
    }
    print $1, $2, $3, $4, $5, $6, $7, $8, $9;
}'  augfin.gff3 > augfin.uniqueID.gff3


#gffを標準化
agat_convert_sp_gxf2gxf.pl -g augfin.uniqueID.gff3 -o augfin.uniqueID.agated.gff3