Genome annotation with MAKER2
Update: I began the annotation process with MAKER2 and took it through a preliminary annotation, but did not continue with multiple rounds of reannotation to produce a stable, well-support set of annotations because FUNannotate was easier to install and run. The analyses are similar, but FUNannotate automates many steps (e.g. transcriptome processing with PASA), more easily incorporates evidence sources, and is more amenable to parallelization and is therefore much faster. I recommend using FUNannotate instead of MAKER2, and will not follow up on the analyses here. Furthermore, I made the mistake of providing protein evidence from too distantly related species to MAKER2 (see section 3.1 and FUNannotate section 4), so continuing with MAKER2 would be restarting the annotation process.
[isg4@farnam ~]$ conda create --name genome
[isg4@farnam ~]$ conda activate genomeI'm going to try a suite of somewhat redundant quality control packages, most of which focus metrics like contig and gap size and Benchmarking Universal Single Copy Orthologs (BUSCOs). Many of these same metrics are remeasured in later annotation and prediction steps (e.g. see Training Augustus), but it's useful to recalculate these (in particular BUSCO stats) as our annotation improves. Furthermore, some of these pipelines have installation processes that are major headaches (for me, circos was a tough one), while others are simple webservers like gVolante. So it's good to know if they produce similar results, when possible, if you don't want to go through superfluous installs and configurations.
0.1 QUAST
QUAST computes various genome quality metrics about
- Contig sizes
- Misassemblies and structural variations
- Genome representation and its functional elements
- Variations of N50 based on aligned blocks
and generates plots for some of these. It can be used with and without a reference genome so that initial quality can be assessed and reassessed after subsequent rounds of curation. One of the most interesting features of QUAST is the Circos output for mapping and visualizing multiple genomic features to a circular representation of the genome. It has to be installed separately, but QUAST will run without it.
QUAST Installation Version: 5.0.2
[isg4@farnam apps]$ wget https://downloads.sourceforge.net/project/quast/quast-5.0.2.tar.gz
[isg4@farnam apps]$ tar -xzf quast-5.0.2.tar.gz
[isg4@farnam quast-5.0.2]$ cd quast-5.0.2
[isg4@farnam quast-5.0.2]$ ./setup.py install_full
[isg4@farnam quast-5.0.2]$ ./setup.py testCircos Installation Version: 0.69-6
[isg4@farnam apps]$ wget http://www.circos.ca/distribution/circos-0.69-6.tgz
[isg4@farnam apps]$ tar -xvzf circos-0.69-6.tgz
# Load dependency libgd
[isg4@farnam apps]$ module load libgd/2.2.3-foss-2016b
[isg4@farnam apps]$ cd /home/isg4/project/apps/circos-0.69-6/bin
[isg4@farnam bin]$ circos -modulesThe final command should list the perl dependencies and their status. I installed everything that was missing with cpanm. Some modules (Clone, GD, GD::Polyline, and List::MoreUtils) are throwing the following error:
/bin/sh: /tmp/build/80754af9/perl_1527832170752/_build_env/bin/x86_64-conda_cos6-linux-gnu-gcc: No such file or directoryAfter digging around different forums I finally found a solution; paths to some compilers seemed to be hardcoded, so they were being searched for in places that didn't exist. Switching to a different perl build fixed this issue.
[isg4@farnam bin]$ conda install -c conda-forge perl=5.26.2=h470a237_0After switching builds, only GD and GD::Polyline remain--which the Circos authors also state are troublesome. I've found that this walk-through on locally compiling is one of the more easy to follow if you're having trouble installing packages without root privledges.
# Install libjpeg
[isg4@farnam ~]$ cd ~/project/apps
[isg4@farnam apps]$ curl -O http://www.ijg.org/files/jpegsrc.v8d.tar.gz
[isg4@farnam apps]$ tar -xzvf jpegsrc.v8d.tar.gz
[isg4@farnam apps]$ cd jpeg-8d/
[isg4@farnam jpeg-8d]$ export HOME=/home/isg4
[isg4@farnam jpeg-8d]$ export PATH=$HOME/bin:$PATH
[isg4@farnam jpeg-8d]$ export LD_LIBRARY_PATH=$HOME/lib:$LD_LIBRARY_PATH
[isg4@farnam jpeg-8d]$ ./configure --prefix=$HOME
[isg4@farnam jpeg-8d]$ make
[isg4@farnam jpeg-8d]$ make install
# Install libpng
[isg4@farnam apps]$ curl -O https://sourceforge.net/projects/libpng/files/libpng16/1.6.36/libpng-1.6.36.tar.gz
[isg4@farnam apps]$ tar -xzvf libpng-1.6.36.tar.gz
[isg4@farnam apps]$ cd libpng-1.6.36/
[isg4@farnam libpng-1.6.36]$ ./configure --prefix=$HOME
[isg4@farnam libpng-1.6.36]$ make
[isg4@farnam libpng-1.6.36]$ make install
# Install lidgd
[isg4@farnam apps]$ curl -O https://github.com/libgd/libgd/releases/download/gd-2.2.5/libgd-2.2.5.tar.gz
[isg4@farnam apps]$ tar -xzvf libgd-2.2.g.tar.gz
[isg4@farnam apps]$ cd libgd-2.2.5/
[isg4@farnam libgd-2.2.5]$ ./configure --prefix=$HOME
** Configuration summary for libgd 2.2.5:
Support for Zlib: yes
Support for PNG library: no
Support for JPEG library: no
Support for WebP library: no
Support for TIFF library: no
Support for Freetype 2.x library: no
Support for Fontconfig library: no
Support for Xpm library: no
Support for liq library: no
Support for pthreads: yesThe circos install guide says we'll need support for Freetype and Fontconfig. But running circos -module showed that everything looks fine. I'll try running the example dataset as a test.
[isg4@farnam ~]$ cd ~/project/apps/circos-0.69-6/example/
[isg4@farnam example]$ ./runThis should create circos.png (below) and circos.svg if everything is in working order.
QUAST command format:
quast.py [input_sequence_data] -t [cores]Now to run QUAST with Circos on our own data. Based on the type of data you start with, you may want to add some additional parameters to your QUAST command. For example,
-
--split-scaffolds: assemblies are scaffolds (rather than contigs) -
--eukaryote: genome is eukaryotic–affects gene finding, conserved orthologs finding and contig alignment -
--large: genome is large (typically > 100 Mbp)–affects some parameters -
--k-mer-stats: compute k-mer-based quality metrics -
--circos: plot Circos version of Icarus contig alignment viewer -
--gene-finding: enables gene finding -
--rna-finding: enables ribosomal RNA gene finding -
--conserved-genes-finding: enables search for Universal Single-Copy Orthologs usingBUSCO
Some of these commands may be redundant (e.g. I think that --large is forced by --eukaryote), but I'm not entirely sure based on the documentation. When I first ran QUAST I received the following warning
WARNING: License period for GeneMark has ended! To update license, please visit http://exon.gatech.edu/GeneMark/license_download.cgi page and fill in the form. You should choose GeneMarkS tool and your operating system (note that GeneMark is free for non-commercial use). Download the license key and replace your ~/.gm_key with the updated version. After that you can restart QUAST.
so I went over to the GeneMark download page and grabbed the GeneMark-ES / ET v.4.38 for Linux 64 key (gm_key_64), which fixed the problem. However, if you're using the conda build of QUAST see issue 84 for a work around. My final QUAST command is rather expensive because of the prediction steps so, after I got everything working in an interactive bash session, I submitted it as a batch file.
#!/bin/bash
#SBATCH --job-name=QUAST
#SBATCH -c 20
#SBATCH --mem=36G
#SBATCH -t 120:00:00
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH -o quast.%A-%a.out
#SBATCH -e quast.%A-%a.err
python /home/isg4/project/apps/quast-5.0.2/quast.py ~/scratch60/P_amilis_genome/Data/portulaca_26Nov2018_oK3Ko.fixed.fasta -t 20 --split-scaffolds --eukaryote --large --k-mer-stats --circos --gene-finding --rna-finding --conserved-genes-findingThe Circos plot isn't being generated here for some reason. I'm thinking because I'm not specifying a reference genome.
- Get
Circosplot working
0.2 gVolante: Completeness Assessment of Genome/Transcriptome Sequences
gVolante let's users assess the completeness of their genomes using contig and BUSCO metrics. It's an easy to use webserver. Here is a link to my results, and below is a summary of the BUSCO results.
| Total number of core genes queried | 1440 |
|---|---|
| Complete core genes detected | 1291 (89.65%) |
| Complete + Partial core genes detected | 1320 (91.67%) |
| Complete and duplicated | 63 (4.375%) |
| Number of missing core genes | 120 (8.33%) |
| Average number of orthologs per core genes | 1.06 |
After the run finished I downloaded a compressed file containing to results because I'm not sure how long the link will last.
The functional annotation of a genome requires a set of expressed sequence tags (ESTs), which usually come in the form of an assembled and annotated transcriptome. I will be following Ya Yang's pipeline in part, but using the reference-guided version of Trinity to build my transcriptome.
Fortunately for me, Portulaca amilis has a leaf transcriptome in the 1KP Initiative (reference ID: ERR2040261), which I accessed using fastaq-dump (Version: 2.9.4) from the SRA-Toolkit.
[isg4@farnam ~]$ cd /home/isg4/scratch60/P_amilis_genome/Data
[isg4@farnam Data]$ module load SRA-Toolkit
[isg4@farnam Data]$ fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files ERR2040261The --defline-seq flag is needed to reformat the headers of the forward and reverse files so that downstream programs recognize them, and --split-files tells fastq-dump to output the forward and reverse reads into separate files. Now we'll clone Dr. Yang's pipeline.
[isg4@farnam apps]$ git clone https://bitbucket.org/yanglab/phylogenomic_dataset_construction.gitThe pipeline has some prequisites that we need to install.
-
Rcorrector
Version: 1.0.3.1
[isg4@farnam apps]$ git clone https://github.com/mourisl/Rcorrector.git
[isg4@farnam apps]$ cd Rcorrector
[isg4@farnam Rcorrector]$ make-
TrimmomaticVersion: 0.036
[isg4@farnam ~]$ module load Trimmomatic-
TrinityVersion: 2.6.5
[isg4@farnam ~]$ module load Trinity/2.6.5-foss-2016b-
Bowtie2Version: 2.3.4.1
[isg4@farnam ~]$ module load Bowtie2/2.3.4.1-foss-2016b-
FastQCVersion: 0.11.5
[isg4@farnam ~]$ module load FastQC/0.11.5-Java-1.8.0_121-
TransrateVersion: 1.0.3
I had trouble using curl to get his to work so I downloaded a copy from the install page and then copied to to my cluster account using rsync.
-
TransDecoderVersion: 5.5.0
[isg4@farnam apps]$ git clone https://github.com/TransDecoder/TransDecoder.git-
BLASTVersion: 2.5.0+
[isg4@farnam ~]$ conda install -c bioconda blast-
SAMtools
Version: 1.9
[isg4@farnam ~]$ module load SAMtools/1.9-foss-2016b-
HISAT2
Version: 2.1.0
[isg4@farnam ~]$ module load HISAT2/2.1.0-foss-2016bThe Yang pipeline specifies that read files needs to be in the format taxonID_1.fq and taxonID_2.fq. I renamed my files ERR2040261_1.fq and ERR2040261_2.fq.
filter_fq.py command format:
python filter_fq.py [forward_reads] [reverse_reads] [order_name] [organelle] [cores] [output_dir]-
clean: removes intermediate files.
#!/bin/bash
#SBATCH --job-name=filter_fq
#SBATCH -c 20
#SBATCH --mem=36G
#SBATCH -t 60:00:00
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH -o filter_fq.%A-%a.out
#SBATCH -e filter_fq.%A-%a.err
conda activate genome
python ~/project/apps/phylogenomic_dataset_construction/scripts/filter_fq.py ../Data/1kp_transcriptome/ERR2040261_1.fq ../Data/1kp_transcriptome/ERR2040261_2.fq Caryophyllales both 20 /home/isg4/scratch60/P_amilis_genome/Transcriptomics clean > filter_fq.logThe above command will
- Correct random sequencing error with
Rcorrector - Remove read pairs that cannot be corrected
- Remove sequencing adapters and low quality sequences with
Trimmomatic - Filter organelle reads (cpDNA, mtDNA or both) with
Bowtie2. Files containing only organelle reads will be produced which can be use to assemble for example the plastomes with Fast-Plast - Run
FastQCto check read quality and detect over-represented reads - Remove over-represented sequences
Downstream analyses with Trinity require a coordinate-sorted .bam file, which just means that the reads are sorted in the file by their position on the genome rather than randomly. These coordinates will be used to bin reads based on their location. We'll use hisat2 to fix split out genome into multiple indices and then map reads to those indices.
hisat2-build command format:
hisat2-build [genome_fasta] [index_prefix][isg4@farnam HISAT-indices]$ hisat2-build ../Data/portulaca_26Nov2018_oK3Ko.fasta P_amilisThis results in 8 indices that look like P_amilis.1.ht2. Now we'll map reads from the transcriptome to these indices to generate a .sam file. The following batch script will creates a path variable (HISAT2_INDEXES) that hisat2 recognizes, map the reads using hisat2, and convert the output to the required .bam format using samtools.
hisat2 command format:
hisat2 -p [cores] --rg-id [read_group_ID] --rg [read_group_info] --summary-file [summary_output] --new-summary -x [index_prefix] -U [forward_reads,reverse_reads] | samtools view -bhuS - | samtools sort - -m [max_memory] -o [output_bam]-
--rg: This flag lets users to associate metadata with read groups. For example, the sequencing platform can be specified as Illumina withPU:ILLUMINA. See tags in thesamtoolsmanual for help -
--new-summary: Print output in a machine-friendly format. -
-x: Note that the index prefix is in the format of a path, for example/home/isg4/scratch60/P_amilis_genome/HISAT2/indices/P_amilisfor a prefix ofP_amilis
#!/bin/bash
#SBATCH --job-name=hisat2bam
#SBATCH -c 8
#SBATCH --mem=48G
#SBATCH -t 48:00:00
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH -o hisat2bam.%A-%a.out
#SBATCH -e hisat2bam.%A-%a.err
HISAT2_INDEXES=/home/isg4/scratch60/P_amilis_genome/HISAT2/indices
hisat2 -p 20 --rg-id "ERR2040261" --rg "LB:ERR2040261" --rg "PU:ILLUMINA" --rg "SM:ERR2040261_leaf_mature" --summary-file /home/isg4/scratch60/P_amilis_genome/hisat2 --rg-id "ERR2040261" --rg "LB:ERR2040261" --rg "PU:ILLUMINA" --rg "SM:ERR2040261_leaf_mature" --summary-file /home/isg4/scratch60/P_amilis_genome/HISAT2/ERR2040261-alnstats.txt --new-summary -x /home/isg4/scratch60/P_amilis_genome/HISAT2/indices/P_amilis -U /home/isg4/scratch60/P_amilis_genome/Transcriptomics/ERR2040261_1.overep_filtered.fq,/home/isg4/scratch60/P_amilis_genome/Transcriptomics/ERR2040261_2.overep_filtered.fq | samtools view -bhuS - | samtools sort - -m 2G -o ERR2040261.bam
samtools index ERR2040261.bamI renamed the ERR2040261.bam file ERR2040261.coordSorted.bam to explicitly denote that it is sorted. Taking a look at the hisat2 summary stats:
| Total reads | 17281410 |
|---|---|
| Aligned 0 time | 752513 (4.35%) |
| Aligned 1 time | 16261258 (94.10%) |
| Aligned >1 times | 267639 (1.55%) |
| Overall alignment rate | 95.65% |
we can see that the majority of our reads mapped to unique regions of the genome.
Because I have an assembled genome I'll deviate from Ya Yang's pipeline here. We'll use the P. amilis genome scaffolds for a guided Trinity assembly. The guided version of Trinity takes a hybrid approach to assembly: reads are first partitioned into bins based on their physical location on the genome, and then those bins undergo local de novo assemblies. This means that the final transcripts are not stitched together genomic sequences, as in some methods like cufflinks. This is beneficial for two main reasons. First, my ESTs and reference genome sequences come from two different organisms, so they may not exactly match in sequence content. Second, this method has potential to retain more paralogs and isoforms, which can be collapsed into consensus sequences when constructing ESTs from genomic sequences. While some of these isoforms may artifacts of sequencing error, misassembly, etc., recent long-read sequencing of the maize transcriptome has shown that many of these are true ESTs (Wang et al. 2016), and therefore should be retained for later estimation of gene models.
Trinity command format:
Trinity --genome_guided_bam [coordinate_sorted_bam] --genome_guided_max_intron [max_intron_size] --max_memory [max_memory] --CPU [cores]#!/bin/bash
#SBATCH --job-name=trinity_guided
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=120G
#SBATCH -t 72:00:00
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH -o trinity_guided.%A-%a.out
#SBATCH -e trinity_guided.%A-%a.err
conda activate genome
Trinity --genome_guided_bam /home/isg4/scratch60/P_amilis_genome/HISAT2/ERR2040261.coordSorted.bam --genome_guided_max_intron 10000 --max_memory 110G --CPU 20See /P_amilis_genome/Transcriptomics/Trinity-guided/trinity_out_dir/Trinity-GG.fasta for reconstructed transcripts. In Yang's pipeline the output transcripts from trinity are put through three filtering steps:
-
transrate, which filters reads based on internal quality metrics -
run_chimera_detection.py, which Yang and Smith (2013) describe -
corset, extracts a representative sequence from clusters of transcripts.
But I'm going to skip these for now because we guided our assembly with the genome, which should increase the quality of our assembly and decrease the chance of retaining chimeric sequences.
ISSUE Currently my MAKER2 runs are using the A. thaliana and B. vulgaris proteomes as guides for building gene models. This probably isn't a good idea because these are fairly distantly related taxa. I'll need to work on getting a proteome for Portulaca
Fortunately for me, Dr. Yang also works on taxa in the Caryophyllales and uses the A. thaliana and B. vulgaris proteomes to build a blast database. I'll use the same file, AthaBeta.fa, which is packaged with the repository.
[isg4@farnam Transcriptomics]$ mkdir Transdecoder
[isg4@farnam Transcriptomics]$ cd Transdecoder
[isg4@farnam Transdecoder]$ cat ~/project/apps/phylogenomic_dataset_construction/databases/Atha.fa ~/project/apps/phylogenomic_dataset_construction/databases/Beta.fa > AthaBeta.fa
[isg4@farnam Transdecoder]$ makeblastdb -in ./AthaBeta.fa -parse_seqids -dbtype prot -out AthaBeta.dbNow we can run the transdecoder wrapper.
Editting transdecoder_wrapper.py
I hit a number of issues using this wrapper script, which I think stem from python 2 vs 3 incompatibilities.
- I added parens around print statements using a quick search-and-replace Vim statement:
%s/print \(.*\)/print(\1)/c - I commented out line 133 because there was an
AssertionError. This error was being caused by asking if theblastdatabase exists, but this database is referenced as a path to the database and the database prefix (e.g./Transdecoder/mydatabase, which tells the program to look in/Transdecoder/forblastdatabase files with the prefixmydatabase) and therefore is not a real file. - Files aren't being opened and closed correctly in lines 170-192. I substituted a
file.open()andfile.close()method for awith file.open()method. Here are my rewritten lines 170-192:
##ISGinfile = gzip.open(cds,"rb")
##ISGfor line in lines:
##ISG if line[0] == ">":
##ISG cdsids.append((line.split(" ")[0][1:]))
##ISGfor line in infile:
##ISG if line[0] == ">": cdsids.append((line.split(" ")[0])[1:])
##ISGinfile.close()
with gzip.open(cds, "rb") as f:##ISG
cds_lines = f.readlines()##ISG
for line in cds_lines:##ISG
if line[0]==">":##ISG
cdsids.append(line.split()[0][1:])##ISG
assert len(cdsids) > 0, cds+"is not in correct fasta format"
##ISGinfile = gzip.open(blastpout,"rb")
with gzip.open(blastpout, "rb") as b:##ISG
query_list = [] # list of all query ids
hitcov_dict = {} # key is hit seqid, value is the max hit coverage
for line in b.readlines():
spls = line.strip().split("\t")
query,hit,piden = spls[0],spls[1],float(spls[2])
if query not in cdsids:
continue # only the ones get translated matters
# Turns out that the cds file contains all seqs that has blastp hits
# Filtering for ids in cds file doesn't really matter
query_list.append(query)
if piden < min_pident or not hit.startswith(hitID):
continue # only look at highly similar hits to the closest proteome
hitcov = int(spls[9])-int(spls[8])
if hit not in hitcov_dict:
hitcov_dict[hit] = hitcov
else: hitcov_dict[hit] = max(hitcov_dict[hit],hitcov)
##ISGinfile.close()Command format
python transdecoder_wrapper.py [transcripts] [cores] [strandedness] [output_directory]#!/bin/bash
#SBATCH --job-name=transdecoder_wrapper
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=36G
#SBATCH -t 72:00:00
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH -o transdecoder_wrapper.%A-%a.out
#SBATCH -e transdecoder_wrapper.%A-%a.err
conda activate genome
python /home/isg4/project/apps/phylogenomic_dataset_construction/scripts/transdecoder_wrapper.py ../Trinity-guided/trinity_out_dir/Trinity-GG.fasta 20 non-stranded .The output files we're mostly concerned with are Trinity-GG.pep.fa and Trinity-GG.cds.fa. Dr. Yang goes on to do some phylogenetically-aware clustering and paralog pruning but since we mapped everything to the reference genome we shouldn't have to worry as much about these problems. At this point I'll copy the transcriptomics directory with all of these analyses to project space for long term storage and use.
- Not critical: compare this transcriptome assembly with that in the 1KP database for 1kp_code QMUN
- Reassess the use of
TransDecoderhere. Right now I'm not using the translated coding sequences (Trinity-GG.pep.fa) in MAKER2, rather I'm using the Arabidopsis and Beta proteomes from earlier.
I'll be following Daren Card's MAKER annotation pipeline and the MAKER Tutorial for WGS Assembly and Annotation.
The first step in de novo genome annotation is building a species specific repeat library for masking purposes. This is done in the MAKER pipeline with RepeatModeler. RepeatModeler collects sequences based on a threshold copy number and then classifies them based on similarity to known transposable elements (TEs). Note that low copy TEs will not be included by RepeatModeler. First we need to install RepeatMasker and RepeatModeler.
Installing RepeatMasker | GitHub
First, download RepeatMasker from Github.
[isg4@farnam ~]$ cd $APPS
[isg4@farnam ~]$ git clone https://github.com/rmhubley/RepeatMasker.gitVersion: 1.332
Then install/download some dependencies.
`conda install -c bioconda trf`Version: 4.09
-
RepBase
Download:
RepBaseRepeatMaskerEdition-20181026.tar.gz, or whatever the latest version is, from the Genetic Information Research Institute and move it to theRepeatMaskerdirectory.
[isg4@farnam ~]$ cd $APPS/RepeatMasker
[isg4@farnam RepeatMasker]$ gunzip RepBaseRepeatMaskerEdition-20181026.tar.gz
[isg4@farnam RepeatMasker]$ tar -xvf RepBaseRepeatMaskerEdition-20181026.tar.gz-
Dfam Download
Dfam.hmm.gzfrom here and save it toRepeatMasker/LibrariesVersion: 23-Sep-2015 05:23Notes: Unfortunately, the Dfam 2.0 database only contains data from human, mouse, zebrafish, fruit fly, and nematode. Dfam 3.0 was released in March 2019 and contains "partial libraries for four additional species", although I can't find out what taxa these are. Since I had runRepeatModelerandRepeatMaskerbefore Dfam 3.0 was released, I'll need to go back and update the repeating masking.
[isg4@farnam ~]$ conda install -c bioconda trfVersion: 4.09
Lastly, we need to configure RepeatMasker once all the dependencies are installed.
[isg4@farnam apps]$ cd RepeatMasker
[isg4@farnam apps]$ perl ./configureRepeatModeler | GitHub
RepeatModeler also has some dependencies we'll need to install but first we'll clone it from GitHub.
[isg4@farnam ~]$ cd $APPS
[isg4@farnam ~]$ git clone https://github.com/rmhubley/RepeatModeler.gitVersion: 1.0.11
Now to the dependencies.
-
RECON(available at RepeatModeler)
[isg4@farnam apps]$ curl -O http://www.repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz
[isg4@farnam apps]$ gunzip RECON-1.08.tar.gz
[isg4@farnam apps]$ cd RECON-1.08
[isg4@farnam RECON-1.08]$ make
[isg4@farnam RECON-1.08]$ make installVersion: 1.08
-
RepeatScout(available at RepeatModeler)
[isg4@farnam apps]$ curl -O http://www.repeatmasker.org/RepeatScout-1.0.5.tar.gz
[isg4@farnam apps]$ gunzip RepeatScout-1.0.5.tar.gz
[isg4@farnam apps]$ cd RepeatScout-1
[isg4@farnam RepeatScout-1]$ make
`Version: 1.0.5`Now we're ready to build a repeat database using the P. amilis genome.
[isg4@farnam RepeatModeler]$ BuildDatabase -name Portulaca_amilis -engine ncbi ../Data/portulaca_26Nov2018_oK3Ko.fasta
>>> Building database Portulaca_amilis:
>>> Adding ../Data/portulaca_26Nov2018_oK3Ko.fasta to database
>>> Number of sequences (bp) added to database: 4053 ( 403885173 bp )Next, as a batch job, run RepeatModeler using the Portulaca_amilis database.
#!/bin/bash
#SBATCH --job-name=repeat_modeler
#SBATCH --time=144:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=36G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=repeatmodeler-%A-%a.out
#SBATCH --error=repeatmodeler-%A-%a.err
#SBATCH --ntasks-per-node=1
/home/isg4/project/apps/RepeatModeler/RepeatModeler -pa 36 -engine ncbi -database Portulaca_amilisThe main outputs from this run of RepeatModeler are in the directory /home/isg4/scratch60/P_amilis_genome/RepeatModeler/RM_28676.MonMar41925222019, which contains the following files of interest:
-
consensi.fa: Consensus sequences for each family identified -
families.stk: Seed alignments for each family identified
- Make sure to add data to above libraries when done
After identification and classification of repeats with RepeatModeler, 1464/2130 (68.7%) repeats were still unclassified. TEclass classifies unknown transpsosable elements into DNA transposons, LTRs, LINEs, and SINEs, and identified ORFs.
The consensi.fa.classified file from RepeatModeler was input into the TEclass online tool. Here is a link to the output.
Repeat statistics:
| Total | 2118 |
|---|---|
| DNA transposons | 686 |
| LTRs | 610 |
| LINEs | 292 |
| SINEs | 67 |
| Unclear | 463 |
The output file P_amilis.consensi.fa.classified-lib.txt was reformatted for consistency with RepeatModeler output using regex:
- find:
^>(\w*.*)#(Unknown)(.*)\|TEclass result: (\w+)\|.*; replace:>$1#$4 $3 - find:
unclear; replace:Unknown
It was saved as P_amilis.consensi.classified.fa.
For an indepth overview of RepeatMasker, see this chapter from Sebastian Tempel. As in the Boa annotation, we'll do multiple rounds of repeat masking, with the output of a given run then fed into the next round of masking as input.
For this first run I'll use the P. amilis custom library and output to directory P_amilis_mask.
Command format:
RepeatMasker -a -gff -x -pa [N cores] -e [search engine] -lib [repeat library] [genome fasta]
-
-a: Writes alignments in .align output file -
-gff: Creates an additional Gene Feature Finding format output -
-x: Returns repetitive regions masked with Xs rather than Ns
RepeatMasker -a -gff -x -pa 20 -e ncbi -lib /home/isg4/scratch60/P_amilis_genome/TEclass/P_amilis.consensi.classified.fa ../Data/portulaca_26Nov2018_oK3Ko.fasta
#!/bin/bash
#SBATCH --job-name=repeat_masker
#SBATCH --time=72:00:00
#SBATCH -c 16
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=repeatmasker-%A-%a.out
#SBATCH --error=repeatmasker-%A-%a.err
#SBATCH --ntasks-per-node=1
/home/isg4/project/apps/RepeatMasker/RepeatMasker -a -gff -x -pa 20 -e ncbi -lib /home/isg4/scratch60/P_amilis_genome/TEclass/P_amilis.consensi.classified.fa ../Data/portulaca_26Nov2018_oK3Ko.fastaHere is the resulting .tbl file statistics:
==================================================
file name: portulaca_26Nov2018_oK3Ko.fasta
sequences: 4053
total length: 403885173 bp (400444544 bp excl N/X-runs)
GC level: 37.05 %
bases masked: 186495552 bp ( 46.18 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
SINEs: 32632 4910134 bp 1.22 %
ALUs 0 0 bp 0.00 %
MIRs 0 0 bp 0.00 %
LINEs: 90978 24862301 bp 6.16 %
LINE1 23598 10541551 bp 2.61 %
LINE2 154 44815 bp 0.01 %
L3/CR1 0 0 bp 0.00 %
LTR elements: 175015 75606386 bp 18.72 %
ERVL 0 0 bp 0.00 %
ERVL-MaLRs 0 0 bp 0.00 %
ERV_classI 0 0 bp 0.00 %
ERV_classII 0 0 bp 0.00 %
DNA elements: 259282 53048834 bp 13.13 %
hAT-Charlie 359 90896 bp 0.02 %
TcMar-Tigger 0 0 bp 0.00 %
Unclassified: 104129 21755512 bp 5.39 %
Total interspersed repeats:180183167 bp 44.61 %
Small RNA: 3162 598313 bp 0.15 %
Satellites: 552 406101 bp 0.10 %
Simple repeats: 132909 6308671 bp 1.56 %
Low complexity: 17954 863838 bp 0.21 %
==================================================To keep repeat libraries clearly labeled I've renamed most output files by removing .fasta. or .fasta.align. to Pamilis or Pamilis_.
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align.landscape.Div.Rclass.tab portulaca_26Nov2018_oK3Ko.Pamilis_landscape.Div.Rclass.tab
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align.landscape.Div.Rfam.tab portulaca_26Nov2018_oK3Ko.Pamilis_landscape.Div.Rfam.tab
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align.landscape.Div.Rname.tab portulaca_26Nov2018_oK3Ko.Pamilis_landscape.Div.Rname.tab
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align.parseRM.all-repeats.tab portulaca_26Nov2018_oK3Ko.Pamilis_parseRM.all-repeats.tab
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align.parseRM.summary.tab portulaca_26Nov2018_oK3Ko.Pamilis_parseRM.summary.tab
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked portulaca_26Nov2018_oK3Ko.Pamilis_masked.fasta
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align portulaca_26Nov2018_oK3Ko.Pamilis_masked.align
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.cat.gz portulaca_26Nov2018_oK3Ko.Pamilis_masked.cat.gz
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.masked portulaca_26Nov2018_oK3Ko.Pamilis_masked.masked
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.out portulaca_26Nov2018_oK3Ko.Pamilis_masked.out
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.out.gff portulaca_26Nov2018_oK3Ko.Pamilis_masked.out.gff
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.tbl portulaca_26Nov2018_oK3Ko.Pamilis_masked.tbl
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.out portulaca_26Nov2018_oK3Ko.Pamilis.out
[isg4@farnam P_amilis_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.align portulaca_26Nov2018_oK3Ko.Pamilis.alignWe'll feed the masked output from custom library back into repeatmasker for the next round of masking instead of doing two separate masks on the same input scaffolds. For this run I'll use a Viridiplantae library assembled from plants on RepBase and output to directory Pamilis_Viridiplantae_mask. The repeat libraries for Viridiplantae come from grasses and Arabidopsis thaliana, so I don't expect too much to be added to our ab initio library. Note that -lib has been substituted for -species to tell RepeatMasker to search RepBase for all Viridiplantae repeats instead of a custom library.
#!/bin/bash
#SBATCH --job-name=repeat_masker
#SBATCH --time=72:00:00
#SBATCH -c 16
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=repeatmasker-%A-%a.out
#SBATCH --error=repeatmasker-%A-%a.err
#SBATCH --ntasks-per-node=1
/home/isg4/project/apps/RepeatMasker/RepeatMasker -a -gff -x -pa 20 -e ncbi -species viridiplantae /home/isg4/scratch60/P_amilis_genome/RepeatMasker/P_amilis_mask/portulaca_26Nov2018_oK3Ko.Pamilis_masked.fastaThe results show that almost all of the repeats masked in the Viridiplantae only mask were already detected in the de novo library.
==================================================
file name: portulaca_26Nov2018_oK3Ko.fasta.masked
sequences: 4053
total length: 403885173 bp (214147364 bp excl N/X-runs)
GC level: 37.24 %
bases masked: 2402689 bp ( 0.59 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
Retroelements 7149 1413325 bp 0.35 %
SINEs: 42 3347 bp 0.00 %
Penelope 0 0 bp 0.00 %
LINEs: 1777 258286 bp 0.06 %
CRE/SLACS 2 77 bp 0.00 %
L2/CR1/Rex 0 0 bp 0.00 %
R1/LOA/Jockey 0 0 bp 0.00 %
R2/R4/NeSL 0 0 bp 0.00 %
RTE/Bov-B 84 12027 bp 0.00 %
L1/CIN4 1689 246220 bp 0.06 %
LTR elements: 5330 1151692 bp 0.29 %
BEL/Pao 0 0 bp 0.00 %
Ty1/Copia 2902 740890 bp 0.18 %
Gypsy/DIRS1 2202 350576 bp 0.09 %
Retroviral 0 0 bp 0.00 %
DNA transposons 3308 707766 bp 0.18 %
hobo-Activator 741 232569 bp 0.06 %
Tc1-IS630-Pogo 66 23441 bp 0.01 %
En-Spm 0 0 bp 0.00 %
MuDR-IS905 0 0 bp 0.00 %
PiggyBac 0 0 bp 0.00 %
Tourist/Harbinger 368 68250 bp 0.02 %
Other (Mirage, 3 147 bp 0.00 %
P-element, Transib)
Rolling-circles 0 0 bp 0.00 %
Unclassified: 711 132855 bp 0.03 %
Total interspersed repeats: 2253946 bp 0.56 %
Small RNA: 706 93632 bp 0.02 %
Satellites: 127 14073 bp 0.00 %
Simple repeats: 558 35176 bp 0.01 %
Low complexity: 119 6682 bp 0.00 %
==================================================The outputs were again renamed to reflect both masks.
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.align
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align.landscape.Div.Rclass.tab portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_landscape.Div.Rclass.tab
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align.landscape.Div.Rfam.tab portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_landscape.Div.Rfam.tab
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align.landscape.Div.Rname.tab portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_landscape.Div.Rname.tab
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align.parseRM.all-repeats.tab portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_parseRM.all-repeats.tab
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.align.parseRM.summary.tab portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_parseRM.summary.tab
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.cat.gz portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.cat.gz
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.masked portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.fasta
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.out portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.out
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.out.gff portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.out.gff
[isg4@farnam Pamilis_Viridiplantae_mask]$ mv portulaca_26Nov2018_oK3Ko.fasta.masked.tbl portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.tblNow we'll combine the results from the two repeat libraries and polish them with ProcessRepeats. From Sebastian Tempel
ProcessRepeats follows the RepeatMasker application. It organizes and processes the results by assembling and sorting the fragmented repeats. Many repeat subfamilies are similar to each other and could be close in the sequence, thus the main difficulty of ProcessRepeats algorithm is to merge fragmented repeats and to assign appropriate subfamily names to the search engine output.
This process involves copying the final repeatmasker output (Pamilis_Viridiplantae) .fasta and .out files to a new directory that I called Full_mask, unzipping the .cat mask files, and combining the unzipped .cat files, which ProcessRepeats takes as an input.
[isg4@farnam Full_mask]$ cp ../Pamilis_Viridiplantae_mask/portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.fasta ./portulaca_26Nov2018_oK3Ko.full_mask.fasta
[isg4@farnam Full_mask]$ cp ../Pamilis_Viridiplantae_mask/portulaca_26Nov2018_oK3Ko.Pamilis_Viridiplantae_masked.out ./portulaca_26Nov2018_oK3Ko.full_mask.out
[isg4@farnam Full_mask]$ gunzip ../Pamilis_Viridiplantae_mask/*.cat.gz ../P_amilis_mask/*.cat.gz
[isg4@farnam Full_mask]$ cat ../P_amilis_mask/*.cat ../Pamilis_Viridiplantae_mask/*.cat > portulaca_26Nov2018_oK3Ko.full_mask.catHere is the final ProcessRepeats command which
#!/bin/bash
#SBATCH --job-name=process_repeats
#SBATCH --time=48:00:00
#SBATCH -c 16
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=processrepeats-%A-%a.out
#SBATCH --error=processrepeats-%A-%a.err
#SBATCH --ntasks-per-node=1
conda activate genome
ProcessRepeats -species viridiplantae ./portulaca_26Nov2018_oK3Ko.full_mask.catThe final, combined, processed repeat library stats are below.
==================================================
file name: portulaca_26Nov2018_oK3Ko.full_mask
sequences: 4053
total length: 403885173 bp (218727942 bp excl N/X-runs)
GC level: Unknown %
bases masked: 190195790 bp ( 47.09 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
Retroelements 305777 106783617 bp 26.44 %
SINEs: 32675 4913584 bp 1.22 %
Penelope 610 67079 bp 0.02 %
LINEs: 92758 25120119 bp 6.22 %
CRE/SLACS 2 77 bp 0.00 %
L2/CR1/Rex 154 44815 bp 0.01 %
R1/LOA/Jockey 0 0 bp 0.00 %
R2/R4/NeSL 1361 199934 bp 0.05 %
RTE/Bov-B 22907 4438791 bp 1.10 %
L1/CIN4 25289 10787434 bp 2.67 %
LTR elements: 180344 76749914 bp 19.00 %
BEL/Pao 887 227343 bp 0.06 %
Ty1/Copia 34581 15478925 bp 3.83 %
Gypsy/DIRS1 46749 33227704 bp 8.23 %
Retroviral 0 0 bp 0.00 %
DNA transposons 262562 53752744 bp 13.31 %
hobo-Activator 37465 8597777 bp 2.13 %
Tc1-IS630-Pogo 35361 6202845 bp 1.54 %
En-Spm 0 0 bp 0.00 %
MuDR-IS905 0 0 bp 0.00 %
PiggyBac 0 0 bp 0.00 %
Tourist/Harbinger 10210 2510657 bp 0.62 %
Other (Mirage, 3 147 bp 0.00 %
P-element, Transib)
Rolling-circles 0 0 bp 0.00 %
Unclassified: 104842 21888128 bp 5.42 %
Total interspersed repeats: 182424489 bp 45.17 %
Small RNA: 3871 691903 bp 0.17 %
Satellites: 679 419983 bp 0.10 %
Simple repeats: 133456 6343506 bp 1.57 %
Low complexity: 18042 870498 bp 0.22 %
==================================================These repeats now need to be converted to .gff3 format and simple repeats need to be removed. The reason that simple repeats need to be removed is that simple repeats can be part of valid gene models, therefore we do not want to hard mask them before we train our gene models later. Complex repeats, however, usually derive from transposable elements that can be inserted into intronic regions. These confuse gene prediction software downstream and lead to many missed or mis-annotated gene models.
# convert to .gff3
[isg4@farnam Full_mask]$ ~/project/apps/RepeatMasker/util/rmOutToGFF3.pl portulaca_26Nov2018_oK3Ko.full_mask.out > portulaca_26Nov2018_oK3Ko.full_mask.out.gff3
# isolate complex repeats
[isg4@farnam Full_mask]$ grep -v -e "Satellite" -e ")n" -e "-rich" portulaca_26Nov2018_oK3Ko.full_mask.out.gff3 > portulaca_26Nov2018_oK3Ko.full_mask.complex.gff3
# reformat complex reads for MAKER
[isg4@farnam Full_mask]$ cat portulaca_26Nov2018_oK3Ko.full_mask.complex.gff3 | perl -ane '$id; if(!/^\#/){@F = split(/\t/, $_); chomp $F[-1];$id++; $F[-1] .= "\;ID=$id"; $_ = join("\t", @F)."\n"} print $_' > portulaca_26Nov2018_oK3Ko.full_mask.complex.reformat.gff3Finally, we can reformat the ouputs so that we can generate a repeat landscape. I had some trouble directly using the accessory scripts included with repeatmasker (calcDivergenceFromAlign.pl and createRepeatLandscape.pl) so I'm using this Parsing-RepeatMasker-Outputs, which can be cloned and then we can use the scripts inside it without any additional installation steps.
[isg4@farnam ~]$ cd $APPS
[isg4@farnam apps]$ git clone https://github.com/4ureliek/Parsing-RepeatMasker-Outputs.gitFirst we'll parse the .align file to generate relavent stats (using the --parse flag) and then build the files needed for plotting (using the --land flag). My first run through this with 120G of memory quit because it exceeded the memory limit. I'm going to put it on the bigmem node instead.
#!/bin/bash
#SBATCH --job-name=parseRM
#SBATCH --time=72:00:00
#SBATCH -p bigmem
#SBATCH -N 1
#SBATCH -c 16
#SBATCH --mem 200G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=parseRM-%A-%a.out
#SBATCH --error=parseRM-%A-%a.err
perl ~/project/apps/Parsing-RepeatMasker-Outputs/parseRM.pl -i portulaca_26Nov2018_oK3Ko.full_mask.align -g 403885173 --parse
perl ~/project/apps/Parsing-RepeatMasker-Outputs/parseRM.pl -i portulaca_26Nov2018_oK3Ko.full_mask.align --land 50,1 -v- Generate final repeat landscape plot
I tried installing MAKER2 multiple times using conda builds and from source. The developers do not suggest using the conda build because there have been many issues posted, but the creator of the build claims it works when installed using the bioconda/label/cf201901 channel versions of MAKER2 and its dependencies. In the end, I ended up using an installation of MAKER2 from the Yale HPC team (thanks Jason), although it wasn't configured for MPI use. Furthermore, I repeatedly had trouble with symbol lookup errors when running MAKER2 that stemmed from different packages used by MAKER2 trying to use different perls. I worked around this by doing unset PERL5LIB before I'd call MAKER2.
MAKER2 can utilize a number of different programs for identifying different genomic features. We'll install some of these later, like SNAP and Augustus for ab initio gene finding, but I'll install tRNA-scan now. There are others that you may want to use depending on the data available and the organisms you're annotating.
Installing tRNA-scan
First, tRNA-scan requires a package called infernal that contains cmscan. I installed that with conda before
[isg4@farnam ~]$ conda install -c bioconda infernal=1.1.2
[isg4@farnam ~]$ export APPS=/home/isg4/apps
[isg4@farnam ~]$ cd $APPS
[isg4@farnam apps]$ curl -O http://trna.ucsc.edu/software/trnascan-se-2.0.2.tar.gz
[isg4@farnam apps]$ gunzip trnascan-se-2.0.2.tar.gz
[isg4@farnam apps]$ cd tRNAscan-SE-2.0Then I had to edit the infernal_dir: {bin_dir} line of tRNAscan-SE.conf and tRNAscan-SE.conf.src to reflect the path to my infernal installation as in this post.
[isg4@farnam apps]$ export HOME=/home/isg4
[isg4@farnam apps]$ ./configure --prefix=$HOME
[isg4@farnam apps]$ make
[isg4@farnam apps]$ make installWe first need ot set up the maker control file. A template can be created with maker -CTL. This creates
-
maker_exe.ctl: contains the path information for the underlying executables -
maker_bopt.ctl: contains filtering statistics for BLAST and Exonerate -
maker_opt.ctl: contains all other information for MAKER, including the location of the input genome file
Here are the sources of evidence we'll give to MAKER,
-
portulaca_26Nov2018_oK3Ko.fasta: the scaffolds of the P. amilis genome -
Trinity-GG.fasta: a P. amilis transcriptome resulting from a reference-guided assembly inTrinity AthaBeta.fa: a concatenated.fastafile of the Arabidopsis thaliana and Beta vulgaris proteomes from Ya Yang's transcriptomic pipeline
Update In this analysis I provide proteins from two, relatively unrelated species to MAKER2. This is likely a bad idea. See my FUNannotate notes on extracting high quality gene models. TD;DR Accuracy drops off sharply as protein evidence comes from more distant relatives than sister species.
Below is the final maker_opt.ctl file, which I renamed maker_opts.R1.ctl,
#-----Genome (these are always required)
genome=/home/isg4/scratch60/P_amilis_genome/Data/portulaca_26Nov2018_oK3Ko.fasta #genome sequence (fasta file or fasta embeded in GFF3 file)
organism_type=eukaryotic #eukaryotic or prokaryotic. Default is eukaryotic
#-----Re-annotation Using MAKER Derived GFF3
maker_gff= #MAKER derived GFF3 file
est_pass=0 #use ESTs in maker_gff: 1 = yes, 0 = no
altest_pass=0 #use alternate organism ESTs in maker_gff: 1 = yes, 0 = no
protein_pass=0 #use protein alignments in maker_gff: 1 = yes, 0 = no
rm_pass=0 #use repeats in maker_gff: 1 = yes, 0 = no
model_pass=0 #use gene models in maker_gff: 1 = yes, 0 = no
pred_pass=0 #use ab-initio predictions in maker_gff: 1 = yes, 0 = no
other_pass=0 #passthrough anyything else in maker_gff: 1 = yes, 0 = no
#-----EST Evidence (for best results provide a file for at least one)
est=/home/isg4/scratch60/P_amilis_genome/Transcriptomics/Trinity-guided/trinity_out_dir/Trinity-GG.fasta #set of ESTs or assembled mRNA-seq in fasta format
altest=/home/isg4/scratch60/P_amilis_genome/Transcriptomics/Transdecoder/Trinity-GG.cds.fa #EST/cDNA sequence file in fasta format from an alternate organism
est_gff= #aligned ESTs or mRNA-seq from an external GFF3 file
altest_gff= #aligned ESTs from a closly relate species in GFF3 format
#-----Protein Homology Evidence (for best results provide a file for at least one)
protein=/home/isg4/scratch60/P_amilis_genome/Transcriptomics/Transdecoder/AthaBeta.fa #protein sequence file in fasta format (i.e. from mutiple oransisms)
protein_gff= #aligned protein homology evidence from an external GFF3 file
#-----Repeat Masking (leave values blank to skip repeat masking)
model_org=viridiplantae #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff=/home/isg4/scratch60/P_amilis_genome/RepeatMasker/Full_mask/portulaca_26Nov2018_oK3Ko.full_mask.complex.reformat.gff3 #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)
#-----Gene Prediction
snaphmm= #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species= #Augustus gene prediction species model
fgenesh_par_file= #FGENESH parameter file
pred_gff= #ab-initio predictions from an external GFF3 file
model_gff= #annotated gene models from an external GFF3 file (annotation pass-through)
est2genome=1 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=1 #infer predictions from protein homology, 1 = yes, 0 = no
trna=0 #find tRNAs with tRNAscan, 1 = yes, 0 = no
snoscan_rrna= #rRNA file to have Snoscan find snoRNAs
unmask=0 #also run ab-initio prediction programs on unmasked sequence, 1 = yes, 0 = no
#-----Other Annotation Feature Types (features MAKER doesn't recognize)
other_gff= #extra features to pass-through to final MAKER generated GFF3 file
#-----External Application Behavior Options
alt_peptide=C #amino acid used to replace non-standard amino acids in BLAST databases
cpus=1 #max number of cpus to use in BLAST and RepeatMasker (not for MPI, leave 1 when using MPI)
#-----MAKER Behavior Options
max_dna_len=100000 #length for dividing up contigs into chunks (increases/decreases memory usage)
min_contig=1 #skip genome contigs below this length (under 10kb are often useless)
pred_flank=200 #flank for extending evidence clusters sent to gene predictors
pred_stats=1 #report AED and QI statistics for all predictions as well as models
AED_threshold=1 #Maximum Annotation Edit Distance allowed (bound by 0 and 1)
min_protein=0 #require at least this many amino acids in predicted proteins
alt_splice=1 #Take extra steps to try and find alternative splicing, 1 = yes, 0 = no
always_complete=0 #extra steps to force start and stop codons, 1 = yes, 0 = no
map_forward=0 #map names and attributes forward from old GFF3 genes, 1 = yes, 0 = no
keep_preds=0 #Concordance threshold to add unsupported gene prediction (bound by 0 and 1)
split_hit=10000 #length for the splitting of hits (expected max intron size for evidence alignments)
single_exon=0 #consider single exon EST evidence when generating annotations, 1 = yes, 0 = no
single_length=250 #min length required for single exon ESTs if 'single_exon is enabled'
correct_est_fusion=0 #limits use of ESTs in annotation to avoid fusion genes
tries=2 #number of times to try a contig if there is a failure for some reason
clean_try=0 #remove all data from previous run before retrying, 1 = yes, 0 = no
clean_up=0 #removes theVoid directory with individual analysis files, 1 = yes, 0 = no
TMP= #specify a directory other than the system default temporary directory for temporary filesand the maker_exe.ctl file, which I renamed maker_exe.ycrc.ctl.
#-----Location of Executables Used by MAKER/EVALUATOR
makeblastdb=/gpfs/ysm/project/isg4/conda_envs/genome/bin/makeblastdb #location of NCBI+ makeblastdb executable
blastn=/gpfs/ysm/project/isg4/conda_envs/genome/bin/blastn #location of NCBI+ blastn executable
blastx=/gpfs/ysm/project/isg4/conda_envs/genome/bin/blastx #location of NCBI+ blastx executable
tblastx=/gpfs/ysm/project/isg4/conda_envs/genome/bin/tblastx #location of NCBI+ tblastx executable
formatdb= #location of NCBI formatdb executable
blastall= #location of NCBI blastall executable
xdformat= #location of WUBLAST xdformat executable
blasta= #location of WUBLAST blasta executable
RepeatMasker=/home/isg4/project/apps/RepeatMasker/RepeatMasker #location of RepeatMasker executable
exonerate=/ysm-gpfs/apps/software/Exonerate/2.4.0-foss-2016b/bin/exonerate #location of exonerate executable
#-----Ab-initio Gene Prediction Algorithms
snap=/home/isg4/project/apps/SNAP/snap #location of snap executable
gmhmme3= #location of eukaryotic genemark executable
gmhmmp= #location of prokaryotic genemark executable
augustus=/ysm-gpfs/apps/software/Augustus/3.2.3-foss-2016b/bin/augustus #location of augustus executable
fgenesh= #location of fgenesh executable
tRNAscan-SE=/home/isg4/project/apps/tRNAscan-SE-2.0/tRNAscan-SE #location of trnascan executable
snoscan=/gpfs/ysm/project/isg4/conda_envs/genome/bin/snoscan #location of snoscan executable
#-----Other Algorithms
probuild= #location of probuild executable (required for genemark)After testing MAKER2 in an interactive session, I ran MAKER2 as a batch file.
#!/bin/bash
#SBATCH --job-name=MAKER
#SBATCH --time=120:00:00
#SBATCH --nodes=1
#SBATCH --cpus-per-task=1
#SBATCH --ntasks=1
#SBATCH --partition=general
#SBATCH --mem=60G
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=MAKER-%A-%a.out
#SBATCH --error=MAKER-%A-%a.err
/ysm-gpfs/project/jsi3/maker2/maker/bin/maker -base P_amilis_R1 maker_opts.R1.ctl maker_bopts.ctl maker_exe.ycrc.ctlMAKER2 took about 5 days too finish, see below.
| Job metric | Usage |
|---|---|
| Cores | 1 |
| CPU Utilized | 4-05:48:06 |
| CPU Efficiency | 91.44% of 4-15:20:10 core-walltime |
| Job Wall-clock time | 4-15:20:10 |
| Memory Utilized | 2.24 GB |
| Memory Efficiency | 3.74% of 60.00 GB |
So it took up most of the time I allotted but didn't strain the memory limits I set. This is in line with previous results that show MAKER2's speed increasing proportionally to the number of cores given. Sorry for wasting space Farnam users...
The MAKER2 outputs consists of the following for each scaffold:
-
run.log: a log of themakerrun for this scaffold -
<Scaffold>.maker.proteins.fasta: protein sequences for MAKER annotations -
<Scaffold>.maker.transcripts.fasta: transcript sequences for MAKER annotations -
<Scaffold>.gff: scaffold annotations in.gffformat -
theVoid.<Scaffold>/: a directory containing raw output files from all the programs MAKER ran
Not all scaffolds will have annotations, so to check out some that do, navigate to your <project>.maker.output directory and do find . -type f -name *.proteins.fasta to list some places where protein files exist.
Now we need to combine all of the MAKER annotations together using the accessory scripts gff3_merge and fasta_merge. Lastly, we'll calculate the number of predicted gene models and the mean length of each.
[isg4@farnam P_amilis_R1.maker.output]$ gff3_merge -s -d P_amilis_R1_master_datastore_index.log > P_amilis_R1.all.maker.gff
[isg4@farnam P_amilis_R1.maker.output]$ fasta_merge -d P_amilis_R1_master_datastore_index.log
# This last one creates a .gff without the sequences
[isg4@farnam P_amilis_R1.maker.output]$ gff3_merge -n -s -d P_amilis_R1_master_datastore_index.log > P_amilis_R1.all.maker.noseq.gff
[isg4@farnam P_amilis_R1.maker.output]$ cat P_amilis_R1.all.maker.noseq.gff | awk '{ if ($3 == "gene") print $0 }' | awk '{ sum += ($5 - $4) } END { print NR, sum / NR }'
23893 3661.38As noted in the Boa annotation walkthrough, the output of the initial MAKER run, which maps our empirical data to our scaffolds, is also an input for training gene prediction software like snap and augustus. The output from these predictors will then be used in subsequent MAKER analyses to refine gene models. For both trainings I'll work inside of the MAKER2 output directory, for me this is P_amilis_R1.maker.output.
Before going forward I should note that I hit a roadblock for a few days trying to get SNAP to work. The initial weird behavior was the genome.ann and genome.dna files were missing the sequence data, and only contained metadata (headers, etc.). If you end up with a similar problem, take a look at a possible bug with a possible solution.
Installing SNAP Version: 2006-07-28
[isg4@farnam P_amilis_R1.maker.output]$ cd ~/project/apps
[isg4@farnam apps]$ git clone https://github.com/KorfLab/SNAP.git
[isg4@farnam SNAP]$ cd SNAP
[isg4@farnam SNAP]$ make
[isg4@farnam SNAP]$ export ZOE=$PWD/ZoeI'll train SNAP using only those annotations with Annotation Edit Distance (AED) less than 0.25 (-x 0.25) and who's length is 50+ amino acids (-l 50). As noted in Yandell and Ence (2012)
An AED of 0 indicates that the annotation is in perfect agreement with its evidence, whereas an AED of 1 indicates a complete lack of evidence support for the annotation.
So we're being a bit more strict than the standard maker2zff settings (no minimum length and AED
[isg4@farnam SNAP]$ cd ~/scratch60/P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output
[isg4@farnam P_amilis_R1.maker.output]$ export MAKERBIN=/ysm-gpfs/project/jsi3/maker2/maker/bin/
[isg4@farnam P_amilis_R1.maker.output]$ mkdir -p snap/round1
[isg4@farnam P_amilis_R1.maker.output]$ cd snap/round1
[isg4@farnam round1]$ $MAKERBIN/maker2zff -x 0.25 -l 50 ../../P_amilis_R1.all.maker.fixed.gff
[isg4@farnam round1]$ ~/project/apps/SNAP/fathom genome.ann genome.dna -gene-stats > gene-stats.log 2>&1
[isg4@farnam round1]$ ~/project/apps/SNAP/fathom genome.ann genome.dna -validate > validate.log 2>&1
[isg4@farnam round1]$ ~/project/apps/SNAP/fathom genome.ann genome.dna -categorize 1000 > categorize.log 2>&1
[isg4@farnam round1]$ ~/project/apps/SNAP/fathom uni.ann uni.dna -export 1000 -plus > uni-plus.log 2>&1This will result in
-
alt.ann,alt.dna(genes with alternative splicing) -
err.ann,err.dna(genes that have errors) -
olp.ann,olp.dna(genes that overlap other genes) -
wrn.ann,wrn.dna(genes with warnings) -
uni.ann,uni.dna(single gene per sequence)
Now we'll assemble the training paramers.
[isg4@farnam round1]$ mkdir params
[isg4@farnam round1]$ cd params
[isg4@farnam params]$ ~/project/apps/SNAP/forge ../export.ann ../export.dna > ../forge.log 2>&1This will result in files that looks like Acceptor*, Coding*, Donor*, etc.
[isg4@farnam params]$ cd ..
[isg4@farnam round1]$ ~/project/apps/SNAP/hmm-assembler.pl genome params > P_amilis.length50.aed0.25.hmmIf you're having the a problem where genome.ann and genome.dna files were missing the sequence data, this section may help resolve that issue. If not, continue to the installing and training Augustus. This problem can occur if your gff from MAKER doesn't contain any information about genes. I checked this by using a few greps.
[isg4@farnam P_amilis_R1.maker.output]$ grep "CDS" P_amilis_R1.all.maker.gff | wc -l
163117
[isg4@farnam P_amilis_R1.maker.output]$ grep "mRNA" P_amilis_R1.all.maker.gff | wc -l
365503
[isg4@farnam P_amilis_R1.maker.output]$ grep "gene" P_amilis_R1.all.maker.gff | wc -l
389396
[isg4@farnam P_amilis_R1.maker.output]$ grep "exon" P_amilis_R1.all.maker.gff | wc -l
389396So that wasn't the issue. Empty files can also result from a lack of features that pass filtering thresholds, so I tried using the -n flag that tells maker2zff to grab everything, regardless of quality,
$ /ysm-gpfs/project/jsi3/maker2/maker/bin/maker2zff -n ../../P_amilis_R1.all.maker.gff
Use of uninitialized value in print at /ysm-gpfs/project/jsi3/maker2/maker/bin/maker2zff line 171, <GFF> line 10251786.but this did not fix the issue either. This is when I realized that my "empty" files are not "empty" as other users have described--mine actually contain .fasta headers with an empty line each, not a completely empty file. This suggested that there is a bug in maker2zff that was preventing sequences from being appended after headers, which, in light of the warning made sense because the "uninitialized value" in line 171 prints the header and sequence to the output file. So, after checking that everything in maker2zff looked okay, I thought that there must be something about how sequences are referenced that prevented them from being associated with headers (or IDs). The most common way to handle sequence IDs and sequences is with some dictionary-like object that uses the IDs as keys and sequences as values. And, in fact, when I took a look as my .gff from gff3_merge, I noticed it refers to scaffolds like Scaffold_27%3BHRSCAF%3D30, but my actual scaffolds are identified like this Scaffold_27;HRSCAF=30. However, after reading documentation on proper GFF3 format, this list of escaped characters suggests that the original file may have been valid.
In addition, the following characters have reserved meanings in column 9 and must be escaped when used in other contexts:
- ; semicolon (%3B)
- = equals (%3D)
- & ampersand (%26)
- , comma (%2C)
Because I wasn't sure if everything was internally valid (Scaffold_27%3BHRSCAF%3D30 in the ##gff-version 3 portion, but Scaffold_27;HRSCAF=30 in the ##FASTA portion of P_amilis_R1.all.maker.gff), so I converted all of the headers to the latter format in iPython. I then used gff3validator from genometools to check if the resulting file was properly formatted.
[isg4@farnam P_amilis_R1.maker.output]$ conda install -c bioconda genometools-genometools
[isg4@farnam P_amilis_R1.maker.output]$ gt gff3validator P_amilis_R1.all.maker.fixed.gff
gt gff3validator: error: illegal uppercase attribute "HRSCAF" on line 2 in file "P_amilis_R1.all.maker.fixed.gff" (uppercase attributes are reserved)Finally, after looking back at the GFF3 format guide, I thought that the ; in the sequence IDs may be screwing the delimitation of various fields.
Column 1: "seqid" The ID of the landmark used to establish the coordinate system for the current feature. IDs may contain any characters, but must escape any characters not in the set [a-zA-Z0-9.:^*$@!+_?-|]. In particular, IDs may not contain unescaped whitespace and must not begin with an unescaped ">".
I also checked the original .fasta headers from the genome scaffolds and they have they two unique identifiers: >Scaffold_865;HRSCAF=1003, so I dropped the second (search: ;HRSCAF=[\d]+, replace: nothing), which made all my sequence IDs look like >Scaffold_865--perfectly kosher IDs. So, finally, all of my SNAP commands appear to be in working order.
Installing Augustus
These next commands install and configure Augustus and BUSCO (this module contains both the BUSCO libraries and the script run_BUSCO.py), which may look a little different for me because I'm using a combination of Yale HPC- and self-installed software. This requires SAMtools, BamTools, BCFtools, and HTSlib, which we'll have to configure Augustus to find when we run our make command.
[isg4@farnam ~]$ module load BUSCO/3.0.2-foss-2016b-Python-2.7.13 SAMtools/1.9-foss-2016b BamTools/2.5.1-foss-2016b HTSlib/1.5-foss-2016b BCFtools/1.5-foss-2016b
[isg4@farnam ~]$ export BUSCO=/ysm-gpfs/apps/software/BUSCO/3.0.2-foss-2016b-Python-2.7.13
[isg4@farnam ~]$ cd $APPS
[isg4@farnam apps]$ git clone https://github.com/Gaius-Augustus/Augustus.git
[isg4@farnam apps]$ cd AugustusThe Augustus documentation tells you to run a make command here, but as this installation guide from Jia-Xing Yue points out, you'll often run into problems with how Augustus references bamtools. Here are all the edits I made to various Makefiles and the final installation.
-
Augustus/auxprogs/bam2hintsq/Makefile- Old:
INCLUDES = /usr/include/bamtools - New:
INCLUDES = $(BAMTOOLS)/include/bamtools - Old:
LIBS = -lbamtools -lz - New:
LIBS = $(BAMTOOLS)/lib64/libbamtools.a -lz
- Old:
-
Augustus/auxprogs/filterBam/src/Makefile- Old:
BAMTOOLS = /usr/include/bamtools - Comment out:
# BAMTOOLS = /usr/include/bamtools - Old:
INCLUDES = -I$(BAMTOOLS) -Iheaders -I./bamtools - New:
INCLUDES = -I$(BAMTOOLS)/include/bamtools -Iheaders -I./bamtools - Old:
LIBS = -lbamtools -lz - New:
LIBS = $(BAMTOOLS)/lib64/libbamtools.a -lz
- Old:
-
Augustus/auxprogs/bam2wig/Makefile- Old:
SAMTOOLS=$(TOOLDIR)/samtools - Comment out:
# SAMTOOLS=$(TOOLDIR)/samtools - Old:
HTSLIB=$(TOOLDIR)/htslib - Comment out:
# HTSLIB=$(TOOLDIR)/htslib - Old:
BCFTOOLS=$(TOOLDIR)/bcftools - Comment out:
# BCFTOOLS=$(TOOLDIR)/bcftools
- Old:
Now when we run make we'll export the paths we just configured.
[isg4@farnam Augustus]$ make BAMTOOLS=/ysm-gpfs/apps/software/BamTools/2.5.1-foss-2016b SAMTOOLS=/ysm-gpfs/apps/software/SAMtools/1.9-foss-2016b/bin/ BCFTOOLS=/ysm-gpfs/apps/software/BCFtools/1.5-foss-2016b/bin/ HTSLIB=/ysm-gpfs/apps/software/HTSlib/1.5-foss-2016b/bin/Even after all of this additional configuration I still had a warning that bam2wig could not be installed because sam.h couldn't be found. This is annoying, but it worked. When I load Augustus it resets my AUGUSTUS_CONFIG_PATH, so I have to set it after I load it.
[isg4@farnam Augustus]$ module load Augustus
[isg4@farnam Augustus]$ cd config
[isg4@farnam config]$ export AUGUSTUS_CONFIG_PATH=$PWD
[isg4@farnam config]$ cd ~/scratch60/P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output/Augustus/round1/Similar to Augustus, run_BUSCO.py requires modifiable configuration files (config.ini), but these do not require a particular directory hierarchy. I made a copy of it in the working directory and edited it there before exporting a new BUSCO_CONFIG_FILE path variable.
[isg4@farnam round1]$ cp $BUSCO/config.ini .
[isg4@farnam round1]$ export BUSCO_CONFIG_FILE=$PWD/config.iniFinally, we're ready to train Augustus. From Daren Card:
First, we must put together training sequences using the gene models we created in our first run of MAKER. We do this by issuing the following command to excise the regions that contain mRNA annotations based on our initial MAKER run (with 1000bp on each side).
Daren uses a series of awk commands to parse his .gff before passing it to bedtools, which extracts a .fasta file from the coordinates given. I'm splitting this command up for readability and debugging.
Command format:
bedtools getfasta -fi [input_fasta] -bed [output_bedfile] -fo [output_fasta][isg4@farnam P_amilis_R1.maker.output]$ mkdir -p Augustus/round1
[isg4@farnam P_amilis_R1.maker.output]$ cd Augustus/round1
[isg4@farnam round1]$ module load BEDTools
[isg4@farnam round1]$ awk -v OFS="\t" '{ if ($3 == "mRNA") print $1, $4, $5 }' ../P_amilis_R1.all.maker.noseq.gff | awk -v OFS="\t" '{ if ($2 < 1000) print $1, "0", $3+1000; else print $1, $2-1000, $3+1000 }' > P_amilis.mRNA1000.txt
[isg4@farnam round1]$ head P_amilis.mRNA1000.txt
Scaffold_2854 0 1667
Scaffold_2993 69460 73795
Scaffold_2993 71785 76503
Scaffold_2993 5439 8586
Scaffold_2993 37880 44766
Scaffold_2993 76460 80062
Scaffold_2993 108811 114899
Scaffold_2993 90342 99079
Scaffold_2993 114848 118034
Scaffold_2993 121951 132264
[isg4@farnam round1]$ module load BEDTools/2.27.1-foss-2016b
[isg4@farnam round1]$ bedtools getfasta -fi ../../../../Data/portulaca_26Nov2018_oK3Ko.fixed.fasta -bed P_amilis.mRNA1000.txt -fo P_amilis.mRNA1000.fastaThere will likely be some warnings like
Feature (Scaffold_3845:0-1981) beyond the length of Scaffold_3845 size (1032 bp). Skipping.that indicate that a certain feature couldn't be extracted because $$\pm$$1000bp on either end extends past the end of the scaffold. We could have more cleverly selected fragments to avoid this, but this is just for training Augustus, and as long as we end up with 100s of features we should be fine. Now, as Daren Card notes, we have a .fasta but no associated gene model information, which is what Augustus needs to train:
While we've only provided sequences from regions likely to contain genes, we've totally eliminated any existing annotation data about the starts/stops of gene elements. Augustus would normally use this as part of the training process. However, BUSCO will essentially do a reannotation of these regions using BLAST and built-in HMMs for a set of conserved genes (hundreds to thousands). This has the effect of recreating some version of our gene models for these conserved genes. We then leverage the internal training that BUSCO can perform (the
--longargument) to optimize the HMM search model to train Augustus and produce a trained HMM for MAKER.
Using run_BUSCO.py here additionally tells us some information about the quality of our initial MAKER2 annotation. We know our original BUSCO metrics from running our raw scaffolds through gVolante. Now we can compare how many of those have moderately high quality evidence (AED Augustus training set.
Command format:
python run_BUSCO.py -i [input_fasta] -c [cores] -o [output_prefix] -m [analysis_mode] -l [taxon_database] --long -sp [starting_species_parameters] -z --augustus_parameters='--progress=true'-
--long: use the initial gene models it creates to optimize the HMM settings -
-sp: species for initialAugustusparameterization. I used tomato (Solanum lycopersicum) because it is the closest relative to P. amilis listed in theAugustus README. -
-z: create tarball of output -
--augustus_parameters=='--progress=true': show a progressmeter
This is giving me trouble as a batch file so I'll go the tmux+srun route. For some reason, which I haven't yet looked into, this order of operations works, but the batch file wouldn't.
[isg4@farnam round1]$ tmux new -s busco
[isg4@farnam round1]$ srun --pty -p interactive bash
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam round1]$ module load BUSCO/3.0.2-foss-2016b-Python-2.7.13
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam round1]$ export BUSCO=/ysm-gpfs/apps/software/BUSCO/3.0.2-foss-2016b-Python-2.7.13/
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam round1]$ module load SAMtools/1.9-foss-2016b BamTools/2.5.1-foss-2016b HTSlib/1.5-foss-2016b BCFtools/1.5-foss-2016b Augustus
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam round1]$ export AUGUSTUS_CONFIG_PATH=~/project/apps/Augustus/config/
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam round1]$ export BUSCO_CONFIG_FILE=$PWD/config.ini
[isg4@farnam round1]$ python $BUSCO/bin/run_BUSCO.py -i P_amilis.mRNA1000.fasta -c 20 -o P_amilis_round1 -m genome -l $BUSCO/datasets/embryophyta_odb9 --long -sp tomato -z --augustus_parameters='--progress=true' 2>&1 | tee BUSCO.logYou'll see a -r flag in there that tells run_BUSCO.py to pick up from a run that did not complete. The way I typically work through pipelines like this is to build up to a working final command in an interactive bash session on our cluster and then submit the working command as a batch file. Including -r here is convenient because it allows the progress that run_BUSCO.py had already made to be utilized.
| BUSCO | Initial MAKER2 | Input scaffolds |
|---|---|---|
| Complete BUSCOs | 1031 (71.6%) | 1291 (89.7%) |
| Complete and single-copy BUSCOs | 814 (56.5%) | 1228 (85.3%) |
| Complete and duplicated BUSCOs | 217 (15.1%) | 63 (4.4%) |
| Fragmented BUSCOs | 158 (11.0%) | 29 (2.0%) |
| Missing BUSCOs | 251 (17.4%) | 120 (8.33%) |
| Total BUSCO groups searched | 1440 | 1440 |
It isn't surprising that our BUSCO scores were lower this time around; as I mentioned above, we only passed a subset of relatively high quality features to run_BUSCO.py+Augustus. These numbers are a bit on the low side though, and we'd ideally see the complete BUSCOs up around 80%. We also found almost four-fold as many duplicated BUSCOs than in our scaffolds. This is most likely due to isoforms that were annotated in MAKER2 but not filtered out before running run_BUSCO.py.
All of my relevant output files are in P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output/Augustus/round1/sample_data/run_P_amilis_round1; I'm not sure why they are in a directory called sample_data, which is empty except for run_P_amilis_round1, so I'll move it back a directory and get rid of sample_data. Then I'll rename the retraining parameters files (which look like BUSCO_P_amilis_round1_4212217375_*) in run_P_amilis_round1/augustus_output/retraining_parameters/ for readability. Some of the retraining files reference other retraining files, so we need to change the file names within those parameter files using sed. If you don't do this last step you'll get error from MAKER2 saying that Augustus can't find things like rate matrices, etc.
[isg4@farnam1 round1]$ mv sample_data/run_P_amilis_round1/ .
[isg4@farnam1 round1]$ rm -r sample_data
[isg4@farnam1 round1]$ cd run_P_amilis_round1/augustus_output/retraining_parameters/
[isg4@farnam1 retraining_parameters]$ rename "BUSCO_P_amilis_round1_4212217375" "Portulaca_amilis" *
[isg4@farnam1 retraining_parameters]$ sed -i 's/BUSCO_P_amilis_round1_4212217375/Portulaca_amilis/g' *Now we need to copy the retraining files into our Augustus config directory. Specifically, we're going to create a directory for a new species, just as there are pre-loaded parameters for dozens of species that come with Augustus; these directories live in Augustus/config/species.
[isg4@farnam1 retraining_parameters]$ export AUGUSTUS_CONFIG_PATH=/home/isg4/project/apps/Augustus/config
[isg4@farnam1 retraining_parameters]$ mkdir $AUGUSTUS_CONFIG_PATH/species/Portulaca_amilis
(base) [isg4@farnam1 retraining_parameters]$ cp * $AUGUSTUS_CONFIG_PATH/species/Portulaca_amilis- Figure out the order of operations with
Augustus--overall this whole section feels sloppy and hard to reproduce
With parameters for ab initio gene prediction in hand, we're ready for another round of MAKER2. This should consist of simply pointing MAKER2 towards our previously created .gff file from the initial run, our SNAP .hmm file, and the name of taxon for Augustus parameterization.
[isg4@farnam1 MAKER_single_core]$ mkdir Round2
[isg4@farnam1 MAKER_single_core]$ cp *.ctl Round2/
[isg4@farnam1 MAKER_single_core]$ cd Round2/
[isg4@farnam1 MAKER_single_core]$ mv maker_opts.R1.ctl maker_opts.R2.ctlI made sure to save my entire environment set up this time so that I can run MAKER2 in the future without having to figure out how I set it up originally; you can find copies of these in the Appendix
[isg4@c18n12 Round2]$ module save MAKER
[isg4@c18n12 Round2]$ env > env_MAKER.txt
[isg4@c18n12 Round2]$ less env_MAKER.txt
[isg4@c18n12 Round2]$ conda list > conda_MAKER.txtIn his tutorial, Daren Card said that MAKER2 wasn't recognizing his .gff from the original run, so he instead lifted the EST, protein, and repeat alignments using the commands below. Then he substituted them in the appropriate places in the MAKER2 opts file. I hit the same problem, and my initial run of MAKER2 with ab initio gene prediction found no proteins or ESTs, so I followed his lead.
# transcript alignments
[isg4@farnam1 P_amilis_R1.maker.output]$ awk '{ if ($2 == "est2genome") print $0 }' P_amilis_R1.all.maker.noseq.fixed.gff > P_amilis_R1.all.maker.est2genome.fixed.gff
# protein alignments
[isg4@farnam1 P_amilis_R1.maker.output]$ awk '{ if ($2 == "protein2genome") print $0 }' P_amilis_R1.all.maker.noseq.fixed.gff > P_amilis_R1.all.maker.protein2genome.fixed.gff
# repeat alignments
[isg4@farnam1 P_amilis_R1.maker.output]$ awk '{ if ($2 ~ "repeat") print $0 }' P_amilis_R1.all.maker.noseq.fixed.gff > P_amilis_R1.all.maker.repeats.fixed.gffNow to edit maker_opts.R2.ctl as stated above. Here is my new opts file. You'll notice a few things have changed for me that hopefully will stay the same for you; namely, I'm now pointing to a copy of my genome .fasta file with fixed scaffold names.
#-----Genome (these are always required)
genome=/home/isg4/scratch60/P_amilis_genome/Data/portulaca_26Nov2018_oK3Ko.fixed.fasta #genome sequence (fasta file or fasta embeded in GFF3 file)
organism_type=eukaryotic #eukaryotic or prokaryotic. Default is eukaryotic
#-----Re-annotation Using MAKER Derived GFF3
maker_gff=#MAKER derived GFF3 file
est_pass=0 #use ESTs in maker_gff: 1 = yes, 0 = no
altest_pass=0 #use alternate organism ESTs in maker_gff: 1 = yes, 0 = no
protein_pass=0 #use protein alignments in maker_gff: 1 = yes, 0 = no
rm_pass=0 #use repeats in maker_gff: 1 = yes, 0 = no
model_pass=0 #use gene models in maker_gff: 1 = yes, 0 = no
pred_pass=0 #use ab-initio predictions in maker_gff: 1 = yes, 0 = no
other_pass=0 #passthrough anyything else in maker_gff: 1 = yes, 0 = no
#-----EST Evidence (for best results provide a file for at least one)
est= #set of ESTs or assembled mRNA-seq in fasta format
altest= #EST/cDNA sequence file in fasta format from an alternate organism
est_gff=/home/isg4/scratch60/P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output/P_amilis_R1.all.maker.est2genome.fixed.gff #aligned ESTs or mRNA-seq from an external GFF3 file
altest_gff= #aligned ESTs from a closly relate species in GFF3 format
#-----Protein Homology Evidence (for best results provide a file for at least one)
protein= #protein sequence file in fasta format (i.e. from mutiple oransisms)
protein_gff= #aligned protein homology evidence from an external GFF3 file
#-----Repeat Masking (leave values blank to skip repeat masking)
model_org= #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff=/home/isg4/scratch60/P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output/P_amilis_R1.all.maker.repeats.fixed.gff #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)
#-----Gene Prediction
snaphmm=/home/isg4/scratch60/P_amilis_genome/MAKER_single_core/P_amilis_R1.maker.output/SNAP/round1/P_amilis.length50.aed0.25.hmm #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species=Portulaca_amilis #Augustus gene prediction species model
fgenesh_par_file= #FGENESH parameter file
pred_gff= #ab-initio predictions from an external GFF3 file
model_gff= #annotated gene models from an external GFF3 file (annotation pass-through)
est2genome=0 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=0 #infer predictions from protein homology, 1 = yes, 0 = no
trna=1 #find tRNAs with tRNAscan, 1 = yes, 0 = no
snoscan_rrna= #rRNA file to have Snoscan find snoRNAs
unmask=0 #also run ab-initio prediction programs on unmasked sequence, 1 = yes, 0 = no
#-----Other Annotation Feature Types (features MAKER doesn't recognize)
other_gff= #extra features to pass-through to final MAKER generated GFF3 file
#-----External Application Behavior Options
alt_peptide=C #amino acid used to replace non-standard amino acids in BLAST databases
cpus=1 #max number of cpus to use in BLAST and RepeatMasker (not for MPI, leave 1 when using MPI)
#-----MAKER Behavior Options
max_dna_len=100000 #length for dividing up contigs into chunks (increases/decreases memory usage)
min_contig=1 #skip genome contigs below this length (under 10kb are often useless)
pred_flank=200 #flank for extending evidence clusters sent to gene predictors
pred_stats=1 #report AED and QI statistics for all predictions as well as models
AED_threshold=1 #Maximum Annotation Edit Distance allowed (bound by 0 and 1)
min_protein=0 #require at least this many amino acids in predicted proteins
alt_splice=1 #Take extra steps to try and find alternative splicing, 1 = yes, 0 = no
always_complete=0 #extra steps to force start and stop codons, 1 = yes, 0 = no
map_forward=0 #map names and attributes forward from old GFF3 genes, 1 = yes, 0 = no
keep_preds=0 #Concordance threshold to add unsupported gene prediction (bound by 0 and 1)
split_hit=10000 #length for the splitting of hits (expected max intron size for evidence alignments)
single_exon=0 #consider single exon EST evidence when generating annotations, 1 = yes, 0 = no
single_length=250 #min length required for single exon ESTs if 'single_exon is enabled'
correct_est_fusion=0 #limits use of ESTs in annotation to avoid fusion genes
tries=2 #number of times to try a contig if there is a failure for some reason
clean_try=0 #remove all data from previous run before retrying, 1 = yes, 0 = no
clean_up=0 #removes theVoid directory with individual analysis files, 1 = yes, 0 = no
TMP= #specify a directory other than the system default temporary directory for temporary filesAs noted in Daren Card's Boa walk though, MAKER2 should ideally be run multiple times to improve gene models. This involves combining the MAKER2 output from the last round (section 2.4.2) retraining SNAP and Augustus (section 2.5) to improve predictions, and then running MAKER2 again (section 2.6). I also had a problem with an outdated install of tRNA-scan (I was using 1.4 but MAKER2 needs 2.x), so I needed to rerun MAKER2 to pick out tRNAs.
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@c18n12 Round2]$ cd P_amilis_R2.maker.output
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@c23n09 P_amilis_R2.maker.output]$ export MAKERBIN=/ysm-gpfs/project/jsi3/maker2/maker/bin
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@farnam1 MAKER_single_core]$ mkdir Round3
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@c23n09 MAKER_single_core]$ mkdir -p Round2/P_amilis_R2.maker.output/snap/round2
(/gpfs/ysm/project/isg4/conda_envs/genome) [isg4@c23n09 MAKER_single_core]$ cd Round2/P_amilis_R2.maker.output/snap/round2
After I finished one pass with MAKER2 and moved forward with a new set of transcriptomes and funannoate (see Genome-annotation-with-FUNannotate) I made a compressed archive of this work.
(base) [isg4@farnam2 project]$ tar -czvf P-amilis-MAKER2.tar.gz P_amilis_genome/
(base) [isg4@farnam2 project]$ rm -r P_amilis_genomeFor pieces of software that link multiple underlying analyses I'll only provide references to those pieces of software that I used, installed, and/or mentioned directly.
- Abrusán G., Grundmann N., DeMester L., Makalowski W. 2009. TEclass--a tool for automated classification of unknown eukaryotic transposable elements. Bioinformatics. 25:1329–1330. | TEclass
- Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc | FastQC
- Bao W., Kojima K.K., Kohany O. 2015. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 6:11.
- Benson G. 1999. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 27:573–580. | Tandem repeats finder
- Besemer J., Borodovsky M. 2005. GeneMark: web software for gene finding in prokaryotes, eukaryotes and viruses. Nucleic Acids Res. 33:W451–4. | GeneMark
- Bolger A.M., Lohse M., Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30:2114–2120. | Trimmomatic
- Camacho C., Coulouris G., Avagyan V., Ma N., Papadopoulos J., Bealer K., Madden T.L. 2009. BLAST+: architecture and applications. BMC Bioinformatics. 10:421. | BLAST+
- Chan, P.P. and Lowe, T. M. (2019) tRNAscan-SE: Searching for tRNA Genes in Genomic Sequences. Methods Mol Biol. 1962:1-14. | tRNA-scan 2.0
- Grabherr M.G., Haas B.J., Yassour M., Levin J.Z., Thompson D.A., Amit I., Adiconis X., Fan L., Raychowdhury R., Zeng Q., Chen Z., Mauceli E., Hacohen N., Gnirke A., Rhind N., di Palma F., Birren B.W., Nusbaum C., Lindblad-Toh K., Friedman N., Regev A. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 29:644–652. | Trinity
- Gurevich A., Saveliev V., Vyahhi N., Tesler G. 2013. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 29:1072–1075. | QUAST
- Holt C., Yandell M. 2011. MAKER2: an annotation pipeline and genome- database management tool for second- generation genome projects. BMC Bioinformatics. 12:491. | MAKER2
- Hubley R., Finn R.D., Clements J., Eddy S.R., Jones T.A., Bao W., Smit A.F.A., Wheeler T.J. 2016. The Dfam database of repetitive DNA families. Nucleic Acids Res. 44:D81–9.
- Kim D., Langmead B., Salzberg S.L. 2015. HISAT: a fast spliced aligner with low memory requirements. Nature Methods. 12:357–360. | HISAT
- Konig S., Lars Romoth, and Mario Stanke. (2017). Comparative Genome Annotation. In Setubal J.C., Stoye J., and Stadler P.F. Comparative Genomics. New York, NY: Humana Press.
- Korf I. 2004. Gene finding in novel genomes. BMC Bioinformatics. 5:59. | SNAP
- Krzywinski M., Schein J., Birol I., Connors J., Gascoyne R., Horsman D., Jones S.J., Marra M.A. 2009. Circos: an information aesthetic for comparative genomics. Genome Research. 19:1639–1645. | Circos
- Langmead B., Salzberg S.L. 2012. Fast gapped-read alignment with Bowtie 2. Nature Methods. 9:357–359. | Bowtie 2
- Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R., 1000 Genome Project Data Processing Subgroup. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 25:2078–2079. | SAMtools
- Matasci N., Hung L.-H., Yan Z., Carpenter E.J., Wickett N.J., Mirarab S., Nguyen N., Warnow T., Ayyampalayam S., Barker M., Burleigh J.G., Gitzendanner M.A., Wafula E., Der J.P., dePamphilis C.W., Roure B., Philippe H., Ruhfel B.R., Miles N.W., Graham S.W., Mathews S., Surek B., Melkonian M., Soltis D.E., Soltis P.S., Rothfels C., Pokorny L., Shaw J.A., DeGironimo L., Stevenson D.W., Villarreal J.C., Chen T., Kutchan T.M., Rolf M., Baucom R.S., Deyholos M.K., Samudrala R., Tian Z., Wu X., Sun X., Zhang Y., Wang J., Leebens-Mack J., Wong G.K.-S. 2014. Data access for the 1,000 Plants (1KP) project. Gigascience. 3:17.
- Nishimura O., Hara Y., Kuraku S. 2017. gVolante for standardizing completeness assessment of genome and transcriptome assemblies. Bioinformatics. 33:3635–3637. | gVolante
- Simão F.A., Waterhouse R.M., Ioannidis P., Kriventseva E.V., Zdobnov E.M. 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 31:3210–3212. | BUSCO
- Smit A. F. A., R. Hubley, P. Green. RepeatMasker Open-4.0. 2013-2015 http://www.repeatmasker.org. | Repeat masker
- Smit A. F. A., R. Hubley. RepeatModeler Open-1.0. 2008-2015 http://www.repeatmasker.org. | Repeat modeler
- Smith-Unna R., Boursnell C., Patro R., Hibberd J.M., Kelly S. 2016. TransRate: reference-free quality assessment of de novo transcriptome assemblies. Genome Research. 26:1134–1144. | TransRate
- Song L., Florea L. 2015. Rcorrector: efficient and accurate error correction for Illumina RNA-seq reads. Gigascience. 4:48. | Rcorrector
- Stanke M., Diekhans M., Baertsch R., Haussler D. 2008. Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics. 24:637–644. | Augustus
- Trapnell C., Williams B.A., Pertea G., Mortazavi A., Kwan G., van Baren M.J., Salzberg S.L., Wold B.J., Pachter L. 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnol. 28:511–515. | Cufflinks
- Wang B., Tseng E., Regulski M., Clark T.A., Hon T., Jiao Y., Lu Z., Olson A., Stein J.C., Ware D. 2016. Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing. Nature Communications. 7:1–13.
- Yandell M., Ence D. 2012. A beginner's guide to eukaryotic genome annotation. Nat Rev Genet. 13:1–14.
- Yang Y., Smith S.A. 2014. Orthology inference in nonmodel organisms using transcriptomes and low-coverage genomes: improving accuracy and matrix occupancy for phylogenomics. Mol. Biol. Evol. 31:3081–3092. | Phylogenomic dataset construction
# Name Version Build Channel
backcall 0.1.0 py37_0
blast 2.5.0 hc0b0e79_3 bioconda
boost 1.57.0 4
bzip2 1.0.6 h14c3975_5
ca-certificates 2019.1.23 0
certifi 2019.3.9 py37_0
decorator 4.4.0 py37_1
expat 2.2.5 hf484d3e_1002 conda-forge
fontconfig 2.13.1 he4413a7_1000 conda-forge
freetype 2.10.0 he983fc9_0 conda-forge
genometools-genometools 1.5.10 h470a237_1 bioconda
gettext 0.19.8.1 hd7bead4_3
giflib 5.1.7 h516909a_1 conda-forge
icu 58.2 h9c2bf20_1
infernal 1.1.2 h14c3975_2 bioconda
ipython 7.4.0 py37h39e3cac_0
ipython_genutils 0.2.0 py37_0
jedi 0.13.3 py37_0
jpeg 9c h14c3975_1001 conda-forge
libedit 3.1.20181209 hc058e9b_0
libffi 3.2.1 hd88cf55_4
libgcc-ng 8.2.0 hdf63c60_1
libgd 2.2.5 h0d07dcb_1005 conda-forge
libiconv 1.15 h516909a_1005 conda-forge
libpng 1.6.36 h84994c4_1000 conda-forge
libstdcxx-ng 8.2.0 hdf63c60_1
libtiff 4.0.10 h648cc4a_1001 conda-forge
libuuid 2.32.1 h14c3975_1000 conda-forge
libwebp 1.0.2 h99fbfcb_2 conda-forge
libxml2 2.9.8 h143f9aa_1005 conda-forge
ncurses 6.1 he6710b0_1
openssl 1.1.1b h7b6447c_1
parso 0.3.4 py37_0
perl 5.26.2 h470a237_0 conda-forge
pexpect 4.6.0 py37_0
pickleshare 0.7.5 py37_0
pip 19.0.3 py37_0
prompt_toolkit 2.0.9 py37_0
ptyprocess 0.6.0 py37_0
pygments 2.3.1 py37_0
python 3.7.3 h0371630_0
readline 7.0 h7b6447c_5
rmblast 2.2.28 h21aa3a5_4 bioconda
setuptools 40.8.0 py37_0
six 1.12.0 py37_0
snoscan 0.9b pl526ha92aebf_2 bioconda
sqlite 3.27.2 h7b6447c_0
tk 8.6.8 hbc83047_0
traitlets 4.3.2 py37_0
wcwidth 0.1.7 py37_0
wheel 0.33.1 py37_0
xz 5.2.4 h14c3975_4
zlib 1.2.11 h7b6447c_3