Mikado models to Gold standard models

Mikado models to Gold standard models for Augustus training

The input and output files for the commands below already exist in the directory below

cd /home/train/Annotation_workshop/Augustus/Training/Inputs/Gold_models

1. Setup the PATH

activate_reat

2. Extract CDS fasta from the Mikado GFF3 file using the reference genome

(
cat Mikado_all.loci.Chr3.gff3 \
    | gffread --stream -S -g Athaliana_447_TAIR10.Chr3.fa \
    -x proteins.fna
)

⚠️ Windows users:

cat Mikado_all.loci.Chr3.gff3 | gffread --stream -S -g Athaliana_447_TAIR10.Chr3.fa -x proteins.fna

3. Convert CDS fasta to protein fasta using appropriate translation tables/genetic codes.

(
cat proteins.fna \
    | seqkit translate -T 1 \
    | strip_stop_codon \
    > proteins.faa
)

⚠️ Windows users:

cat proteins.fna | seqkit translate -T 1 | strip_stop_codon > proteins.faa

4. Run blastp of our query protein fasta against the curated database using diamond [~2 minutes]

Download curated database and build DIAMOND database

Download curated database (Uniprot, in this case) and create diamond database for blastp of our query protein fasta

(
cd /home/train/Annotation_workshop/Augustus/Training/Inputs
gunzip uniprot_sp_tr.fasta.gz
diamond makedb \
    -d uniprot_sp_tr.fasta.dmnd \
    --in uniprot_sp_tr.fasta
gzip uniprot_sp_tr.fasta
)

⚠️ Windows users:

cd /home/train/Annotation_workshop/Augustus/Training/Inputs; gunzip uniprot_sp_tr.fasta.gz; diamond makedb -d uniprot_sp_tr.fasta.dmnd --in uniprot_sp_tr.fasta; gzip uniprot_sp_tr.fasta

(
cd /home/train/Annotation_workshop/Augustus/Training/Inputs/Gold_models
/usr/bin/time -v diamond blastp \
    -p 30 \
    -d ../uniprot_sp_tr.fasta.dmnd \
    -q proteins.faa \
    -f6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore ppos btop \
    > diamond.hits.tsv
)

⚠️ Windows users:

cd /home/train/Annotation_workshop/Augustus/Training/Inputs/Gold_models; /usr/bin/time -v diamond blastp -p 30 -d ../uniprot_sp_tr.fasta.dmnd -q proteins.faa -f6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore ppos btop > diamond.hits.tsv

5. Format Mikado GFF3 into clusters using gffread for downstream classify_transcripts command.

(
gffread \
    -g Athaliana_447_TAIR10.Chr3.fa \
    --cluster-only \
    --keep-genes \
    -P Mikado_all.loci.Chr3.gff3 \
    > all_models.clustered.gff
)

⚠️ Windows users:

gffread -g Athaliana_447_TAIR10.Chr3.fa --cluster-only --keep-genes -P Mikado_all.loci.Chr3.gff3 > all_models.clustered.gff

6. Run the classify_transcripts command to create the Gold standard models gold.gff set for Augustus training. [~1 mins]

(
/usr/bin/time -v classify_transcripts \
    --evalue_filter 1.0E-6 \
    --min_pct_cds_fraction 0.5 \
    --max_tp_utr_complete 1 \
    --max_tp_utr 2 \
    --min_tp_utr 1 \
    --max_fp_utr_complete 2 \
    --max_fp_utr 3 \
    --min_fp_utr 1 \
    --query_start_hard_filter_distance 10 \
    --query_start_score 5 \
    --query_start_scoring_distance 30 \
    --query_end_hard_filter_distance 10 \
    --query_end_score 5 \
    --query_end_scoring_distance 30 \
    --target_start_hard_filter_distance 10 \
    --target_start_score 5 \
    --target_start_scoring_distance 30 \
    --target_end_hard_filter_distance 10 \
    --target_end_score 5 \
    --target_end_scoring_distance 30 \
    --min_query_coverage_hard_filter 90 \
    --min_query_coverage_score 5 \
    --min_query_coverage_scoring_percentage 30 \
    --min_target_coverage_hard_filter 90 \
    --min_target_coverage_score 5 \
    --min_target_coverage_scoring_percentage 30 \
    --max_single_gap_hard_filter 20 \
    --max_single_gap_score 5 \
    --max_single_gap_scoring_length 30 \
    -b diamond.hits.tsv \
    -t all_models.clustered.gff
)

⚠️ Windows users:

/usr/bin/time -v classify_transcripts --evalue_filter 1.0E-6 --min_pct_cds_fraction 0.5 --max_tp_utr_complete 1 --max_tp_utr 2 --min_tp_utr 1 --max_fp_utr_complete 2 --max_fp_utr 3 --min_fp_utr 1 --query_start_hard_filter_distance 10 --query_start_score 5 --query_start_scoring_distance 30 --query_end_hard_filter_distance 10 --query_end_score 5 --query_end_scoring_distance 30 --target_start_hard_filter_distance 10 --target_start_score 5 --target_start_scoring_distance 30 --target_end_hard_filter_distance 10 --target_end_score 5 --target_end_scoring_distance 30 --min_query_coverage_hard_filter 90 --min_query_coverage_score 5 --min_query_coverage_scoring_percentage 30 --min_target_coverage_hard_filter 90 --min_target_coverage_score 5 --min_target_coverage_scoring_percentage 30 --max_single_gap_hard_filter 20 --max_single_gap_score 5 --max_single_gap_scoring_length 30 -b diamond.hits.tsv -t all_models.clustered.gff