Genome annotation with funannotate
I opted to use funannotate, over MAKER for the final P. amilis genome annotation because it is much more user-friendly. At the time of the genome annotation, funannotate was not available on conda and we had a few hiccups running it on our HPC. Support from the YCRC's Ben Evans was key in a working installation of funannotate, which I access through the conda environment super_funannotate. This environment contained funannotate v1.6.0. Since that time the funannotate developers have created a conda package that is easily accessible.
To test to see if all of the components of funannotate are installed and working, try funannotate check --show-versions. When I did this I had and error stating that gmes_petap.pl couldn't be found. This is because the GeneMark suite of perl modules has a hardcoded shebang that calls /usr/bin/perl. All of the perl modules can be changed by using the change_path_in_perl_scripts.pl script.
(super_funannotate) [isg4@c14n10 ~]$ cd super_funannotate/funannotate_deps/gmes_petap
(super_funannotate) [isg4@c14n10 gmes_petap]$ perl change_path_in_perl_scripts.pl "super_funannotate/bin/perl"Genomic scaffolds will be made available through Phytozome upon release of our manuscript, but I'll begin with the data given to us by Dovetail Genomics.
(super_funannotate) [isg4@c15n09 scratch60]$ mkdir Pamilis_funannotate
(super_funannotate) [isg4@c15n09 scratch60]$ cd Pamilis_funannotate/
(super_funannotate) [isg4@c15n09 Pamilis_funannotate]$ mkdir dataAfter adding all the files, here is what data contains.
(super_funannotate) [isg4@farnam1 data]$ ls -lah
36G Apr 30 2018 DTG-SG-110_R1_001.fastq.gz
39G Apr 30 2018 DTG-SG-110_R2_001.fastq.gz
32G May 2 2018 DTG-SG-112-NS_S1_R1_001.fastq.gz
33G May 2 2018 DTG-SG-112-NS_S1_R2_001.fastq.gz
30G Apr 22 14:21 lib_001.sorted.md.bam
53G Apr 22 14:26 lib_002.sorted.md.bam
1.3K Apr 22 10:43 manifest.txt
181 Apr 22 10:43 md5sums.txt
390M Apr 22 10:43 Pamilis_denovo.fasta
1.7K Apr 22 10:43 Pamilis_denovo.input_breaks.txt
1.4M Apr 22 10:43 Pamilis_denovo.table.txt
390M Apr 22 10:43 Pamilis_Final_HiRise.fasta
2.4K Apr 22 10:43 Pamilis_Final_HiRise.input_breaks.txt
321K Apr 22 10:43 Pamilis_Final_HiRise.table.txtDovetail includes this description of some files in manifest.txt.
Portulaca Manifest | November 26, 2018
-
Pamilis_Final_HiRise.fasta.gz- FASTA file containing the final Hirise scaffolds.
-
Pamilis_Final_HiRise.input_breaks.txt- Tab-delimited table describing positions of breaks made in the input assembly scaffolds. Follows the same format as the table file described below.
-
Pamilis_Final_HiRise.report.pdf- Hirise scaffolding report, contains Dovetail library insert distribution(s), and contiguity plots and statistics for the input and output assemblies.
-
Pamilis_Final_HiRise.table.txt- Tab-delimited table describing positions of input assembly scaffolds in the Hirise scaffolds. The table has the following format:
- HiRise scaffold name
- Input sequence name
- Starting base (zero-based) of the input sequence
- Ending base of the input sequence
- Strand (- or +) of the input sequence in the scaffold, where '-' in the strand column indicates that the sequence is reverse complemented relative to the input assembly
- Starting base (zero-based) in the HiRise scaffold
- Ending base in the HiRise scaffold
- Tab-delimited table describing positions of input assembly scaffolds in the Hirise scaffolds. The table has the following format:
-
lib_001.sorted.md.bam- Sorted
bamfile of the original reads
- Sorted
-
lib_002.sorted.md.bam- Sorted
bamfile of the HiRise reads
- Sorted
-
md5sums.txt- Hash table for
lib_001.sorted.md.bam,lib_002.sorted.md.bam, andPamilis_Final_HiRise.fasta.gz
- Hash table for
Note that I changed the original file names to make them more discernible for what they contain with rename "portulaca_26Nov2018_oK3Ko" "Pamilis_Final_HiRise" *.
I've taken excerpts from the authors of funannotate, and purge_haplotigs that describe some of the problems faced when taking an assembly in the form of contigs and preparing it for haploid genome annotation.
Problems
-
funannotate: When working with haploid assemblies, sometimes you want to remove some repetitive contigs that are contained in other scaffolds of the assembly. If the repeats are indeed unique, then we want to keep them in the assembly. -
purge_haplotigs: Some parts of a genome may have a very high degree of heterozygosity. This causes contigs for both haplotypes of that part of the genome to be assembled as separate primary contigs, rather than as a contig and an associated haplotig. This can be an issue for downstream analysis whether you're working on the haploid or phased-diploid assembly.
Solutions
-
funannotate: Use a “leave one out” method (minimap2 or mummer (nucmer)), where the the shortest contigs/scaffolds are aligned to the rest of the assembly to determine if it is repetitive. -
purge_haplotigs: Identify pairs of contigs that are syntenic and move one of them to the haplotig 'pool'.
I'm not sure how Dovetail did this during the assembly process, but the previous run-through of MAKER produced sensible results, so I'm not overly concerned with the downstream influence heterozygosity, especially since P. amilis is highly cleistogamous. Regardless, it would be ideal to understand the degree of heterozygosity and separate haplotigs if they are identifiable.
In my data directory I ran the follow batch file.
#!/bin/bash
#SBATCH --job-name=funannotate_clean
#SBATCH --time=120:00:00
#SBATCH --nodes=1
#SBATCH --cpus-per-task=12
#SBATCH --ntasks=1
#SBATCH --partition=general
#SBATCH --mem=60G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=funannotate_clean-%j.out
#SBATCH --error=funannotate_clean-%j.err
funannotate clean --input Pamilis_Final_HiRise.fasta --out Pamilis_Final_HiRise.clean.fastaNext I used the sort and rename function of funannotate to rename scaffolds that have names that make Augustus unhappy and sorted scaffolds by size. Then, to keep track of which files were generated by which steps, I put all of the clean-related files into a directory called CLEAN
(super_funannotate) [isg4@c20n01 data]$ funannotate sort -b scaffold -i Pamilis_Final_HiRise.clean.fasta -o Pamilis_Final_HiRise.clean.sorted.fasta
(super_funannotate) [isg4@c20n01 Pamilis_funannotate]$ mkdir -p ../CLEAN
(super_funannotate) [isg4@c20n01 Pamilis_funannotate]$ mv data/*clean* ../CLEANNow we have all our scaffolds cleaned and organized from largest to smallest. No scaffolds were removed using this strategy.
Following the purge_haplotigs walkthrough, I created a new conda environment so that I don't interfere with the super_funannotate environment.
(super_funannotate) [isg4@c20n01 data]$ conda activate base
(base) [isg4@c20n01 data]$ conda create -n purge_haplotigs
(base) [isg4@c20n01 data]$ conda activate purge_haplotigs
(purge_haplotigs) [isg4@c20n01 data]$ conda config --add channels defaults
(purge_haplotigs) [isg4@c20n01 data]$ conda config --add channels bioconda
(purge_haplotigs) [isg4@c20n01 data]$ conda config --add channels conda-forge
(purge_haplotigs) [isg4@c20n01 data]$ conda install purge_haplotigs
(purge_haplotigs) [isg4@c20n01 data]$ mkdir ../Purge_Haplotigs
(purge_haplotigs) [isg4@c20n01 data]$ cd ../Purge_HaplotigsThe initial step of purge_haplotigs is to generate a coverage histogram, which I did by recycling some files used for another analysis, Referee, which is not presented here. I first needed to create and index and map the genomic reads into cooridnate sorted files.
#!/bin/bash
#SBATCH --job-name=bwa_mem
#SBATCH --time=120:00:00
#SBATCH --nodes=1
#SBATCH --cpus-per-task=20
#SBATCH --ntasks=1
#SBATCH --partition=general
#SBATCH --mem=120G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=bwa_mem-%j.out
#SBATCH --error=bwa_mem-%j.err
bwa index Pamilis_Final_HiRise.clean.sorted.fasta
bwa mem Pamilis_Final_HiRise.clean.sorted.fasta ../data/DTG-SG-112-NS_S1_R1_001.fastq ../data/DTG-SG-112-NS_S1_R2_001.fastq > Pamilis_Final_HiRise.clean.sorted-reads-mapped.samThe resulting .sam file was over 300Gb! Checking out the contents of the .sam
(super_funannotate) [isg4@c13n11 Referee]$ head Pamilis_Final_HiRise.clean.sorted-reads-mapped.sam
@SQ SN:scaffold_1 LN:53436919
@SQ SN:scaffold_2 LN:45464428
@SQ SN:scaffold_3 LN:45200182
@SQ SN:scaffold_4 LN:44979772
@SQ SN:scaffold_5 LN:42597560
@SQ SN:scaffold_6 LN:42146318
@SQ SN:scaffold_7 LN:41375687
@SQ SN:scaffold_8 LN:41005119
@SQ SN:scaffold_9 LN:39183218
@SQ SN:scaffold_10 LN:560248showed it was properly formatted (first line didn't start with @HQ). Next onto .bam conversion, sorting, and creation of the pileup file.
(super_funannotate) [isg4@c13n11 Referee]$ samtools view -b Pamilis_Final_HiRise.clean.sorted-reads-mapped.sam > Pamilis_Final_HiRise.clean.reads-mapped.bam
(super_funannotate) [isg4@c19n05 Referee]$ samtools sort Pamilis_Final_HiRise.clean.reads-mapped.bam -o Pamilis_Final_HiRise.clean.sorted-reads-mapped.bam
(super_funannotate) [isg4@c19n05 Referee]$ samtools mpileup -d 999999999 -f Pamilis_Final_HiRise.clean.sorted.fasta -Q 0 -s -B -o Pamilis_Final_HiRise.clean.sorted-reads-mapped.pileup Pamilis_Final_HiRise.clean.sorted-reads-mapped.bamI used the resulting .bam to create the read histogram in purge_haplotigs.
#!/bin/bash
#SBATCH --job-name=readhist
#SBATCH --time=120:00:00
#SBATCH --nodes=1
#SBATCH --cpus-per-task=20
#SBATCH --ntasks=1
#SBATCH --partition=general
#SBATCH --mem=36G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=readhist-%j.out
#SBATCH --error=reahist-%j.err
purge_haplotigs readhist -b ../Referee/Pamilis_Final_HiRise.clean.sorted-reads-mapped.bam -g ../CLEAN/Pamilis_Final_HiRise.clean.sorted.fasta -t 20Taking a look at the coverage plot below it does not have a bimodal distribution, as expected when assembling reads from a diploid organism into a haploid genome, therefore I did not continue with purple_haplotigs, as this appears to be an accurate haploid representation.
I will use the RNAseq data generated by a CAM induction experiment (see the Portulaca-coexpression repo) for annotating the P. amilis genome. These data will be availe under the NCBI's BioProject PRJNA732408 upon publication of our manuscript.
CAM-induction experiment Plants were grown in a growth chamber with 14 hour days and a maximum PAR of ~385 (mol)(m^-2)(s^-1). Temperature and lighting were as follows:
- 00:00 lights off, 22.0°
- 06:00 lights on, 25°C
- 20:00 lights off, 22.0°
Plants were allowed to grow for 4-6 weeks until they had begun to flower, fruit, and had developed sufficient branching and leaves for repeated sampling. At which point RNA was sampled from six replicates at the following time points:
- 10:00, a few hours into the light period plants are rapidly photosynthesizing; CAM plants are decarboxylating over-night malate stores
- 16:00, in the late afternoon photosynthetic rates tend to slow, especially in CAM plants that have used much of their stored malate
- 22:00, a few hours into the dark period photosynthesis has halted and CAM plants ramp up malate production
- 04:00, malate production is showing in CAM plants
Water was then withheld for ~10 days and induction of CAM was monitored in a sample of designated indicator plants that were subject to repeated tests of nocturnal acid accumulation. After significant changes in nocturnal acid accumulation were found in the indicator samples, all test plants were considered to be exhibiting CAM, and again sampled at the above time points. This resulted in 4 time points x 6 replicates x 2 conditions (well-water vs. drought) = 48 samples. We also included two populations of P. oleracea in this experiment, but those data won't be discussed here.
Library preparation and sequencing was done at the Yale Center for Genome Analysis (YCGA). Our libraries were done with poly-A tailing and sequenced on an Illumina HiSeq (100bp PE) in pools such that each sample received ~25M reads.
FastQC analyses metrics like per base/sequence quality, GC content, N content, sequence length distribution, duplication levels, and adapter content. It nicely summarizes these in html files with visuals to discern between low, moderate, and high quality data.
(super_funannotate) [isg4@c14n03 drought-diff-exp]$ mkdir -p quality_control/FastQC
(super_funannotate) [isg4@c14n03 drought-diff-exp]$ cd quality_control/Then, in an array-formatted batch file, I ran my FastQC command:
fastqc -t 12 --outdir ./FastQC/ ../data/*/*.fastq.gzI found that my forward (_R1_) reads frequently had poly-A as an over-represented sequence.
I initially tried removing poly-A/T tails with Trimmomatic, but it wasn't very succesful. Instead, I found PRINSEQ worked well, although was less user friendly. Because PRINSEQ does not take compressed files, I first uncompressed the sequence data and then ran them through and array batch job, below.
#!/bin/bash
#SBATCH --job-name=prinseq
#SBATCH --time=24:00:00
#SBATCH -c 5
#SBATCH -p general
#SBATCH --mem=24G
#SBATCH --array=1-48
#SBATCH --ntasks-per-node=1
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=prinseq-%A-%a.out
#SBATCH --error=prinseq-%A-%a.err
module load PRINSEQ
sleep $((SLURM_ARRAY_TASK_ID*60))
SAMPLES=(
P-amilis_002_174_019
P-amilis_004_150_043
...
P-amilis_096_295_282
)
SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin $SAMPLE"
mkdir -p ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Unzipping"
for f in ~/scratch60/drought-diff-exp/data/${SAMPLE}/*.gz
do STEM=$(basename "${f}" .gz)
echo $STEM
gunzip -c "${f}" > ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/"${STEM}"
done
cd ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin PRINSEQ"
prinseq-lite.pl -verbose -fastq *_R1_001.fastq -fastq2 *_R2_001.fastq -out_good ${SAMPLE}.G -out_bad null -log ${SAMPLE}.log -graph_data ${SAMPLE}.gd -out_format 3 -min_qual_mean 25 -trim_tail_left 12 -trim_tail_right 12 -trim_qual_left 20 -trim_qual_right 20 -ns_max_p 1
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished PRINSEQ"
for g in ./*.gd
do STEM2=$(basename "${g}" .gd)
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin graphing $STEM2"
prinseq-graphs.pl -log -html_all -png_all -i ${STEM2}.gd -o ${STEM2}.gd
done
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished graphing"For annotation purposes, we're only concerned with capturing diversity of transcripts to provide evidence to gene models. To this end, I concatenated all of the .fastq files that made it through filter with PRINSEQ and use those for annotation. However, when these reads were input into Trinity v2.5.1 as part as funannotate train (below) I hit an issue with header formatting, which I solved by using the reformatting suggested by Brian Haas and implemented by izabelcavassim. I used her solution to reformat each set of reads in a batch script that cra the following code block.
for f in *amilis*/*.G_1.fastq
do
STEM=$(basename "${f}" .fastq)
DIR=$(basename "${STEM}" .G_1)
echo $STEM
echo $DIR
cat $f | sed '/^@/ s/ .*//' > ./${DIR}/${STEM}.RFT.fastq.tmp
awk '{ if (NR%4==1) { print $1""$2"/1" } else { print } }' ./${DIR}/${STEM}.RFT.fastq.tmp > ./${DIR}/${STEM}.RFT.fastq
rm ./${DIR}/${STEM}.RFT.fastq.tmp
cat ./${DIR}/${STEM}.RFT.fastq >> P-amilis-concat.RFT.G_1.fastq
done
for f in *amilis*/*.G_2.fastq
do
STEM=$(basename "${f}" .fastq)
DIR=$(basename "${STEM}" .G_2)
echo $STEM
echo $DIR
cat $f | sed '/^@/ s/ .*//' > ./${DIR}/${STEM}.RFT.fastq.tmp
awk '{ if (NR%4==1) { print $1""$2"/2" } else { print } }' ./${DIR}/${STEM}.RFT.fastq.tmp > ./${DIR}/${STEM}.RFT.fastq
rm ./${DIR}/${STEM}.RFT.fastq.tmp
cat ./${DIR}/${STEM}.RFT.fastq >> P-amilis-concat.RFT.G_2.fastq
doneI initially installed RepeatModeler, RepeatMasker, and dependencies with Dfam 2.0, so first I updated the database to 3.0 within the directory RepeatMasker/Libraries. Check the release notes for changes from 2.0 to 3.0.
(super_funannotate) [isg4@c20n01 Libraries]$ curl -O http://www.dfam.org/releases/Dfam_3.0/families/Dfam.hmm.gz
(super_funannotate) [isg4@c20n01 Libraries]$ curl -O http://www.dfam.org/releases/Dfam_3.0/families/Dfam.embl.gz
(super_funannotate) [isg4@c20n01 Libraries]$ gunzip Dfam.hmm.gz
(super_funannotate) [isg4@c20n01 Libraries]$ gunzip Dfam.embl.gzBack in my funannotate directory, I made a MASK working directory for repeat-related steps.
(super_funannotate) [isg4@c20n01 Pamilis_funannotate]$ mkdir MASK
(super_funannotate) [isg4@c20n01 Pamilis_funannotate]$ cd MASKFirst, I built the database for a custom RepeatModeler run.
(super_funannotate) [isg4@c20n01 MASK]$ BuildDatabase -name Portulaca_amilis -engine ncbi ../CLEAN/Pamilis_Final_HiRise.clean.sorted.fasta
>>> Building database Portulaca_amilis:
>>> Adding ../CLEAN/Pamilis_Final_HiRise.clean.sorted.fasta to database
>>> Number of sequences (bp) added to database: 4053 ( 403885173 bp )This produced the expected output files
<prefix>.nhr<prefix>.nnd<prefix>.nog<prefix>.translation<prefix>.nin<prefix>.nni<prefix>.nsq- something that looks like
RM_18698.WedMay80847052019
There was no update to RepeatModeler (still v1.0.11) but a new RepeatMasker (v4.0.9-p2) release contained new databases so I installed a local verion outside of the funannotate pipeline.
(base) [isg4@c13n11 ~]$ cd project/apps
(base) [isg4@c13n11 apps]$ mv RepeatMasker/ RepeatMasker4.0.9
(base) [isg4@c13n11 apps]$ curl -O http://www.repeatmasker.org/RepeatMasker-open-4-0-9-p2.tar.gz
(base) [isg4@c13n11 apps]$ gunzip RepeatMasker-open-4-0-9-p2.tar.gz
(base) [isg4@c13n11 apps]$ tar xvf RepeatMasker-open-4-0-9-p2.tar
(base) [isg4@c13n11 apps]$ mv RepeatMasker RepeatMasker4.0.9-p2
(base) [isg4@c13n11 apps]$ cd RepeatMasker4.0.9-p2/
(base) [isg4@c13n11 RepeatMasker4.0.9-p2]$ curl -O https://www.girinst.org/server/RepBase/protected/repeatmaskerlibraries/RepBaseRepeatMaskerEdition-20181026.tar.gz
(base) [isg4@c13n11 RepeatMasker4.0.9-p2]$ gunzip RepBaseRepeatMaskerEdition-20181026.tar.gz
(base) [isg4@c13n11 RepeatMasker4.0.9-p2]$ tar xvf RepBaseRepeatMaskerEdition-20181026.tar
(base) [isg4@c13n11 RepeatMasker4.0.9-p2]$ conda activate super_funannotate
(super_funannotate) [isg4@c13n11 RepeatMasker4.0.9-p2]$ perl ./configure
(super_funannotate) [isg4@c13n11 RepeatMasker4.0.9-p2]$ cd ../RepeatModeler
(super_funannotate) [isg4@c13n11 RepeatMasker4.0.9-p2]$ perl ./configureThen I built the database,
(super_funannotate) [isg4@c13n11 RepeatMasker4.0.9-p2]$ cd ~/scratch60/Pamilis_funannotate/MASK/
(super_funannotate) [isg4@c13n11 MASK]$ /home/isg4/project/apps/RepeatModeler/BuildDatabase -name Portulaca_amilis -engine ncbi ../CLEAN/Pamilis_Final_HiRise.clean.sorted.fasta
Building database Portulaca_amilis:
Reading ../CLEAN/Pamilis_Final_HiRise.clean.sorted.fasta...
Number of sequences (bp) added to database: 4053 ( 403885173 bp )and submitted RepeatModeler as a batch file.
#!/bin/bash
#SBATCH --job-name=repeat_modeler
#SBATCH --time=144:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=36G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=repeatmodeler-%j.out
#SBATCH --error=repeatmodeler-%j.err
#SBATCH --ntasks-per-node=1
/home/isg4/project/apps/RepeatModeler/RepeatModeler -pa 36 -engine ncbi -database Portulaca_amilisChecking the results,
(super_funannotate) [isg4@c19n05 RM_7749.ThuMay160205052019]$ grep ">" consensi.fa.classified | wc
2092 37625 240749
(super_funannotate) [isg4@c19n05 RM_7749.ThuMay160205052019]$ grep "Unknown" consensi.fa.classified | wc
1419 25515 162009~67.8% of our repeats were not classified. In my previous run of RepeatModeler as part of the MAKER pipeline I recovered slightly more repeat families (2130) but classified the same percentage of them. Next, I'll put the consensus file (consensi.fa.classified) into TEclass to aid in repeat family classification. TEclass classifies unknown transpsosable elements into DNA transposons, LTRs, LINEs, and SINEs, and identified ORFs.
-
ID: -
Title: RM_7749.ThuMay160205052019 -
File: consensi.fa.classified -
Text:
Resulting ID is TEclass 155809850939466, and stats are below.
| Total | 2089 |
|---|---|
| DNA transposons | 661 |
| LTRs | 617 |
| LINEs | 288 |
| SINEs | 71 |
| Unclear | 452 |
I downloaded the resulting file (RM_7749.ThuMay160205052019-lib.txt) and opened it in a text editor so that I could change the format of the sequence names to match RepeatModeler.
- find:
^>(\w*.*)#(Unknown)(.*)\|TEclass result: (\w+)\|.*; replace:>$1#$4 $3 - find:
unclear; replace:Unknown
Then I pushed it back to the cluster, and renamed it.
# LOCAL MACHINE
➜ RM_7749.ThuMay160205052019 $ rsync -av RM_7749.ThuMay160205052019-lib.txt isg4@farnam.hpc.yale.edu:/home/isg4/scratch60/Pamilis_funannotate/MASK/RM_7749.ThuMay160205052019/
# CLUSTER
(super_funannotate) [isg4@c19n05 RM_7749.ThuMay160205052019]$ mv RM_7749.ThuMay160205052019-lib.txt consensi.TEclass.fa.classifiedI first masked repeats using the custom library I made above and then used the Viridiplantae repeat library. I made a new directory for each analysis and run RepeatMasker in batch files.
(super_funannotate) [isg4@c19n05 MASK]$ mkdir Pamilis_MASK
(super_funannotate) [isg4@c19n05 MASK]$ cd Pamilis_MASK/Note: PASA cannot take hardmasked files, therefore RepeatMasker must be run with -xsmall, which softmasks the genome, rather than -x.
#!/bin/bash
#SBATCH --job-name=repeat_masker
#SBATCH --time=72:00:00
#SBATCH -c 16
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=repeatmasker-%j.out
#SBATCH --error=repeatmasker-%j.err
#SBATCH --ntasks-per-node=1
~/project/apps/RepeatMasker4.0.9-p2/RepeatMasker -a -gff -xsmall -pa 20 -e ncbi -lib ../RM_7749.ThuMay160205052019/consensi.TEclass.fa.classified -dir ./ ../../CLEAN/Pamilis_Final_HiRise.clean.sorted.fastaThe results are below, and follow very closely to the previous analysis.
==================================================
file name: Pamilis_Final_HiRise.clean.sorted.fasta
sequences: 4053
total length: 403885173 bp (400444544 bp excl N/X-runs)
GC level: 37.05 %
bases masked: 185395132 bp ( 45.90 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
SINEs: 21169 2585805 bp 0.64 %
ALUs 0 0 bp 0.00 %
MIRs 0 0 bp 0.00 %
LINEs: 112223 30401101 bp 7.53 %
LINE1 23651 10545765 bp 2.61 %
LINE2 363 192115 bp 0.05 %
L3/CR1 1054 921651 bp 0.23 %
LTR elements: 166932 70282446 bp 17.40 %
ERVL 0 0 bp 0.00 %
ERVL-MaLRs 0 0 bp 0.00 %
ERV_classI 191 48453 bp 0.01 %
ERV_classII 0 0 bp 0.00 %
DNA elements: 274638 56733054 bp 14.05 %
hAT-Charlie 227 259411 bp 0.06 %
TcMar-Tigger 0 0 bp 0.00 %
Unclassified: 93162 19071981 bp 4.72 %
Total interspersed repeats:179074387 bp 44.34 %
Small RNA: 1137 154279 bp 0.04 %
Satellites: 305 351303 bp 0.09 %
Simple repeats: 126037 5319059 bp 1.32 %
Low complexity: 17962 857397 bp 0.21 %Because the file names kept growing in length, I renamed all of the RepeatMasker output files with a truncated prefix and indicated that this is a custom library mask.
(super_funannotate) [isg4@c13n10 Pamilis_MASK]$ rename Pamilis_Final_HiRise.clean.sorted.fasta Pamilis_FHR.clean.sort.PAMmask.fasta *Here, FHR stands for Final_HiRise and PAMmask refers to the P. amilis custom library mask. I then used the output from the custom library repeat masking (Pamilis_FHR.clean.sort.PAMmask.fasta.masked) as the input into the next round of RepeatMasker with the Viridiplantae library.
(super_funannotate) [isg4@c13n10 MASK]$ mkdir Viridiplantae_MASK
(super_funannotate) [isg4@c13n10 MASK]$ cd Viridiplantae_MASK/I then altered the RepeatMasker command in the batch file above,
~/project/apps/RepeatMasker4.0.9-p2/RepeatMasker -a -gff -xsmall -pa 20 -e ncbi -species viridiplantae -dir ./ ../Pamilis_MASK/Pamilis_FHR.clean.sort.PAMmask.fasta.maskedand ran it. Once it finished I renamed the results.
(super_funannotate) [isg4@c13n10 Viridiplantae_MASK]$ rename PAMmask VIRmask *
(super_funannotate) [isg4@c13n10 Viridiplantae_MASK]$ rename .masked. . *We can see that very little additional masking was done.
==================================================
file name: Pamilis_FHR.clean.sort.PAMmask.fasta.masked
sequences: 4053
total length: 403885173 bp (215249207 bp excl N/X-runs)
GC level: 37.28 %
bases masked: 2760236 bp ( 0.68 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
Retroelements 8009 1742061 bp 0.43 %
SINEs: 42 3621 bp 0.00 %
Penelope 1 101 bp 0.00 %
LINEs: 2050 322990 bp 0.08 %
CRE/SLACS 1 45 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 90 12702 bp 0.00 %
L1/CIN4 1954 309936 bp 0.08 %
LTR elements: 5917 1415450 bp 0.35 %
BEL/Pao 0 0 bp 0.00 %
Ty1/Copia 3319 912973 bp 0.23 %
Gypsy/DIRS1 2358 449750 bp 0.11 %
Retroviral 0 0 bp 0.00 %
DNA transposons 3471 749055 bp 0.19 %
hobo-Activator 897 281623 bp 0.07 %
Tc1-IS630-Pogo 60 23185 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 387 95104 bp 0.02 %
Other (Mirage, 2 95 bp 0.00 %
P-element, Transib)
Rolling-circles 559 111953 bp 0.03 %
Unclassified: 89 13760 bp 0.00 %
Total interspersed repeats: 2504876 bp 0.62 %
Small RNA: 976 112690 bp 0.03 %
Satellites: 50 2964 bp 0.00 %
Simple repeats: 467 23079 bp 0.01 %
Low complexity: 115 6305 bp 0.00 %
==================================================Now we'll combine the results from the two repeat libraries and polish them with ProcessRepeats. From Sebastian Tempel
ProcessRepeatsfollows 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 ofProcessRepeatsalgorithm is to merge fragmented repeats and to assign appropriate subfamily names to the search engine output.
This process involves copying the final repeatmasker output .fasta and .out files to a new directory, unzipping the .cat mask files, and combining the unzipped .cat files, which ProcessRepeats takes as an input. After that's done we can parse the repeats to build a repeat landscape.
(super_funannotate) [isg4@c13n10 MASK]$ mkdir Final_MASK
(super_funannotate) [isg4@c13n10 MASK]$ cd Final_MASK/
(super_funannotate) [isg4@c13n10 Final_MASK]$ cp ../Viridiplantae_MASK/Pamilis_FHR.clean.sort.VIRmask.fasta.masked ./Pamilis_FHR.clean.sort.FINALmask.fasta
(super_funannotate) [isg4@c13n10 Final_MASK]$ cp ../Viridiplantae_MASK/Pamilis_FHR.clean.sort.VIRmask.fasta.out ./Pamilis_FHR.clean.sort.FINALmask.out
(super_funannotate) [isg4@c13n10 Final_MASK]$ gunzip ../Viridiplantae_MASK/*.cat.gz ../Pamilis_MASK/*.cat.gz
(super_funannotate) [isg4@c13n10 Final_MASK]$ cat ../Viridiplantae_MASK/*.cat ../Pamilis_MASK/*.cat > Pamilis_FHR.clean.sort.FINALmask.catFinally, I input the .cat files into ProcessRepeats.
#!/bin/bash
#SBATCH --job-name=process_repeats
#SBATCH --time=48:00:00
#SBATCH -c 16
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=processrepeats-%j.out
#SBATCH --error=processrepeats-%j.err
#SBATCH --ntasks-per-node=1
ProcessRepeats -species viridiplantae Pamilis_FHR.clean.sort.FINALmask.cat -aRemember to include the -a so ProcessRepeats outputs a .align file for use in parseRM below. The final output looks almost identical to the original set of RepeatMasker runs.
==================================================
file name: Pamilis_FHR.clean.sort.FINALmask
sequences: 4053
total length: 403885173 bp (400656769 bp excl N/X-runs)
GC level: Unknown %
bases masked: 189675681 bp ( 46.96 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
Retroelements 308416 104987661 bp 25.99 %
SINEs: 21212 2589408 bp 0.64 %
Penelope 313 310102 bp 0.08 %
LINEs: 114290 30723714 bp 7.61 %
CRE/SLACS 1 45 bp 0.00 %
L2/CR1/Rex 1417 1113766 bp 0.28 %
R1/LOA/Jockey 0 0 bp 0.00 %
R2/R4/NeSL 2427 336211 bp 0.08 %
RTE/Bov-B 25497 4664931 bp 1.16 %
L1/CIN4 25604 10855295 bp 2.69 %
LTR elements: 172914 71674539 bp 17.75 %
BEL/Pao 1253 454421 bp 0.11 %
Ty1/Copia 33964 15102607 bp 3.74 %
Gypsy/DIRS1 47222 32995813 bp 8.17 %
Retroviral 191 48453 bp 0.01 %
DNA transposons 278077 57477188 bp 14.23 %
hobo-Activator 35562 8219511 bp 2.04 %
Tc1-IS630-Pogo 37369 6253034 bp 1.55 %
En-Spm 0 0 bp 0.00 %
MuDR-IS905 0 0 bp 0.00 %
PiggyBac 0 0 bp 0.00 %
Tourist/Harbinger 6678 1926232 bp 0.48 %
Other (Mirage, 2 95 bp 0.00 %
P-element, Transib)
Rolling-circles 0 0 bp 0.00 %
Unclassified: 99562 20464466 bp 5.07 %
Total interspersed repeats: 182929315 bp 45.29 %
Small RNA: 2117 266904 bp 0.07 %
Satellites: 356 354267 bp 0.09 %
Simple repeats: 126494 5341955 bp 1.32 %
Low complexity: 18049 863668 bp 0.21 %
==================================================Now we'll do a little post-processing to get them ready for downstream analyses.
(super_funannotate) [isg4@c13n10 Final_MASK]$ ~/project/apps/RepeatMasker4.0.9-p2/util/rmOutToGFF3.pl Pamilis_FHR.clean.sort.FINALmask.out > Pamilis_FHR.clean.sort.FINALmask.out.gff3
(super_funannotate) [isg4@c13n10 Final_MASK]$ grep -v -e "Satellite" -e ")n" -e "-rich" Pamilis_FHR.clean.sort.FINALmask.out.gff3 > Pamilis_FHR.clean.sort.FINALmask.complex.out.gff3
(super_funannotate) [isg4@c13n10 Final_MASK]$ cat Pamilis_FHR.clean.sort.FINALmask.complex.out.gff3 | perl -ane '$id; if(!/^\#/){@F = split(/\t/, $_); chomp $F[-1];$id++; $F[-1] .= "\;ID=$id"; $_ = join("\t", @F)."\n"} print $_' > Pamilis_FHR.clean.sort.FINALmask.complex.out.reformat.gff3And finally I'll run parseRM to generate the data for a repeat landscape.
#!/bin/bash
#SBATCH --job-name=parseRM
#SBATCH --time=48:00:00
#SBATCH -p bigmem
#SBATCH -N 1
#SBATCH -c 16
#SBATCH --mem 200G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=parseRM-%j.out
#SBATCH --error=parseRM-%j.err
perl ~/project/apps/Parsing-RepeatMasker-Outputs/parseRM.pl -i Pamilis_FHR.clean.sort.FINALmask.align -g 403885173 --parse
perl ~/project/apps/Parsing-RepeatMasker-Outputs/parseRM.pl -i Pamilis_FHR.clean.sort.FINALmask.align --land 50,1 -vI created some code to create repeat landscapes that is available in the Genome statistics notebook. Here's the Portulaca amilis repeat landscape, where the main figure shows the repeat families and the inset shows the broader repeat classes. The top left color in the legends refers to the top color of each stacked bar and proceed down each legend column from left to right. I find the easiest portion to see how the legend reads is from MULE-MuDR to Gypsy.
The purpose of most genome annotation projects is ultimately, on some level, to produce a functional genome annotation; that is, to identify and name homologous genomic features. In particular, gene models are often what we are most interested in. Therefore, choosing which evidence to provide to the annotation pipeline is one of the most important decisions a user makes. At a minimum, a researcher should have a transcriptome from the same organism, but the addition of ESTs and protein models will greatly improve final gene models.
Most non-model systems will not have EST data or sets of manually curated proteins for annotation. It's tempting to use the giant repositories online to pull ESTs and protein models from other organisms, and many papers, walkthroughs, and software docs say it's fine to use data from "closely related organisms". Unfortunately, virtually no one states what "closely related" means and what the short comings are of using less closely related organisms. As far as I am aware, at the time of writing this (early spring 2019), Stefanie Konig et al. (2017) provide the only assessment of accuracy when annotating a target genome (Drosphila melanogaster) with protein models from increasing distantly related taxa. Their Figure 2 is below.
A few notes are needed to fully explain this figure out of context. See also Box 4 of Yandell et al. (2012) for more information on annotation accuracy.
- these metrics are usually measured on a nucleotide basis, but they do not specify that here (alternatively they can be measured at the exon level)
- "exon F1" = "accuracy", where accuracy is the harmonic mean of specificity and sensitivity
- sensitivity = (true positives)/(true positives + false negatives); a measure of what fraction of the reference was predicted
- specificity = (true positives)/(true positives + false positives); a measure of the overlap between prediction and reference annotation
- Ding et al. (2015) find a divergence time of ~65 MY between Drosphila and Musca using mitochondial data
- divergence dates between memebers of Drosphila are given in the table below from Obbard et al. (2012), who suggest the Hawaiian Island divergence models may be less accurate than the previously published dates or those based on mutation rate
The authors of funannotate suggest avoiding proteins from closely related organisms (see the docs) because they will more often then not hurt your predictions, as can be seen by the drop off in accuracy when using more distantly related Drosophila species, above.
I'll be using the PASA pipeline wrapper in funannotate, the train command. First I'll create a new directory called PASA and then use the cleaning script provided.
(super_funannotate) [isg4@c14n06 Pamilis_funannotate]$ mkdir TRAIN
(super_funannotate) [isg4@c14n06 Pamilis_funannotate]$ cd TRAIN/When I tried to run the test command, below, I got an error saying Missing Dependencies: fasta.
(super_funannotate) [isg4@c13n01 TRAIN]$ funannotate train --input ../MASK/Final_MASK/Pamilis_FHR.clean.sort.FINALmask.fasta --out ./PASA --left ../data/1kp_transcriptome/ERR2040261_R1_.fastq --right ../data/1kp_transcriptome/ERR2040261_R2_.fastq --species "Portulaca amilis" --cpus 20 --memory 100GI checked that I had fasta installed, and I did, but the program is called fasta36 on my machine. Following this suggestion I created a symlink that solved the problem.
(super_funannotate) [isg4@c13n01 TRAIN]$ cd super_funannotate/funannotate_deps/fasta-36/bin/
(super_funannotate) [isg4@c13n05 bin]$ ln -s fasta36 fasta
(super_funannotate) [isg4@c13n05 bin]$ cd ~/scratch60/Pamilis_funannotate/TRAINHowever, when I tried to run the train command with my raw .fastq files I repeatedly received an error from BioPython: ValueError: Sequence and quality captions differ. I think this is being thrown during trimmomatic, which I've run before on these files. So I tried to use my previously generated Trinity-guided transcripts. This did not solve the problem so I took the advice of peterjc from this issue and commented out
raise ValueError("Sequence and quality captions differ.")and replaced it with
warnings.warn("Sequence and quality captions differ.", BiopythonParserWarning)which allowed funannotate to continue. After that, an error stating TransDecoder.LongOrfs and TransDecoder.Predict could not be found.
* Running CMD: super_funannotate/funannotate_deps/PASApipeline/pasa-plugins/transdecoder/TransDecoder.Predict -t __all_transcripts.fasta
Cant exec "super_funannotate/funannotate_deps/PASApipeline/pasa-plugins/transdecoder/TransDecoder.Predict": No such file or directory at super_funannotate/funannotate_deps/PASApipeline/PerlLib/Pipeliner.pm line 175.
Error, cmd: super_funannotate/funannotate_deps/PASApipeline/pasa-plugins/transdecoder/TransDecoder.Predict -t __all_transcripts.fasta died with ret -1 No such file or directory at super_funannotate/funannotate_deps/PASApipeline/PerlLib/Pipeliner.pm line 186.
Pipeliner::run(Pipeliner=HASH(0x55baaff45428)) called at super_funannotate/funannotate_deps/PASApipeline/Launch_PASA_pipeline.pl line 1044For some reason, although I had a transdecoder directory installed in $PASAHOME/pasa-plugins/, there was nothing inside it. I created a symlink to a previous TransDecoder install and that fixed the issue.
(super_funannotate) [isg4@c14n03 transdecoder]$ ln -s ~/project/apps/TransDecoder/ $PASAHOME/pasa-plugins/transdecoderFinally, the train command above ran to completion without any errors in a batch file.
#!/bin/bash
#SBATCH --job-name=TRAIN
#SBATCH --time=24:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=64G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=train-%j.out
#SBATCH --error=train-%j.err
#SBATCH --ntasks-per-node=1
funannotate train --input ../MASK/Final_MASK/Pamilis_FHR.clean.sort.FINALsoftmask.fasta --out ./Portulaca-amilis.v0-FA1.6.0/ --left ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-amilis-concat.RFT.G_1.fastq --right ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-amilis-concat.RFT.G_2.fastq --max_intronlen 10000 --species "Portulaca amilis" --cpus 20 --memory 50G --no_trimmomaticThis is again yielded an error during read normalization that I opened as issue #377 on GitHub. The issue traces back to an unrecognized flag (--min_cov 5) that funannotate is passing to Trinity's read normalization script (insilico_read_normalization.pl). I got around this by running normalization (from Trinity v2.8.5) in a standalone script
#!/bin/bash
#SBATCH --job-name=normalize
#SBATCH --time=3-00:00:00
#SBATCH -c 20
#SBATCH -p bigmem
#SBATCH --mem=512G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=normalize-%j.out
#SBATCH --error=normalize-%j.err
#SBATCH --ntasks-per-node=1
/gpfs/ysm/project/edwards/isg4/conda_envs/funannotate/opt/trinity-2.8.5/util/insilico_read_normalization.pl --PARALLEL_STATS --JM 500G --max_cov 50 --seqType fq --output ./normalize/ --CPU 20 --pairs_together --left ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-amilis-concat.RFT.G_1.fastq --right ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-amilis-concat.RFT.G_2.fastqOne quick thing beore we run train again. I originally thought that producing these normalized reads exactly where funannotate does would be the best idea because it would preserve the directory structure so that funannotate knows where to find all of the ingredients and leave output files. However, when you pass reads to funannotate it creates a symlink to the original data such that if things don't go the way you expect, or if you want to rerun analyses, the orginal files are preserved. So, when I originally placed my new left.norm.fq and right.norm.fq in training/normalize, it created a symlink to itself, which throws an error. So, to fix this mistake I moved the normalized reads to a new dirctory: training/backup/. Then I ran train again, but now with a slightly amended set of flags that skip normalization by declaring my reads are already normalized.
funannotate train --input ../MASK/Final_MASK/Pamilis_FHR.clean.sort.FINALsoftmask.fasta --out ./Portulaca-amilis.v0-FA1.6.0/ --left_norm ./Portulaca-amilis.v0-FA1.6.0/training/backup/P-amilis-concat.RFT.G_1.fastq.normalized_K25_maxC50_minC0_maxCV10000.fq --right_norm ./Portulaca-amilis.v0-FA1.6.0/training/backup/P-amilis-concat.RFT.G_2.fastq.normalized_K25_maxC50_minC0_maxCV10000.fq --max_intronlen 10000 --species "Portulaca amilis" --cpus 20 --memory 55GThe HISAT2 stats (see funannotate-train.log) show that we had an overall alignment rate of 83.82% from the normalized RNAseq reads to the genome.
31285193 reads; of these:
31285193 (100.00%) were paired; of these:
7279748 (23.27%) aligned concordantly 0 times
23872104 (76.30%) aligned concordantly exactly 1 time
133341 (0.43%) aligned concordantly >1 times
----
7279748 pairs aligned concordantly 0 times; of these:
447852 (6.15%) aligned discordantly 1 time
----
6831896 pairs aligned 0 times concordantly or discordantly; of these:
13663792 mates make up the pairs; of these:
10126049 (74.11%) aligned 0 times
3446476 (25.22%) aligned exactly 1 time
91267 (0.67%) aligned >1 times
83.82% overall alignment rateThe aligned reads were grouped into 45,522 clusters during genome-guided assembly by Trinity, and assembled into 223,935 trancripts.
The next step in the funannotate pipeline is predict. It's worth reading Jon's comments on how funannotate handles information from repeats and RNAseq data here. First I'll make the funannotate_database, which maybe should have been done earlier, but wasn't necessary. When I tried to setup the database using
(super_funannotate) [isg4@c20n11 Pamilis_funannotate]$ funannotate setup --install all --busco_db embryophyta --database funannotate_databaseThen I ran predict. I altered the weights a bit to reflect what I've found in the EVM manual and what other users have tried. It seems like most users give a lot of weight (i.e., 10) to PASA models because they are direct evidence of gene models, while ab initio predictors are weighted lightly (usually 1).
#!/bin/bash
#SBATCH --job-name=PREDICT
#SBATCH --time=2-00:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=predict-%j.out
#SBATCH --error=predict-%j.err
#SBATCH --ntasks-per-node=1
export FUNANNOTATE_DB=/gpfs/ysm/scratch60/isg4/Pamilis_FUNannotate/TRAIN/funannotate_database
funannotate predict --input ../MASK/Final_MASK/Pamilis_FHR.clean.sort.FINALsoftmask.fasta --pasa_gff ./Portulaca-amilis.v0-FA1.6.0/training/funannotate_train.pasa.gff3 --rna_bam ./Portulaca-amilis.v0-FA1.6.0/training/funannotate_train.coordSorted.bam --transcript_evidence ./Portulaca-amilis.v0-FA1.6.0/training/funannotate_train.trinity-GG.fasta --out ./Portulaca-amilis.v0-FA1.6.0/ --species "Portulaca amilis v0" --cpus 20 --genemark_mode ET --busco_seed_species embryophyta --optimize_augustus --busco_db embryophyta --organism other --weights augustus:1 hiq:1 genemark:1 pasa:10 codingquarry:0 snap:1 glimmerhmm:1 proteins:1 transcripts:1Looking at the output, we recovered ~300k models.
(super_funannotate) [isg4@c15n03 TRAIN]$ cat predict-5128113.out
Feature Specificity Sensitivity
-------------------------------------------------------
nucleotides 96.1% 92.5%
exons 84.1% 82.5%
genes 38.2% 36.1%
Feature Specificity Sensitivity
-------------------------------------------------------
nucleotides 96.3% 92.9%
exons 85.1% 83.6%
genes 42.1% 40.6%
Source Weight Count
-------------------------------------------------------
Augustus 1 38091
GeneMark 1 73646
GlimmerHMM 1 84494
HiQ 1 11410
pasa 10 27667
snap 1 68073
Total - 303381
-------------------------------------------------------
FUN_026572 Feature begins or ends in gap starting at 4468020
FUN_026574 Feature begins or ends in gap starting at 4469294
FUN_047027 Feature begins or ends in gap starting at 36649341Now I'll fix or drop the 3 gene models that were flagged for errors. To do this I downloaded the Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/predict_results/Portulaca_amilis.gff3 and looked at the flagged models in mygenomebrowser. Below FUN_025728 (model numbers vary among analyses) runs into a gap spanned with Ns, so I'll drop that model.
This was true of the other two models, so I dropped those too. In each analysis directory the list of models with errors are in predict_results/Portulaca_amilis.models-need-fixing.txt, which must contain only one gene model per line with no header or second column.
(super_funannotate) [isg4@c13n05 predict_results]$ cp Portulaca_amilis.models-need-fixing.txt Portulaca_amilis.models-to-drop.txtThen I edited this file to look like
FUN_025243
FUN_025245
FUN_044740and ran fix
(super_funannotate) [isg4@c13n05 GeneMarkES]$ funannotate fix -i ./predict_results/Portulaca_amilis.gbk -t ./predict_results/Portulaca_amilis.tbl --drop ./predict_results/Portulaca_amilis.models-to-drop.txtAfter removing the problematic models I updated the models using our RNAseq data,
#!/bin/bash
#SBATCH --job-name=UPDATE
#SBATCH --time=24:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=48G
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=update-%j.out
#SBATCH --error=update-%j.err
#SBATCH --ntasks-per-node=1
export FUNANNOTATE_DB=/gpfs/ysm/scratch60/isg4/Pamilis_FUNannotate/TRAIN/funannotate_database
funannotate update -i Portulaca-amilis.v0-FA1.6.0/ --cpus 20 --left_norm Portulaca-amilis.v0-FA1.6.0/training/normalize/left.norm.fq --right_norm Portulaca-amilis.v0-FA1.6.0/training/normalize/left.norm.fq --trinity Portulaca-amilis.v0-FA1.6.0/training/funannotate_train.trinity-GG.fasta --species "Portulaca amilis" --out Portulaca-amilis.v0-FA1.6.0/ --memory 60GAfter update finished I recoved ~60,000 transcripts that correspond to ~53,00 gene models
-------------------------------------------------------
Total Gene Models: 53,094
Total transcripts: 59,789
New Gene Models: 364
No Change: 35,328
Update UTRs: 17,344
Exons Changed: 33
Exons/CDS Changed: 25
Dropped Models: 0
CDS AED: 0.005
mRNA AED: 0.051
-------------------------------------------------------For 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
- Bao W., Kojima K.K., Kohany O. 2015. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 6:11.
- Besemer J., Borodovsky M. 2005. GeneMark: web software for gene finding in prokaryotes, eukaryotes and viruses. Nucleic Acids Res. 33:W451–4. | GeneMark
- Ding S., Li X., Wang N., Cameron S.L., Mao M., Wang Y., Xi Y., Yang D. 2015. The Phylogeny and Evolutionary Timescale of Muscoidea (Diptera: Brachycera: Calyptratae) Inferred from Mitochondrial Genomes. PLoS ONE. 10:e0134170.
- Haas B.J., Delcher A.L., Mount S.M., Wortman J.R., Smith R.K., Hannick L.I., Maiti R., Ronning C.M., Rusch D.B., Town C.D., Salzberg S.L., White O. 2003. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res. 31:5654–5666. | PASA
- Haas B.J., Salzberg S.L., Zhu W., Pertea M., Allen J.E., Orvis J., White O., Robin C.R., Wortman J.R. 2008. Automated eukaryotic gene structure annotation using EVidenceModeler and the Program to Assemble Spliced Alignments. Genome Biology. 9:R7. | Evidence Modeler
- 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.
- Huerta-Cepas J., Forslund K., Coelho L.P., Szklarczyk D., Jensen L.J., Mering von C., Bork P. 2017. Fast Genome-Wide Functional Annotation through Orthology Assignment by eggNOG-Mapper. Mol. Biol. Evol. 34:2115–2122. | eggNOG-mapper
- Jones P., Binns D., Chang H.-Y., Fraser M., Li W., McAnulla C., McWilliam H., Maslen J., Mitchell A., Nuka G., Pesseat S., Quinn A.F., Sangrador-Vegas A., Scheremetjew M., Yong S.-Y., Lopez R., Hunter S. 2014. InterProScan 5: genome-scale protein function classification. Bioinformatics. 30:1236–1240. | InterProScan 5
- Käll L., Krogh A., Sonnhammer E.L.L. 2004. A combined transmembrane topology and signal peptide prediction method. J. Mol. Biol. 338:1027–1036. | Phobius
- 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
- Krogh A., Larsson B., Heijne von G., Sonnhammer E.L. 2001. Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. J. Mol. Biol. 305:567–580. | TMHMM
- Majoros W.H., Pertea M., Salzberg S.L. 2004. TigrScan and GlimmerHMM: two open source ab initio eukaryotic gene-finders. Bioinformatics. 20:2878–2879. | GlimmerHMM
- 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.
- Petersen T.N., Brunak S., Heijne von G., Nielsen H. 2011. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nature Methods. 8:785–786. | SignalP 4.0
- Roach M.J., Schmidt S.A., Borneman A.R. 2018. Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies. BMC Bioinformatics. 19:1–10. | Purge Haplotigs
- 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
- 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
- Thomas G.W.C., Hahn M.W. 2019. Referee: Reference assembly quality scores. Genome Biol Evol. 11:1483–1486. | Referee
- 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.
# Name Version Build Channel
argtable2 2.13 0 etetoolkit
asn1crypto 0.24.0 py27_1003 conda-forge
atk 2.25.90 hb9dd440_1002 conda-forge
augustus 3.2.3 pl526hcfae127_5 bioconda
backports 1.0 py_2 conda-forge
backports.functools_lru_cache 1.5 py_1 conda-forge
backports_abc 0.5 py_1 conda-forge
bamtools 2.4.1 1 bioconda
bedtools 2.28.0 hdf88d34_0 bioconda
binutils_impl_linux-64 2.31.1 h6176602_1
binutils_linux-64 2.31.1 h6176602_6
biopython 1.73 py27h14c3975_0 conda-forge
blast 2.7.1 h4422958_6 bioconda
blat 36 0 bioconda
boost 1.67.0 py27h3e44d54_0 conda-forge
boost-cpp 1.67.0 h3a22d5f_0 conda-forge
bowtie 1.2.2 py27h2d50403_1 bioconda
bowtie2 2.3.5 py27he860b03_0 bioconda
bwa 0.7.17 h84994c4_5 bioconda
bzip2 1.0.6 h14c3975_1002 conda-forge
ca-certificates 2019.3.9 hecc5488_0 conda-forge
cairo 1.16.0 ha4e643d_1000 conda-forge
certifi 2019.3.9 py27_0 conda-forge
cffi 1.12.3 py27h8022711_0 conda-forge
chardet 3.0.4 py27_1003 conda-forge
clustalo 1.2.4 h4346872_0 etetoolkit
clustalw 2.1 h6bb024c_4 bioconda
codingquarry 2.0 py27h6bb024c_0 bioconda
collectl 4.0.4 2 bioconda
cryptography 2.6.1 py27h72c5cf5_0 conda-forge
curl 7.64.1 hf8cf82a_0 conda-forge
cycler 0.10.0 py_1 conda-forge
dbus 1.13.6 he372182_0 conda-forge
dialigntx 1.0.2 hdce4c0c_0 etetoolkit
diamond 0.9.24 ha87ae23_0 bioconda
eggnog-mapper 1.0.0 py27_0 bioconda
enum34 1.1.6 py27_1001 conda-forge
ete3 3.1.1 pyhf5214e1_0 etetoolkit
ete_toolchain 3.0.0 h73706c9_0 etetoolkit
exonerate 2.4.0 hb9dd440_2 bioconda
expat 2.2.5 hf484d3e_1002 conda-forge
fasta3 36.3.8 h470a237_4 bioconda
fastool 0.1.4 ha92aebf_3 bioconda
fasttree 2.1 hdfd2403_0 etetoolkit
fisher 0.1.4 py27h24bf2e0_1 bioconda
fontconfig 2.13.1 he4413a7_1000 conda-forge
freetype 2.10.0 he983fc9_0 conda-forge
functools32 3.2.3.2 py_3 conda-forge
futures 3.2.0 py27_1000 conda-forge
gcc_impl_linux-64 7.3.0 habb00fd_1 conda-forge
gcc_linux-64 7.3.0 h553295d_6 conda-forge
gdk-pixbuf 2.32.2 1 bioconda
gettext 0.19.8.1 hc5be6a0_1002 conda-forge
gfortran_impl_linux-64 7.3.0 hdf63c60_1
gfortran_linux-64 7.3.0 h553295d_6
giflib 5.1.9 h516909a_0 conda-forge
glib 2.58.3 hf63aee3_1001 conda-forge
gmap 2018.07.04 pl526he4cf2ce_0 bioconda
gmp 6.1.2 hf484d3e_1000 conda-forge
gnutls 3.5.19 h2a4e5f8_1 conda-forge
goatools 0.8.12 py_0 bioconda
gobject-introspection 1.56.1 py27h2da5eee_1002 conda-forge
graphite2 1.3.13 hf484d3e_1000 conda-forge
graphviz 2.40.1 h0dab3d1_0 conda-forge
gsl 2.2.1 h0c605f7_3
gst-plugins-base 1.14.4 hdf3bae2_1001 conda-forge
gstreamer 1.14.4 h66beb1c_1001 conda-forge
gtk2 2.24.31 hb68c50a_1001 conda-forge
harfbuzz 2.4.0 h37c48d4_0 conda-forge
hdf5 1.10.4 nompi_h3c11f04_1106 conda-forge
hisat2 2.1.0 py27h6bb024c_3 bioconda
hmmer 3.2.1 hf484d3e_1 bioconda
icu 58.2 hf484d3e_1000 conda-forge
idna 2.8 py27_1000 conda-forge
infernal 1.1.2 h14c3975_2 bioconda
intel-openmp 2019.3 199
ipaddress 1.0.22 py_1 conda-forge
iqtree 1.5.5 he390d98_0 etetoolkit
jellyfish 2.2.10 h6bb024c_1 bioconda
jemalloc 4.5.0 0 bioconda
jpeg 9c h14c3975_1001 conda-forge
kalign 2.03 h29c49b8_0 etetoolkit
kallisto 0.45.1 hb6a4e58_0 bioconda
kiwisolver 1.1.0 py27hc9558a2_0 conda-forge
krb5 1.16.3 h05b26f9_1001 conda-forge
libblas 3.8.0 8_mkl conda-forge
libcblas 3.8.0 8_mkl conda-forge
libcurl 7.64.1 hda55be3_0 conda-forge
libdb 6.1.26 0 bioconda
libdeflate 1.0 h14c3975_1 bioconda
libedit 3.1.20170329 hf8c457e_1001 conda-forge
libffi 3.2.1 he1b5a44_1006 conda-forge
libgcc 7.2.0 h69d50b8_2 conda-forge
libgcc-ng 8.2.0 hdf63c60_1
libgd 2.2.5 h0d07dcb_1005 conda-forge
libgfortran 3.0.0 1 conda-forge
libgfortran-ng 7.3.0 hdf63c60_0
libiconv 1.15 h516909a_1005 conda-forge
libidn11 1.34 h1cef754_0 conda-forge
liblapack 3.8.0 8_mkl conda-forge
libpng 1.6.37 hed695b0_0 conda-forge
libssh2 1.8.2 h22169c7_2 conda-forge
libstdcxx-ng 8.2.0 hdf63c60_1
libtiff 4.0.10 h648cc4a_1001 conda-forge
libtool 2.4.6 h14c3975_1002 conda-forge
libuuid 2.32.1 h14c3975_1000 conda-forge
libwebp 1.0.2 h576950b_1 conda-forge
libxcb 1.13 h14c3975_1002 conda-forge
libxml2 2.9.9 h13577e0_0 conda-forge
libxslt 1.1.32 h4785a14_1002 conda-forge
lp_solve 5.5.2.5 h14c3975_1001 conda-forge
lxml 4.3.0 py27hefd8a0e_0
mafft 7.407 1 bioconda
matplotlib 2.2.3 py27hb69df0a_0
matplotlib-base 2.2.3 py27h60b886d_1 conda-forge
metis 5.1.0 hf484d3e_1003 conda-forge
minimap2 2.16 h84994c4_1 bioconda
mkl 2019.3 199
mmtf-python 1.0.2 py27_0 bioconda
msgpack-python 0.6.1 py27h6bb024c_0 conda-forge
mummer 3.23 pl526_8 bioconda
muscle 3.8.31 he5e28f3_0 etetoolkit
mysql-connector-c 6.1.6 2 bioconda
natsort 6.0.0 py_0 conda-forge
ncurses 6.1 hf484d3e_1002 conda-forge
nettle 3.3 0 conda-forge
numpy 1.16.3 py27he5ce36f_0 conda-forge
olefile 0.46 py_0 conda-forge
openblas 0.2.20 8 conda-forge
openjdk 11.0.1 h516909a_1015 conda-forge
openssl 1.1.1b h14c3975_1 conda-forge
paml 4.8 h48adae2_0 etetoolkit
pandas 0.24.2 py27hf484d3e_0 conda-forge
pango 1.40.14 h4ea9474_1004 conda-forge
parafly r2013_01_21 1 bioconda
patsy 0.5.1 py_0 conda-forge
pcre 8.41 hf484d3e_1003 conda-forge
perl 5.26.2 h516909a_1006 conda-forge
perl-aceperl 1.92 pl526_2 bioconda
perl-algorithm-diff 1.1903 pl526_2 bioconda
perl-algorithm-munkres 0.08 pl526_1 bioconda
perl-apache-test 1.40 pl526_1 bioconda
perl-app-cpanminus 1.7044 pl526_1 bioconda
perl-appconfig 1.71 pl526_1 bioconda
perl-archive-tar 2.32 pl526_0 bioconda
perl-array-compare 3.0.1 pl526_1 bioconda
perl-autoloader 5.74 pl526_2 bioconda
perl-base 2.23 pl526_1 bioconda
perl-bio-asn1-entrezgene 1.73 pl526_0 bioconda
perl-bio-coordinate 1.007001 pl526_0 bioconda
perl-bio-featureio 1.6.905 pl526_1 bioconda
perl-bio-phylo 0.58 pl526_1 bioconda
perl-bio-samtools 1.43 pl526h1341992_1 bioconda
perl-bio-tools-phylo-paml 1.7.3 pl526_0 bioconda
perl-bio-tools-run-alignment-clustalw 1.7.4 pl526_0 bioconda
perl-bio-tools-run-alignment-tcoffee 1.7.4 pl526_1 bioconda
perl-bioperl 1.7.2 pl526_9 bioconda
perl-bioperl-core 1.7.2 pl526_3 bioconda
perl-bioperl-run 1.007002 pl526_3 bioconda
perl-business-isbn 3.004 pl526_0 bioconda
perl-business-isbn-data 20140910.003 pl526_0 bioconda
perl-cache-cache 1.08 pl526_0 bioconda
perl-capture-tiny 0.48 pl526_0 bioconda
perl-carp 1.38 pl526_3 bioconda
perl-cgi 4.40 pl526h14c3975_1 bioconda
perl-class-data-inheritable 0.08 pl526_1 bioconda
perl-class-inspector 1.34 pl526_0 bioconda
perl-class-load 0.25 pl526_0 bioconda
perl-class-load-xs 0.10 pl526h6bb024c_2 bioconda
perl-class-method-modifiers 2.12 pl526_0 bioconda
perl-clone 0.41 pl526h14c3975_1 bioconda
perl-clone-choose 0.010 pl526_0 bioconda
perl-common-sense 3.74 pl526_2 bioconda
perl-compress-raw-bzip2 2.086 pl526hf484d3e_0 bioconda
perl-compress-raw-zlib 2.086 pl526h6bb024c_0 bioconda
perl-constant 1.33 pl526_1 bioconda
perl-convert-binary-c 0.78 pl526h6bb024c_3 bioconda
perl-convert-binhex 1.125 pl526_1 bioconda
perl-cpan-meta 2.150010 pl526_0 bioconda
perl-cpan-meta-requirements 2.140 pl526_0 bioconda
perl-cpan-meta-yaml 0.018 pl526_0 bioconda
perl-crypt-rc4 2.02 pl526_1 bioconda
perl-data-dumper 2.173 pl526_0 bioconda
perl-data-optlist 0.110 pl526_2 bioconda
perl-data-stag 0.14 pl526_1 bioconda
perl-date-format 2.30 pl526_2 bioconda
perl-db-file 1.852 pl526h14c3975_0 bioconda
perl-dbd-mysql 4.046 pl526h2d50403_0 bioconda
perl-dbd-sqlite 1.60 pl526h470a237_0 bioconda
perl-dbi 1.642 pl526_0 bioconda
perl-devel-globaldestruction 0.14 pl526_0 bioconda
perl-devel-overloadinfo 0.005 pl526_0 bioconda
perl-devel-stacktrace 2.03 pl526_1 bioconda
perl-digest-hmac 1.03 pl526_3 bioconda
perl-digest-md5 2.55 pl526_0 bioconda
perl-digest-perl-md5 1.9 pl526_1 bioconda
perl-digest-sha1 2.13 pl526h6bb024c_1 bioconda
perl-dist-checkconflicts 0.11 pl526_2 bioconda
perl-dynaloader 1.25 pl526_1 bioconda
perl-email-date-format 1.005 pl526_2 bioconda
perl-encode 2.88 pl526_1 bioconda
perl-encode-locale 1.05 pl526_6 bioconda
perl-error 0.17027 pl526_1 bioconda
perl-eval-closure 0.14 pl526h6bb024c_4 bioconda
perl-exception-class 1.44 pl526_0 bioconda
perl-qer 5.72 pl526_1 bioconda
perl-exporter-tiny 1.002001 pl526_0 bioconda
perl-extutils-cbuilder 0.280230 pl526_1 bioconda
perl-extutils-makemaker 7.36 pl526_1 bioconda
perl-extutils-manifest 1.72 pl526_0 bioconda
perl-extutils-parsexs 3.35 pl526_0 bioconda
perl-file-listing 6.04 pl526_1 bioconda
perl-file-path 2.16 pl526_0 bioconda
perl-file-slurp-tiny 0.004 pl526_1 bioconda
perl-file-sort 1.01 pl526_2 bioconda
perl-file-temp 0.2304 pl526_2 bioconda
perl-file-which 1.23 pl526_0 bioconda
perl-font-afm 1.20 pl526_2 bioconda
perl-font-ttf 1.06 pl526_0 bioconda
perl-gd 2.70 pl526he941832_0 bioconda
perl-getopt-long 2.50 pl526_1 bioconda
perl-graph 0.9704 pl526_1 bioconda
perl-graphviz 2.24 pl526h734ff71_0 bioconda
perl-hash-merge 0.300 pl526_0 bioconda
perl-html-element-extended 1.18 pl526_1 bioconda
perl-html-entities-numbered 0.04 pl526_1 bioconda
perl-html-formatter 2.16 pl526_0 bioconda
perl-html-parser 3.72 pl526h6bb024c_5 bioconda
perl-html-tableextract 2.13 pl526_2 bioconda
perl-html-tagset 3.20 pl526_3 bioconda
perl-html-tidy 1.60 pl526_0 bioconda
perl-html-tree 5.07 pl526_1 bioconda
perl-html-treebuilder-xpath 0.14 pl526_1 bioconda
perl-http-cookies 6.04 pl526_0 bioconda
perl-http-daemon 6.01 pl526_1 bioconda
perl-http-date 6.02 pl526_3 bioconda
perl-http-message 6.18 pl526_0 bioconda
perl-http-negotiate 6.01 pl526_3 bioconda
perl-image-info 1.38 pl526_1 bioconda
perl-image-size 3.300 pl526_2 bioconda
perl-io-compress 2.086 pl526hf484d3e_0 bioconda
perl-io-html 1.001 pl526_2 bioconda
perl-io-sessiondata 1.03 pl526_1 bioconda
perl-io-socket-ssl 2.066 pl526_0 bioconda
perl-io-string 1.08 pl526_3 bioconda
perl-io-stringy 2.111 pl526_1 bioconda
perl-io-tty 1.12 pl526_1 bioconda
perl-io-zlib 1.10 pl526_2 bioconda
perl-ipc-cmd 1.02 pl526_0 bioconda
perl-ipc-run 20180523.0 pl526_0 bioconda
perl-ipc-sharelite 0.17 pl526h6bb024c_1 bioconda
perl-jcode 2.07 pl526_2 bioconda
perl-json 4.00 pl526_0 bioconda
perl-json-pp 4.02 pl526_0 bioconda
perl-json-xs 2.34 pl526h2d50403_2 bioconda
perl-lib 0.63 pl526_1 bioconda
perl-libwww-perl 6.36 pl526_1 bioconda
perl-libxml-perl 0.08 pl526_2 bioconda
perl-list-moreutils 0.428 pl526_1 bioconda
perl-list-moreutils-xs 0.428 pl526_0 bioconda
perl-locale-maketext-simple 0.21 pl526_2 bioconda
perl-lwp-mediatypes 6.02 pl526_3 bioconda
perl-lwp-protocol-https 6.07 pl526_4 bioconda
perl-lwp-simple 6.15 pl526h470a237_4 bioconda
perl-mailtools 2.20 pl526_0 bioconda
perl-math-cdf 0.1 pl526h14c3975_5 bioconda
perl-math-derivative 1.01 pl526_0 bioconda
perl-math-random 0.72 pl526h14c3975_2 bioconda
perl-math-spline 0.02 pl526_2 bioconda
perl-mime-base64 3.15 pl526_1 bioconda
perl-mime-lite 3.030 pl526_1 bioconda
perl-mime-tools 5.508 pl526_1 bioconda
perl-mime-types 2.17 pl526_0 bioconda
perl-mldbm 2.05 pl526_1 bioconda
perl-module-build 0.4224 pl526_3 bioconda
perl-module-corelist 5.20181218 pl526_0 bioconda
perl-module-implementation 0.09 pl526_2 bioconda
perl-module-load 0.32 pl526_1 bioconda
perl-module-load-conditional 0.68 pl526_2 bioconda
perl-module-metadata 1.000033 pl526_0 bioconda
perl-module-runtime 0.016 pl526_0 bioconda
perl-module-runtime-conflicts 0.003 pl526_0 bioconda
perl-moo 2.003004 pl526_0 bioconda
perl-moose 2.2011 pl526h2d50403_0 bioconda
perl-mozilla-ca 20180117 pl526_0 bioconda
perl-mro-compat 0.13 pl526_0 bioconda
perl-net-http 6.18 pl526_0 bioconda
perl-net-ssleay 1.86 pl526h90d6eec_0 bioconda
perl-ntlm 1.09 pl526_4 bioconda
perl-ole-storage_lite 0.19 pl526_3 bioconda
perl-package-deprecationmanager 0.17 pl526_0 bioconda
perl-package-stash 0.38 pl526hfc679d8_0 bioconda
perl-package-stash-xs 0.28 pl526hf484d3e_1 bioconda
perl-parallel-forkmanager 2.02 pl526_0 bioconda
perl-params-check 0.38 pl526_1 bioconda
perl-params-util 1.07 pl526h6bb024c_4 bioconda
perl-parent 0.236 pl526_1 bioconda
perl-parse-recdescent 1.967015 pl526_0 bioconda
perl-pathtools 3.75 pl526h14c3975_1 bioconda
perl-pdf-api2 2.033 pl526_0 bioconda
perl-perl-ostype 1.010 pl526_1 bioconda
perl-pod-escapes 1.07 pl526_1 bioconda
perl-pod-usage 1.69 pl526_1 bioconda
perl-postscript 0.06 pl526_2 bioconda
perl-role-tiny 2.000006 pl526_0 bioconda
perl-scalar-list-utils 1.45 pl526h470a237_3 bioconda
perl-scalar-util-numeric 0.40 pl526_1 bioconda
perl-set-scalar 1.29 pl526_2 bioconda
perl-soap-lite 1.19 pl526_1 bioconda
perl-socket 2.027 pl526_1 bioconda
perl-sort-naturally 1.03 pl526_2 bioconda
perl-spreadsheet-parseexcel 0.65 pl526_2 bioconda
perl-spreadsheet-writeexcel 2.40 pl526_2 bioconda
perl-statistics-descriptive 3.0702 pl526_0 bioconda
perl-storable 3.11 pl526_0 bioconda
perl-sub-exporter 0.987 pl526_2 bioconda
perl-sub-exporter-progressive 0.001013 pl526_0 bioconda
perl-sub-identify 0.12 pl526_1 bioconda
perl-sub-install 0.928 pl526_2 bioconda
perl-sub-name 0.21 pl526_1 bioconda
perl-sub-quote 2.006003 pl526_0 bioconda
perl-sub-uplevel 0.2800 pl526h14c3975_2 bioconda
perl-svg 2.84 pl526_0 bioconda
perl-svg-graph 0.02 pl526_3 bioconda
perl-task-weaken 1.06 pl526_0 bioconda
perl-template-toolkit 2.26 pl526_1 bioconda
perl-test 1.26 pl526_1 bioconda
perl-test-deep 1.128 pl526_1 bioconda
perl-test-differences 0.64 pl526_2 bioconda
perl-test-exception 0.43 pl526_2 bioconda
perl-test-harness 3.42 pl526_0 bioconda
perl-test-leaktrace 0.16 pl526h14c3975_2 bioconda
perl-test-most 0.35 pl526_0 bioconda
perl-test-requiresinternet 0.05 pl526_0 bioconda
perl-test-warn 0.36 pl526_1 bioconda
perl-text-abbrev 1.02 pl526_0 bioconda
perl-text-diff 1.45 pl526_0 bioconda
perl-text-parsewords 3.30 pl526_0 bioconda
perl-text-soundex 3.05 pl526_1000 conda-forge
perl-threaded 5.26.0 0 bioconda
perl-tie-ixhash 1.23 pl526_2 bioconda
perl-time-hires 1.9760 pl526h14c3975_1 bioconda
perl-time-local 1.28 pl526_1 bioconda
perl-timedate 2.30 pl526_1 bioconda
perl-tree-dag_node 1.31 pl526_0 bioconda
perl-try-tiny 0.30 pl526_1 bioconda
perl-type-tiny 1.004004 pl526_0 bioconda
perl-types-serialiser 1.0 pl526_2 bioconda
perl-unicode-map 0.112 pl526h6bb024c_3 bioconda
perl-uri 1.76 pl526_0 bioconda
perl-version 0.9924 pl526_0 bioconda
perl-www-robotrules 6.02 pl526_3 bioconda
perl-xml-dom 1.46 pl526_0 bioconda
perl-xml-dom-xpath 0.14 pl526_1 bioconda
perl-xml-filter-buffertext 1.01 pl526_2 bioconda
perl-xml-libxml 2.0132 pl526hbc14f71_0 bioconda
perl-xml-libxslt 1.94 pl526_1 bioconda
perl-xml-namespacesupport 1.12 pl526_0 bioconda
perl-xml-parser 2.44 pl526h3a4f0e9_6 bioconda
perl-xml-regexp 0.04 pl526_2 bioconda
perl-xml-sax 1.00 pl526_0 bioconda
perl-xml-sax-base 1.09 pl526_0 bioconda
perl-xml-sax-expat 0.51 pl526_3 bioconda
perl-xml-sax-writer 0.57 pl526_0 bioconda
perl-xml-simple 2.25 pl526_1 bioconda
perl-xml-twig 3.52 pl526_2 bioconda
perl-xml-writer 0.625 pl526_2 bioconda
perl-xml-xpath 1.44 pl526_0 bioconda
perl-xml-xpathengine 0.14 pl526_2 bioconda
perl-xsloader 0.24 pl526_0 bioconda
perl-yaml 1.27 pl526_0 bioconda
phylobayes 4.1c hac87e47_0 etetoolkit
phyml 20160115.patched hee5dff1_0 etetoolkit
pillow 6.0.0 py27he7afcd5_0 conda-forge
pip 19.1 py27_0 conda-forge
pixman 0.34.0 h14c3975_1003 conda-forge
pmodeltest 1.4 py27hb4d0d57_0 etetoolkit
psutil 5.6.2 py27h516909a_0 conda-forge
pthread-stubs 0.4 h14c3975_1001 conda-forge
pycparser 2.19 py27_1 conda-forge
pydot 1.4.1 py27_1000 conda-forge
pygraphviz 1.3.1 py27_0 bioconda
pyopenssl 19.0.0 py27_0 conda-forge
pyparsing 2.4.0 py_0 conda-forge
pyqt 5.9.2 py27h05f1152_2
pysocks 1.6.8 py27_1002 conda-forge
python 2.7.15 h721da81_1008 conda-forge
python-dateutil 2.8.0 py_0 conda-forge
python-wget 3.2 py27_0 bioconda
pytz 2019.1 py_0 conda-forge
qt 5.9.7 h52cfd70_1 conda-forge
raxml 8.2.11 h6db2ed4_0 etetoolkit
readline 7.0 hf8c457e_1001 conda-forge
recon 1.08 h470a237_1 bioconda
repeatmasker 4.0.8 pl526_14 bioconda
repeatmodeler 1.0.11 pl526_1 bioconda
repeatscout 1.0.5 h470a237_1 bioconda
reportlab 3.5.20 py27h7d98c4e_0 conda-forge
requests 2.21.0 py27_1000 conda-forge
rmblast 2.6.0 h4422958_0 bioconda
samtools 1.9 h8571acd_11 bioconda
scikit-learn 0.20.3 py27ha8026db_1 conda-forge
scipy 1.2.1 py27h09a28d5_1 conda-forge
seaborn 0.9.0 py_1 conda-forge
setuptools 41.0.1 py27_0 conda-forge
singledispatch 3.4.0.3 py27_1000 conda-forge
sip 4.19.8 py27hf484d3e_1000 conda-forge
six 1.12.0 py27_1000 conda-forge
slclust 02022010 2 bioconda
slr 1.4.3 h816ab49_0 etetoolkit
sqlite 3.26.0 h67949de_1001 conda-forge
statsmodels 0.9.0 py27h3010b51_1000 conda-forge
stringtie 1.3.6 h92e31bf_0 bioconda
subprocess32 3.5.3 py27h14c3975_0 conda-forge
suitesparse 4.5.6 heab0a99_1202 conda-forge
t_coffee 11.0.8 py27hea885bf_8 bioconda
tbb 2019.5 hc9558a2_0 conda-forge
tbl2asn 25.6 3 bioconda
tidyp 1.04 1 bioconda
tk 8.6.9 h84994c4_1001 conda-forge
tornado 5.1.1 py27h14c3975_1000 conda-forge
trf 4.09 1 bioconda
trimal 1.4.1 h6bb024c_3 bioconda
trimmomatic 0.39 1 bioconda
trinity 2.5.1 1 bioconda
trnascan-se 2.0 pl526h470a237_1 bioconda
ucsc-pslcdnafilter 357 1 bioconda
urllib3 1.24.2 py27_0 conda-forge
wheel 0.33.1 py27_0 conda-forge
xlsxwriter 1.1.7 py_0 conda-forge
xorg-kbproto 1.0.7 h14c3975_1002 conda-forge
xorg-libice 1.0.9 h516909a_1004 conda-forge
xorg-libsm 1.2.3 h84519dc_1000 conda-forge
xorg-libx11 1.6.7 h14c3975_1000 conda-forge
xorg-libxau 1.0.9 h14c3975_0 conda-forge
xorg-libxdmcp 1.1.3 h516909a_0 conda-forge
xorg-libxext 1.3.4 h516909a_0 conda-forge
xorg-libxpm 3.5.12 h14c3975_1002 conda-forge
xorg-libxrender 0.9.10 h516909a_1002 conda-forge
xorg-libxt 1.1.5 h14c3975_1002 conda-forge
xorg-renderproto 0.11.1 h14c3975_1002 conda-forge
xorg-xextproto 7.3.0 h14c3975_1002 conda-forge
xorg-xproto 7.0.31 h14c3975_1007 conda-forge
xz 5.2.4 h14c3975_1001 conda-forge
zlib 1.2.11 h14c3975_1004 conda-forge# Name Version Build Channel
_r-mutex 1.0.0 anacondar_1
bedtools 2.28.0 hdf88d34_0 bioconda
binutils_impl_linux-64 2.31.1 h6176602_1
binutils_linux-64 2.31.1 h6176602_3
bwidget 1.9.11 1
bzip2 1.0.6 h14c3975_1002 conda-forge
ca-certificates 2019.3.9 hecc5488_0 conda-forge
cairo 1.16.0 ha4e643d_1000 conda-forge
curl 7.64.1 hf8cf82a_0 conda-forge
fontconfig 2.13.1 he4413a7_1000 conda-forge
freetype 2.10.0 he983fc9_0 conda-forge
gcc_impl_linux-64 7.3.0 habb00fd_1 conda-forge
gcc_linux-64 7.3.0 h553295d_3 conda-forge
gettext 0.19.8.1 hc5be6a0_1002 conda-forge
gfortran_impl_linux-64 7.3.0 hdf63c60_1
gfortran_linux-64 7.3.0 h553295d_3
glib 2.58.3 hf63aee3_1001 conda-forge
graphite2 1.3.13 hf484d3e_1000 conda-forge
gsl 2.4 h294904e_1006 conda-forge
gxx_impl_linux-64 7.3.0 hdf63c60_1 conda-forge
gxx_linux-64 7.3.0 h553295d_3 conda-forge
harfbuzz 2.4.0 h37c48d4_0 conda-forge
icu 58.2 hf484d3e_1000 conda-forge
jpeg 9c h14c3975_1001 conda-forge
krb5 1.16.3 h05b26f9_1001 conda-forge
libblas 3.8.0 8_openblas conda-forge
libcblas 3.8.0 8_openblas conda-forge
libcurl 7.64.1 hda55be3_0 conda-forge
libdeflate 1.0 h14c3975_1 bioconda
libedit 3.1.20170329 hf8c457e_1001 conda-forge
libffi 3.2.1 he1b5a44_1006 conda-forge
libgcc-ng 8.2.0 hdf63c60_1
libgfortran-ng 7.3.0 hdf63c60_0
libiconv 1.15 h516909a_1005 conda-forge
libpng 1.6.37 hed695b0_0 conda-forge
libssh2 1.8.2 h22169c7_2 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
libxcb 1.13 h14c3975_1002 conda-forge
libxml2 2.9.9 h13577e0_0 conda-forge
make 4.2.1 h14c3975_2004 conda-forge
minimap2 2.17 h84994c4_0 bioconda
mummer4 4.0.0beta2 pl526hf484d3e_4 bioconda
ncurses 6.1 hf484d3e_1002 conda-forge
openblas 0.3.6 h6e990d7_1 conda-forge
openssl 1.1.1b h14c3975_1 conda-forge
pango 1.40.14 h4ea9474_1004 conda-forge
pcre 8.41 hf484d3e_1003 conda-forge
perl 5.26.2 h516909a_1006 conda-forge
pixman 0.34.0 h14c3975_1003 conda-forge
pthread-stubs 0.4 h14c3975_1001 conda-forge
purge_haplotigs 1.0.4 0 bioconda
r-assertthat 0.2.1 r351h6115d3f_0 conda-forge
r-base 3.5.1 h271c98b_1006 conda-forge
r-cli 1.1.0 r351h6115d3f_0 conda-forge
r-colorspace 1.4_1 r351hcdcec82_0 conda-forge
r-crayon 1.3.4 r351h6115d3f_1001 conda-forge
r-digest 0.6.18 r351h96ca727_1000 conda-forge
r-fansi 0.4.0 r351h96ca727_1000 conda-forge
r-ggplot2 3.1.1 r351h6115d3f_0 conda-forge
r-glue 1.3.1 r351hcdcec82_0 conda-forge
r-gtable 0.3.0 r351h6115d3f_0 conda-forge
r-labeling 0.3 r351h6115d3f_1001 conda-forge
r-lattice 0.20_38 r351h96ca727_1000 conda-forge
r-lazyeval 0.2.2 r351hcdcec82_0 conda-forge
r-magrittr 1.5 r351h6115d3f_1001 conda-forge
r-mass 7.3_51.4 r351hcdcec82_0 conda-forge
r-matrix 1.2_17 r351hcdcec82_0 conda-forge
r-mgcv 1.8_28 r351hcdcec82_0 conda-forge
r-munsell 0.5.0 r351h6115d3f_1001 conda-forge
r-nlme 3.1_139 r351h9bbef5b_0 conda-forge
r-pillar 1.3.1 r351h6115d3f_1000 conda-forge
r-pkgconfig 2.0.2 r351h6115d3f_1001 conda-forge
r-plyr 1.8.4 r351h29659fb_1002 conda-forge
r-r6 2.4.0 r351h6115d3f_0 conda-forge
r-rcolorbrewer 1.1_2 r351h6115d3f_1001 conda-forge
r-rcpp 1.0.1 r351h0357c0b_0 conda-forge
r-reshape2 1.4.3 r351h29659fb_1003 conda-forge
r-rlang 0.3.4 r351hcdcec82_0 conda-forge
r-scales 1.0.0 r351h29659fb_1001 conda-forge
r-stringi 1.4.3 r351h0357c0b_0 conda-forge
r-stringr 1.4.0 r351h6115d3f_0 conda-forge
r-tibble 2.1.1 r351hcdcec82_0 conda-forge
r-utf8 1.1.4 r351h96ca727_1000 conda-forge
r-viridislite 0.3.0 r351h6115d3f_1001 conda-forge
r-withr 2.1.2 r351h6115d3f_1000 conda-forge
readline 7.0 hf8c457e_1001 conda-forge
samtools 1.9 h8571acd_11 bioconda
tk 8.6.9 h84994c4_1001 conda-forge
tktable 2.10 h14c3975_0
xorg-kbproto 1.0.7 h14c3975_1002 conda-forge
xorg-libice 1.0.9 h516909a_1004 conda-forge
xorg-libsm 1.2.3 h84519dc_1000 conda-forge
xorg-libx11 1.6.7 h14c3975_1000 conda-forge
xorg-libxau 1.0.9 h14c3975_0 conda-forge
xorg-libxdmcp 1.1.3 h516909a_0 conda-forge
xorg-libxext 1.3.4 h516909a_0 conda-forge
xorg-libxrender 0.9.10 h516909a_1002 conda-forge
xorg-renderproto 0.11.1 h14c3975_1002 conda-forge
xorg-xextproto 7.3.0 h14c3975_1002 conda-forge
xorg-xproto 7.0.31 h14c3975_1007 conda-forge
xz 5.2.4 h14c3975_1001 conda-forge
zlib 1.2.11 h14c3975_1004 conda-forge-------------------------------------------------------
Checking dependencies for funannotate v1.6.0-dfd805f
-------------------------------------------------------
You are running Python v 2.7.15. Now checking python packages...
biopython: 1.73
goatools: 0.8.12
matplotlib: 2.2.3
natsort: 6.0.0
numpy: 1.16.3
pandas: 0.24.2
psutil: 5.6.2
requests: 2.21.0
scikit-learn: 0.20.3
scipy: 1.2.1
seaborn: 0.9.0
All 11 python packages installed
You are running Perl v 5.026002. Now checking perl modules...
Bio::Perl: 1.007002
Carp: 1.38
Clone: 0.41
DBD::SQLite: 1.60
DBD::mysql: 4.046
DBI: 1.642
DB_File: 1.852
Data::Dumper: 2.173
File::Basename: 2.85
File::Which: 1.23
Getopt::Long: 2.5
Hash::Merge: 0.300
JSON: 4.00
LWP::UserAgent: 6.36
Logger::Simple: 2.0
POSIX: 1.76
Parallel::ForkManager: 2.02
Pod::Usage: 1.69
Scalar::Util::Numeric: 0.40
Storable: 3.11
Text::Soundex: 3.05
Thread::Queue: 3.13
Tie::File: 1.02
URI::Escape: 3.31
YAML: 1.27
threads: 2.21
threads::shared: 1.59
All 27 Perl modules installed
Checking external dependencies...
CodingQuarry: 2.0
RepeatMasker: RepeatMasker 4.0.8
RepeatModeler: RepeatModeler version DEV
Trinity: 2.5.1
augustus: 3.2.3
bamtools: bamtools 2.4.1
bedtools: bedtools v2.28.0
blat: BLAT v36
diamond: diamond 0.9.24
emapper.py: emapper-0.12.7
ete3: 3.1.1
exonerate: exonerate 2.4.0
fasta: no way to determine
gmap: 2018-07-04
gmes_petap.pl: 4.38
hisat2: 2.1.0
hmmscan: HMMER 3.2.1 (June 2018)
hmmsearch: HMMER 3.2.1 (June 2018)
java: 11.0.1
kallisto: 0.45.1
mafft: v7.407 (2018/Jul/23)
makeblastdb: makeblastdb 2.6.0+
minimap2: 2.16-r922
nucmer: 3.1
pslCDnaFilter: no way to determine
rmblastn: rmblastn 2.6.0+
samtools: samtools 1.9
stringtie: 1.3.6
tRNAscan-SE: 2.0 (December 2017)
tbl2asn: unknown, likely 25.3
tblastn: tblastn 2.6.0+
trimal: trimAl v1.4.rev15 build[2013-12-17]
All 32 external dependencies are installed
Checking Environmental Variables...
$FUNANNOTATE_DB=/gpfs/ysm/scratch60/isg4/Pamilis_funannotate/TRAIN/funannotate_database
$PASAHOME=super_funannotate/funannotate_deps/PASApipeline
$TRINITYHOME=super_funannotate/opt/trinity-2.5.1
$EVM_HOME=super_funannotate/funannotate_deps/evidencemodeler
$AUGUSTUS_CONFIG_PATH=super_funannotate/config/
$GENEMARK_PATH=super_funannotate/funannotate_deps/gmes_petap
$BAMTOOLS_PATH=super_funannotate/bin
All 7 environmental variables are set
-------------------------------------------------------