Well, here we go! Initiating full-blown annotation of the Pgenerosa_v070 (FastA; 2.1GB), using MAKER (v.2.31.10) on Mox. This will perform the following:
- one round of MAKER gene model predictions
- two rounds of SNAP gene model training/predictions
- renaming of gene models to NCBI-standardized convention
- functional characterization of protein models (via BLASTp)
- functional characterization of protein domains (via InterProScan5)
I’ve submitted the job to Mox and now we wait. When I last ran this process on the Olympia oyster genome, the job took about two weeks to complete. I’ve made some changes that will allow it to run a bit faster (I think/hope) - primarily by using gene model GFFs for the SNAP gene model trainings, instead of relying on the FastA file file each time, since the FastA file needs to get BLASTed on each round and BLASTing is a time consuming process.
To finish off this post, I’ll provide the SBATCH script used to submit this job to Mox. Be prepared, it’s a doozy and it took me many hours to fully put this together (the script has >350 lines in it!). When the job summary is typed out like above, it seems so easy! Yeesh!
#!/bin/bash
## Job Name
#SBATCH --job-name=maker
## Allocation Definition
#SBATCH --account=srlab
#SBATCH --partition=srlab
## Resources
## Nodes
#SBATCH --nodes=2
## Walltime (days-hours:minutes:seconds format)
#SBATCH --time=25-00:00:00
## Memory per node
#SBATCH --mem=120G
##turn on e-mail notification
#SBATCH --mail-type=ALL
#SBATCH --mail-user=samwhite@uw.edu
## Specify the working directory for this job
#SBATCH --workdir=/gscratch/scrubbed/samwhite/outputs/20190115_geoduck_maker_genome_annotation
# Load Python Mox module for Python module availability
module load intel-python3_2017
# Load Open MPI module for parallel, multi-node processing
module load icc_19-ompi_3.1.2
# SegFault fix?
export THREADS_DAEMON_MODEL=1
# Document programs in PATH (primarily for program version ID)
date >> system_path.log
echo "" >> system_path.log
echo "System PATH for $SLURM_JOB_ID" >> system_path.log
echo "" >> system_path.log
printf "%0.s-" {1..10} >> system_path.log
echo ${PATH} | tr : \\n >> system_path.log
# Add BLAST to system PATH
export PATH=$PATH:/gscratch/srlab/programs/ncbi-blast-2.6.0+/bin
export BLASTDB=/gscratch/srlab/blastdbs/UniProtKB_20181008/
## Establish variables for more readable code
wd=$(pwd)
maker_dir=/gscratch/srlab/programs/maker-2.31.10/bin
snap_dir=/gscratch/srlab/programs/maker-2.31.10/exe/snap
### Paths to Maker binaries
maker=${maker_dir}/maker
gff3_merge=${maker_dir}/gff3_merge
maker2zff=${maker_dir}/maker2zff
fathom=${snap_dir}/fathom
forge=${snap_dir}/forge
hmmassembler=${snap_dir}/hmm-assembler.pl
fasta_merge=${maker_dir}/fasta_merge
map_ids=${maker_dir}/maker_map_ids
map_gff_ids=${maker_dir}/map_gff_ids
map_fasta_ids=${maker_dir}/map_fasta_ids
functional_fasta=${maker_dir}/maker_functional_fasta
functional_gff=${maker_dir}/maker_functional_gff
ipr_update_gff=${maker_dir}/ipr_update_gff
iprscan2gff3=${maker_dir}/iprscan2gff3
blastp_dir=${wd}/blastp_annotation
maker_blastp=${wd}/blastp_annotation/20190108_blastp.outfmt6
maker_prot_fasta=${wd}/snap02/Pgenerosa_v70_snap02.all.maker.proteins.fasta
maker_prot_fasta_renamed=${wd}/snap02/Pgenerosa_v70_snap02.all.maker.proteins.renamed.fasta
maker_transcripts_fasta=${wd}/snap02/Pgenerosa_v70_snap02.all.maker.transcripts.fasta
maker_transcripts_fasta_renamed=${wd}/snap02/Pgenerosa_v70_snap02.all.maker.transcripts.renamed.fasta
snap02_est_gff=${wd}/snap02/Pgenerosa_v70_snap01.maker.all.noseqs.est2genome.gff
snap02_gff=${wd}/snap02/Pgenerosa_v70_snap02.all.gff
snap02_gff_renamed=${wd}/snap02/Pgenerosa_v70_snap02.all.renamed.gff
snap02_protein_gff=${wd}/snap02/Pgenerosa_v70_snap01.maker.all.noseqs.protein2genome.gff
snap02_rm_gff=${wd}/snap02/Pgenerosa_v70_snap01.maker.all.noseqs.repeats.gff
put_func_gff=Pgenerosa_v70_genome_snap02.all.renamed.putative_function.gff
put_func_prot=Pgenerosa_v70_genome_snap02.all.maker.proteins.renamed.putative_function.fasta
put_func_trans=Pgenerosa_v70_genome_snap02.all.maker.transcripts.renamed.putative_function.fasta
put_domain_gff=Pgenerosa_v70_genome_snap02.all.renamed.putative_function.domain_added.gff
ips_dir=${wd}/interproscan_annotation
ips_base=Pgenerosa_v70_maker_proteins_ips
ips_name=Pgenerosa_v70_maker_proteins_ips.tsv
id_map=${wd}/snap02/Pgenerosa_v70_genome.map
ips_domains=Pgenerosa_v70_genome_snap02.all.renamed.visible_ips_domains.gff
## Path to blastp
blastp=/gscratch/srlab/programs/ncbi-blast-2.6.0+/bin/blastp
## Path to InterProScan5
interproscan=/gscratch/srlab/programs/interproscan-5.31-70.0/interproscan.sh
## Store path to options control file
maker_opts_file=./maker_opts.ctl
### Path to genome FastA file
genome=/gscratch/srlab/sam/data/P_generosa/generosa_genomes/Pgenerosa_v070.fa
### Path to transcriptome FastA file
transcriptome=/gscratch/srlab/sam/data/P_generosa/generosa_transcriptomes/20180827_trinity_geoduck.fasta
### Path to Crassotrea gigas NCBI protein FastA
gigas_proteome=/gscratch/srlab/sam/data/C_gigas/gigas_ncbi_protein/GCA_000297895.1_oyster_v9_protein.faa
### Path to Crassostrea virginica NCBI protein FastA
virginica_proteome=/gscratch/srlab/sam/data/C_virginica/virginica_ncbi_protein/GCF_002022765.2_C_virginica-3.0_protein.faa
### Path to Panopea generosa TransDecoder protein FastA
panopea_td_proteome=/gscratch/srlab/sam/data/P_generosa/generosa_proteomes/20180827_trinity_geoduck.fasta.transdecoder.pep
### Path to concatenated NCBI prteins FastA
combined_proteomes=${wd}/combined_proteomes.fasta
### Path to P.generosa-specific repeat library
repeat_library=/gscratch/srlab/sam/data/P_generosa/generosa_repeats/Pgenerosa_v070-families.fa
### Path to SwissProt database for BLASTp
sp_db_blastp=/gscratch/srlab/blastdbs/UniProtKB_20190109/uniprot_sprot.fasta
## Make directories
mkdir blastp_annotation
mkdir interproscan_annotation
## Create Maker control files needed for running Maker, only if it doesn't already exist and then edit it.
### Edit options file
### Set paths to P.generosa genome and transcriptome.
### Set path to combined C. gigas, C.virginica, P.generosa proteomes.
### The use of the % symbol sets the delimiter sed uses for arguments.
### Normally, the delimiter that most examples use is a slash "/".
### But, we need to expand the variables into a full path with slashes, which screws up sed.
### Thus, the use of % symbol instead (it could be any character that is NOT present in the expanded variable; doesn't have to be "%").
if [ ! -e maker_opts.ctl ]; then
$maker -CTL
sed -i "/^genome=/ s% %$genome %" "$maker_opts_file"
sed -i "/^est=/ s% %$transcriptome %" "$maker_opts_file"
sed -i "/^protein=/ s% %$combined_proteomes %" "$maker_opts_file"
sed -i "/^rmlib=/ s% %$repeat_library %" "$maker_opts_file"
sed -i "/^est2genome=0/ s/est2genome=0/est2genome=1/" "$maker_opts_file"
sed -i "/^protein2genome=0/ s/protein2genome=0/protein2genome=1/" "$maker_opts_file"
fi
## Create combined proteome FastA file, only if it doesn't already exist.
if [ ! -e combined_proteomes.fasta ]; then
touch combined_proteomes.fasta
cat "$gigas_proteome" >> combined_proteomes.fasta
cat "$virginica_proteome" >> combined_proteomes.fasta
cat "$panopea_td_proteome" >> combined_proteomes.fasta
fi
## Run Maker
### Specify number of nodes to use.
mpiexec -n 56 $maker
## Merge gffs
${gff3_merge} -d Pgenerosa_v70.maker.output/Pgenerosa_v70_master_datastore_index.log
## GFF with no FastA in footer
${gff3_merge} -n -s -d Pgenerosa_v70.maker.output/Pgenerosa_v70_master_datastore_index.log > Pgenerosa_v70.maker.all.noseqs.gff
## Merge all FastAs
${fasta_merge} -d Pgenerosa_v70.maker.output/Pgenerosa_v70_master_datastore_index.log
## Extract GFF alignments for use in subsequent MAKER rounds
### Transcript alignments
awk '{ if ($2 == "est2genome") print $0 }' Pgenerosa_v70.maker.all.noseqs.gff > Pgenerosa_v70.maker.all.noseqs.est2genome.gff
### Protein alignments
awk '{ if ($2 == "protein2genome") print $0 }' Pgenerosa_v70.maker.all.noseqs.gff > Pgenerosa_v70.maker.all.noseqs.protein2genome.gff
### Repeat alignments
awk '{ if ($2 ~ "repeat") print $0 }' Pgenerosa_v70.maker.all.noseqs.gff > Pgenerosa_v70.maker.all.noseqs.repeats.gff
## Run SNAP training, round 1
mkdir snap01 && cd snap01
${maker2zff} ../Pgenerosa_v70.all.gff
${fathom} -categorize 1000 genome.ann genome.dna
${fathom} -export 1000 -plus uni.ann uni.dna
${forge} export.ann export.dna
${hmmassembler} genome . > Pgenerosa_v70_snap01.hmm
## Initiate second Maker run.
### Copy initial maker control files and
### - change gene prediction settings to 0 (i.e. don't generate Maker gene predictions)
### - use GFF subsets generated in first round of MAKER
### - set location of snaphmm file to use for gene prediction
### Percent symbols used below are the sed delimiters, instead of the default "/",
### due to the need to use file paths.
if [ ! -e maker_opts.ctl ]; then
$maker -CTL
sed -i "/^genome=/ s% %$genome %" maker_opts.ctl
sed -i "/^est2genome=1/ s/est2genome=1/est2genome=0/" maker_opts.ctl
sed -i "/^protein2genome=1/ s/protein2genome=1/protein2genome=0/" maker_opts.ctl
sed -i "/^est_gff=/ s% %../Pgenerosa_v70.maker.all.noseqs.est2genome.gff %" maker_opts.ctl
sed -i "/^protein_gff=/ s% %../Pgenerosa_v70.maker.all.noseqs.protein2genome.gff %" maker_opts.ctl
sed -i "/^rm_gff=/ s% %../Pgenerosa_v70.maker.all.noseqs.repeats.gff %" maker_opts.ctl
sed -i "/^snaphmm=/ s% %Pgenerosa_v70_snap01.hmm %" maker_opts.ctl
fi
## Run Maker
### Set basename of files and specify number of CPUs to use
mpiexec -n 56 $maker \
-base Pgenerosa_v70_snap01
## Merge gffs
${gff3_merge} -d Pgenerosa_v70_snap01.maker.output/Pgenerosa_v70_snap01_master_datastore_index.log
## GFF with no FastA in footer
${gff3_merge} -n -s -d Pgenerosa_v70_snap01.maker.output/Pgenerosa_v70_snap01_master_datastore_index.log > Pgenerosa_v70_snap01.maker.all.noseqs.gff
## Extract GFF alignments for use in subsequent MAKER rounds
### Transcript alignments
awk '{ if ($2 == "est2genome") print $0 }' Pgenerosa_v70_snap01.maker.all.noseqs.gff > Pgenerosa_v70_snap01.maker.all.noseqs.est2genome.gff
### Protein alignments
awk '{ if ($2 == "protein2genome") print $0 }' Pgenerosa_v70_snap01.maker.all.noseqs.gff > Pgenerosa_v70_snap01.maker.all.noseqs.protein2genome.gff
### Repeat alignments
awk '{ if ($2 ~ "repeat") print $0 }' Pgenerosa_v70_snap01.maker.all.noseqs.gff > Pgenerosa_v70_snap01.maker.all.noseqs.repeats.gff
## Run SNAP training, round 2
cd ..
mkdir snap02 && cd snap02
${maker2zff} ../snap01/Pgenerosa_v70_snap01.all.gff
${fathom} -categorize 1000 genome.ann genome.dna
${fathom} -export 1000 -plus uni.ann uni.dna
${forge} export.ann export.dna
${hmmassembler} genome . > Pgenerosa_v70_snap02.hmm
## Initiate third and final Maker run.
### Copy initial maker control files and:
### - change gene prediction settings to 0 (i.e. don't generate Maker gene predictions)
### - use GFF subsets generated in first round of SNAP
### - set location of snaphmm file to use for gene prediction.
### Percent symbols used below are the sed delimiters, instead of the default "/",
### due to the need to use file paths.
if [ ! -e maker_opts.ctl ]; then
$maker -CTL
sed -i "/^genome=/ s% %$genome %" maker_opts.ctl
sed -i "/^est2genome=1/ s/est2genome=1/est2genome=0/" maker_opts.ctl
sed -i "/^protein2genome=1/ s/protein2genome=1/protein2genome=0/" maker_opts.ctl
sed -i "/^est_gff=/ s% %${snap02_est_gff} %" maker_opts.ctl
sed -i "/^protein_gff=/ s% %${snap02_protein_gff} %" maker_opts.ctl
sed -i "/^rm_gff=/ s% %${snap02_rm_gff} %" maker_opts.ctl
sed -i "/^snaphmm=/ s% %Pgenerosa_v70_snap02.hmm %" maker_opts.ctl
fi
## Run Maker
### Set basename of files and specify number of CPUs to use
mpiexec -n 56 $maker \
-base Pgenerosa_v70_snap02
## Merge gffs
${gff3_merge} \
-d Pgenerosa_v70_snap02.maker.output/Pgenerosa_v70_snap02_master_datastore_index.log
## GFF with no FastA in footer
${gff3_merge} -n -s -d Pgenerosa_v70_snap02.maker.output/Pgenerosa_v70_snap02_master_datastore_index.log > Pgenerosa_v70_snap02.maker.all.noseqs.gff
## Merge FastAs
${fasta_merge} \
-d Pgenerosa_v70_snap02.maker.output/Pgenerosa_v70_snap02_master_datastore_index.log
# Create copies of files for mapping
cp ${maker_prot_fasta} ${maker_prot_fasta_renamed}
cp ${maker_transcripts_fasta} ${maker_transcripts_fasta_renamed}
cp ${snap02_gff} ${snap02_gff_renamed}
# Map IDs
## Change gene names
${maker_map_ids} \
--prefix PGEN_ \
--justify 8 \
${snap02_gff} \
> ${id_map}
## Map GFF IDs
${map_gff_ids} \
${id_map} \
${snap02_gff_renamed}
## Map FastAs
### Proteins
${map_fasta_ids} \
${id_map} \
${maker_prot_fasta_renamed}
### Transcripts
${map_fasta_ids} \
${id_map} \
${maker_transcripts_fasta_renamed}
# Run InterProScan 5
## disable-precalc since this requires external database access (which Mox does not allow)
cd ${ips_dir}
${interproscan} \
--input ${maker_prot_fasta_renamed} \
--goterms \
--output-file-base ${ips_base} \
--disable-precalc
# Run BLASTp
cd ${blastp_annotation}
${blastp} \
-query ${maker_prot_fasta_renamed} \
-db ${sp_db_blastp} \
-out ${maker_blastp} \
-max_target_seqs 1 \
-evalue 1e-6 \
-outfmt 6 \
-num_threads 28
# Functional annotations
cd ${wd}
## Add putative gene functions
### GFF
${functional_gff} \
${sp_db_annotations} \
${maker_blastp} \
${snap02_gff_renamed} \
> ${put_func_gff}
### Proteins
${maker_functional_fasta} \
${sp_db_annotations} \
${maker_blastp} \
${maker_prot_fasta_renamed} \
> ${put_func_prot}
### Transcripts
${maker_functional_fasta} \
${sp_db_annotations} \
${maker_blastp} \
${maker_transcripts_fasta_renamed} \
> ${put_func_trans}
## Add InterProScan domain info
### Add searchable tags
${ipr_update_gff} \
${put_func_gff} \
${ips_dir}/${ips_name} \
> ${put_domain_gff}
### Add viewable features for genome browsers (JBrowse, Gbrowse, Web Apollo)
${iprscan2gff3} \
${ips_dir}/${ips_name} \
${snap02_gff_renamed} \
> ${ips_domains}