gaps/SETUP

Setting up GapMind:

These scripts should work on any Linux system, and would probably work
on other Unix or MacOS as well. All of the code is written in perl
5. (I am using perl v5.10.1 or perl v5.16.3.) The other major
requirements are HMMer 3, either usearch or diamond, and sqlite3.

There are three phases for running GapMind: 

1. Building a "curated" database of proteins of known function
2. Converting the pathway descriptions to a steps database
3. Using the curated and steps databases to analyze genomes.

If you just want to use the existing definitions of pathways and steps
to analyze your genomes, then you only need to run the 3rd phase, so
setting up phase 3 is described first.

ANALYZE GENOMES (PHASE 3)

Make sure you have the following perl modules installed:
  DBI FindBin Getopt::Long List::Util

If you want to use the web server CGI (gapView.cgi) to view the
results, you will also need these perl modules:
  CGI HTML::Entities IO::Handle DB_File File::stat Time::HiRes
  LWP::Simple JSON IO::String Bio::SeqIO
  (Bio::SeqIO is part of bioperl; DB_File depends on the BerkeleyDB
   software library)

Make sure sqlite3 is on your path and supported by DBI (this may
require installing the DBD::SQLite perl library). I am currently using
sqlite version 3.7.17.

Put the usearch executable in the bin/ subdirectory. I am currently
using usearch 11.0.667. The 64-bit executables for usearch are now
freely available, see
https://github.com/rcedgar/usearch_old_binaries

Or, if using diamond, put the diamond executable in the bin/
subdirectory. I have tested diamond v2.1.6.160. It gives very similar
results as usearch.

Put the hmmsearch and hmmfetch executables in the bin/
subdirectory. I am currently using HMMER 3.3.1. Source code for that
version is available at
http://eddylab.org/software/hmmer/hmmer-3.3.1.tar.gz

Put the BLAST executables (formatdb, blastall, and fastacmd)
in the bin/blast/ subdirectory. I am currently using NCBI BLAST
2.2.18 (not BLAST+), which is available at
ftp://ftp.ncbi.nlm.nih.gov/blast/executables/legacy.NOTSUPPORTED/2.2.18/
(BLAST is used by gapView.cgi, but is not needed to run GapMind analyses.)

Make sure the steps files are up to date (i.e. from the gaps/aa/
or gaps/carbon/ subdirectory in github).

Create some subdirectories:

mkdir fbrowse_data
mkdir private
mkdir tmp
mkdir tmp/downloaded

Download the curated and steps databases, i.e. from
  https://papers.genomics.lbl.gov/tmp/path.aa/curated.faa
  https://papers.genomics.lbl.gov/tmp/path.aa/curated.db
  https://papers.genomics.lbl.gov/tmp/path.aa/steps.db
for amino acid biosynthesis, or from
  https://papers.genomics.lbl.gov/tmp/path.carbon/curated.faa
  https://papers.genomics.lbl.gov/tmp/path.carbon/curated.db
  https://papers.genomics.lbl.gov/tmp/path.carbon/steps.db
for carbon catabolism, and put them into
tmp/path.aa or tmp/path.carbon
(You will need to create the tmp/path.aa and tmp/path.carbon
directories first.)

Extract the HMM files from the steps database, i.e.
# for amino acid biosynthesis
bin/extractHmms.pl tmp/path.aa/steps.db tmp/path.aa
# for cabon
bin/extractHmms.pl tmp/path.carbon/steps.db tmp/path.carbon

Format the database of characterized proteins for usearch (or diamond) and BLAST,
i.e. for amino acid biosynthesis do:

bin/blast/formatdb -p T -o T -i tmp/path.aa/curated.faa
  and then
bin/usearch -makeudb_ublast tmp/path.aa/curated.faa -output tmp/path.aa/curated.faa.udb
  or
bin/diamond makedb --in tmp/path.aa/curated.faa -d tmp/path.aa/curated.faa.dmnd

(The BLAST database is needed for viewing results on the web site, but
not for running the GapMind analysis itself.)

Streamlined instructions for all the setup steps above, including
getting some of the perl libraries, are described here.

  https://github.com/morgannprice/PaperBLAST/issues/16#issuecomment-2141157340

To run the analysis using the CGI scripts and your web server
(this requires usearch):

Make sure that the tmp and tmp/downloaded directories are writable by
apache (or by everyone). Then visit cgi/gapView.cgi to run the
analysis and view the results.

To run the analysis using the CGI scripts from the command line
(this requires usearch):

The key script is gapView.cgi. The key arguments are the set (aa or
carbon), gdb (which genome database), and gid (which genome
identifier). For example, to analyze assembly GCF_902167245.1 of
Desulfovibrio vulgaris Hildenborough from NCBI, do

cd cgi
./gapView.cgi 'force=1&set=aa&gdb=NCBI&gid=GCF_902167245.1'

You can use the CPU_USE environment variable to control how many
threads gapView.cgi uses.

The results will be in tmp/NCBI__GCF_902167245.1/aa.sum.* --
aa.sum.rules, aa.sum.steps, and aa.sum.cand are tab-delimited tables
with one row per rule, per step, or per candidate for a step, and
aa.sum.db is a sqlite3 database. The first thing I usually look at is
the results for rule=all for each pathway. nHi, nMed, and nLo show the
total number of high-confidence, medium-confidence, and low-confidence
steps. In particular, if nMed = 0 and nLo = 0 that means that every
step in the pathway is high-confidence, and the entire pathway can be
considered high confidence. If there are no low-confidence steps
(nLo=0), or if the only low-confidence steps are known gaps, then the
pathway is probably present. expandedPath has the list of steps for
the highest-confidence path.

aa.sum.steps has one line per step; look there to see what the actual
best candidates are. For more details about the scoring of the
candidates, see aa.sum.cand, which has one row per candidate for each
step.

If you're interested in known gaps, then aa.sum.knownsim will report
which organisms in your data set ("orgId") are similar to organisms
with known gaps ("orgId2"). For the actual known gaps see
gaps/aa/aa,known.steps.gaps and gaps/aa/aa.known.gaps.orgs (both
tab-delimited).

To analyze your own genomes, use bin/buildorgs.pl to convert the
genome(s) of interest to a table of organisms and their proteins, in
the format expected by gapView.cgi. For example:

# Download a protein fasta file
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/167/245/GCF_902167245.1_DSV_KAUST_v1/GCF_902167245.1_DSV_KAUST_v1_protein.faa.gz
gunzip GCF_902167245.1_DSV_KAUST_v1_protein.faa.gz

# Convert it to the form GapMind is set up, with orgs=DvH, so
# that the results will be in tmp/DvH/
mkdir tmp/DvH
bin/buildorgs.pl -out tmp/DvH/orgs -orgs 'file:GCF_902167245.1_DSV_KAUST_v1_protein.faa:Desulfovibrio vulgaris downloaded'

# Run gapView with the orgs argument instead of the gdb/gid arguments
(cd cgi; ./gapView.cgi 'force=1&set=aa&orgs=DvH')

You can also use buildorgs.pl to set up many genomes at once (all in
one directory). The out.org file will have a resulting table of
genomes and the out.faa file will have the combined fasta protein
sequences. If you want to analyze thousands of genomes, or if you're
using a 32-bit version of usearch, you should split them into groups
and then use bin/combineOrgs.pl and bin/combineGaps.pl to combine the
results. (Or use diamond with the step-by-step instructions below.)

To run the analysis from the command line, step by step:

# First set up tmp/DvH/orgs as described above, then

# Use bin/gapsearch.pl to compare the queries (one file per pathway) to
# the proteins.
bin/gapsearch.pl -orgs tmp/DvH/orgs -set aa -out tmp/DvH/aa.hits -nCPU 4

# or if using diamond:
diamond makedb --in tmp/DvH/orgs.faa -d tmp/DvH/orgs.aa.dmnd
bin/gapsearch.pl -diamond -orgs tmp/DvH/orgs -set aa -out tmp/DvH/aa.hits -nCPU 4

# Use bin/gaprevsearch.pl to see if the candidates are similar to
# characterized proteins that have other functions.
bin/gaprevsearch.pl -orgs tmp/DvH/orgs -hits tmp/DvH/aa.hits -curated tmp/path.aa/curated.faa.udb -out tmp/DvH/aa.revhits -nCPU 4

#  or if using diamond:
bin/gaprevsearch.pl -diamond -orgs tmp/DvH/orgs -hits tmp/DvH/aa.hits -curated tmp/path.aa/curated.faa.dmnd -out tmp/DvH/aa.revhits -nCPU 4

# Use bin/gapsummary.pl to score all the candidates and steps and to
# find the best-scoring path for each pathway. This produces three
# tables: one for each pathway and rule (including the rule "all" for
# the entire pathway); one for each step; and one for each candidate for
# each step.
bin/gapsummary.pl -orgs tmp/DvH/orgs -set aa -hits tmp/DvH/aa.hits -rev tmp/DvH/aa.revhits -out tmp/DvH/aa.sum

# Optionally, use bin/checkGapRequirements.pl to check dependencies
# between pathways. The output table will list any warnings.
bin/checkGapRequirements.pl -org DvH -set aa -out tmp/DvH/aa.sum.warn

# Optionally, use orgsVsMarkers.pl to compare the genome to organisms
# with known gaps (for amino acid biosynthesis only).
bin/orgsVsMarkers.pl -orgs tmp/DvH/orgs -vs gaps/aa/aa.known.gaps.markers.faa -out tmp/DvH/aa.sum.knownsim

# Optionally, use buildGapsDb.pl to combine the results into a sqlite3
# database.
bin/buildGapsDb.pl -gaps tmp/DvH/aa.sum -requirements tmp/DvH/aa.sum.warn -steps tmp/path.aa/steps.db -out tmp/DvH/aa.sum.db -markersim tmp/DvH/aa.sum.knownsim

CONVERT PATHWAY DESCRIPTIONS TO A STEPS DATABASE (PHASE 2)

This requires the curated.db and curated.faa files to be in
tmp/path.setName, where setName is aa or carbon.

It also requires two databases of HMMs: PFam and TIGRFam. Put these in
the hmm/ directory (relative to the code base) and set it up to
include Pfam-A.hmm and TIGRFAMs.hmm (in our case, a symbolic link to
TIGRFAMs_15.0_HMM.LIB).  You'll need to run hmmpress on both of these
files. Finally, you'll need tab-delimited tables with metadata about
the HMMs. For TIGRFam, which is no longer updated, you can use the
table in the PaperBLAST code base:

cp gaps/tigrinfo hmm/

For PFam, use this command to make the table:

perl -e 'print "acc\tname\n"; while(<STDIN>) { chomp; $name = $1 if m/^NAME\s+(.*)$/; $acc = $1 if m/^ACC\s+(.*)$/; next unless $_ eq "//"; print "$acc\t$name\n" if $name && $acc; $name = undef; $acc = undef; }' < hmm/Pfam-A.hmm > hmm/pfam.tab

Once you have set up curated.faa*, curated2.faa, pfam, and tigrfam, you
can compute the queries and build the steps database, i.e.

bin/buildStepsDb.pl -set carbon -doquery

You should now be ready to run gapsearch.pl or gapView.cgi.

The optional -doquery option will run gapquery.pl for each
pathway. Or, run gapquery.pl on each pathway before running
buildStepsDb without the -doquery option. Here is an example of
running gapquery.pl for pathway phenylacetate in set = carbon:

  bin/gapquery.pl -hmm hmm -steps gaps/carbon/phenylacetate.steps -out tmp/path.carbon

BUILD THE CURATED DATABASE (PHASE 1)

This requires two protein databases: a list of characterized and curated
proteins (from the characterized subset of PaperBLAST), and a
supplementary set of curated mostly-uncharacterized proteins (from
Swiss-Prot).

If you have the PaperBLAST database (usually in data/) and the
ind/uniprot_sprot.dat.gz file from downloading Swiss-Prot, then you
can set up the various databases that are needed by running:

zcat ind/uniprot_sprot.dat.gz | bin/sprotCharacterized.pl > sprot.curated_parsed
bin/setupGaps.pl  -ind ind -set setName -data data -sprot sprot.curated_parsed

where setName would be aa for the amino acid biosynthesis pathways

DEFINE YOUR OWN SET OF PATHWAYS

Choose a new setName and create a new gaps/setName directory in the
source code base.

For each pathway, create a pathway.steps file. There is some
documentation of the file format in bin/gapquery.pl, or see examples
in gaps/aa/*.steps.

Create a tab-delimited file setName/setName.table, with the fields
pathwayId and desc. It should include a row for each pathwayId that
has a pathwayId.steps file and for pathwayId = all. For an example,
see gaps/aa/aa.table

Optional: create a tab-delimited file of curated gaps in
gaps/setName/setName.curated.gaps.tsv. It must have a header line the
fields gdb gid pathway step class comment. (gdb is the genome
database, such as "NCBI", and gid is the genome identifier.) There
need not be any other rows. For an example, see
gaps/aa/aa.curated.gaps.tsv

Optional: if you will be using gapView.cgi to view the results, or
if you want to note dependencies between pathways, create a
tab-delimited file of dependencies. See gaps/aa/requires.tsv for an
example. The required fields are rule, requires, and
comment. Dependencies are optional but the header line is required.

Once these files are created, you can create a steps database as
described above using buildStepsDb.pl.