dammit!
dammit is an annotation
pipeline written by Camille
Scott. dammit runs a relatively standard annotation
protocol for transcriptomes: it begins by building gene models with Transdecoder,
and then
uses the following protein databases as evidence for annotation:
Pfam-A, Rfam,
OrthoDB,
uniref90 (uniref is optional with
--full).
If a protein dataset is available, this can also be supplied to the
dammit pipeline with --user-databases as optional evidence for
annotation.
In addition, BUSCO v3 is run, which will compare the gene content in your transcriptome
with a lineage-specific data set. The output is a proportion of your
transcriptome that matches with the data set, which can be used as an
estimate of the completeness of your transcriptome based on evolutionary
expectation (Simho et al.
2015).
There are several lineage-specific datasets available from the authors
of BUSCO. We will use the metazoa dataset for this transcriptome.
Annotation necessarily requires a lot of software! dammit attempts to simplify this and make it as reliable as possible, but we still have some dependencies.
sudo apt-get -y install python3-dev hmmer unzip \
infernal ncbi-blast+ liburi-escape-xs-perl emboss liburi-perl \
build-essential libsm6 libxrender1 libfontconfig1 \
parallel libx11-dev python3-venv last-align transdecoder
Create a python 3 environment for dammit:
python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
Install shmlast
pip install -r <(curl https://raw.githubusercontent.com/camillescott/shmlast/master/requirements.txt)
pip install --upgrade pip
pip install shmlast
Install the proper version of GNU parallel:
cd
(wget -O - pi.dk/3 || curl pi.dk/3/ || fetch -o - http://pi.dk/3) | bash
sudo cp $HOME/bin/parallel /usr/bin/parallel
and then BUSCO...
cd
git clone https://gitlab.com/ezlab/busco.git
pushd busco && python setup.py install && popd
export PATH=$HOME/busco/scripts:$PATH
echo 'export PATH=$HOME/busco/scripts:$PATH' >> ~/py3/bin/activate
Finally, install dammit from the refactor/1.0 branch:
pip install https://github.com/camillescott/dammit/archive/refactor/1.0.zip
dammit has two major subcommands: dammit databases and dammit annotate. databases
checks that the databases are installed and prepared, and if run with the --install flag,
will perform that installation and preparation. If you just run dammit databases on its
own, you should get a notification that some database tasks are not up-to-date -- we need
to install them!
Unless we're running short on time, we're going to do a full run. If you want to run a quick
version of the pipeline, add a parameter, --quick, to omit OrthoDB, Uniref, Pfam, and Rfam.
A "full" run will take longer to install and run, but you'll have access to the full annotation pipeline.
dammit databases --install --busco-group metazoa # --quick
We used the "metazoa" BUSCO group. We can use any of the BUSCO databases, so long as we install
them with the dammit databases subcommand. You can see the whole list by running
dammit databases -h. You should try to match your species as closely as possible for the best
results. If we want to install another, for example:
dammit databases --install --busco-group fungi # --quick
Note: if you have limited space on your instance, you can also install these databases in a different location (e.g. on an external volume). You would want to run this command before running the database installs we just ran.
#Run ONLY if you want to install databases in different location.
#To run, remove the `#` from the front of the following command:
# dammit databases --database-dir /path/to/databases
Keep things organized! Let's make a project directory:
export PROJECT=/mnt/work
cd $PROJECT
mkdir -p annotation
cd annotation
You all ran Trinity earlier to generate an assembly, but just in case, we're going to download a version of that assembly to annotate.
curl -OL https://raw.githubusercontent.com/ngs-docs/angus/2017/_static/Trinity.fasta
mv Trinity.fasta trinity.nema.fasta
Now we'll download a custom Nematostella vectensis protein database available from JGI. Here, somebody has already created a proper database for us [1] (it has a reference proteome available through uniprot). If your critter is a non-model organism, you will likely need to create your own with proteins from closely-related species. This will rely on your knowledge of your system!
curl -LO ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Eukaryota/UP000001593_45351.fasta.gz
gunzip -c UP000001593_45351.fasta.gz > nema.reference.prot.faa
Run the command:
dammit annotate trinity.nema.fasta --busco-group metazoa --user-databases nema.reference.prot.faa --n_threads 6 # --quick
While dammit runs, it will print out which tasks its running to the terminal. dammit is written with a library called pydoit, which is a python workflow library similar to GNU Make. This not only helps organize the underlying workflow, but also means that if we interrupt it, it will properly resume!
After a successful run, you'll have a new directory called trinity.nema.fasta.dammit. If you
look inside, you'll see a lot of files:
ls trinity.nema.fasta.dammit/
annotate.doit.db trinity.nema.fasta.dammit.namemap.csv trinity.nema.fasta.transdecoder.pep
dammit.log trinity.nema.fasta.dammit.stats.json trinity.nema.fasta.x.nema.reference.prot.faa.crbl.csv
run_trinity.nema.fasta.metazoa.busco.results trinity.nema.fasta.transdecoder.bed trinity.nema.fasta.x.nema.reference.prot.faa.crbl.gff3
tmp trinity.nema.fasta.transdecoder.cds trinity.nema.fasta.x.nema.reference.prot.faa.crbl.model.csv
trinity.nema.fasta trinity.nema.fasta.transdecoder_dir trinity.nema.fasta.x.nema.reference.prot.faa.crbl.model.plot.pdf
trinity.nema.fasta.dammit.fasta trinity.nema.fasta.transdecoder.gff3
trinity.nema.fasta.dammit.gff3 trinity.nema.fasta.transdecoder.mRNA
The most important files for you are trinity.nema.fasta.dammit.fasta,
trinity.nema.fasta.dammit.gff3, and trinity.nema.fasta.dammit.stats.json.
If the above dammit command is run again, there will be a message:
**Pipeline is already completed!**
Cammille wrote dammit in Python, which includes a library to parse gff3 dammit output. To send this output to a useful table, we will need to open the Python environemnt.
cd trinity.nema.fasta.dammit
python
Then, manually enter each line of code to output a list of gene ID:
import pandas as pd
from dammit.fileio.gff3 import GFF3Parser
gff_file = "trinity.nema.fasta.dammit.gff3"
annotations = GFF3Parser(filename=gff_file).read()
names = annotations.sort_values(by=['seqid', 'score'], ascending=True).query('score < 1e-05').drop_duplicates(subset='seqid')[['seqid', 'Name']]
new_file = names.dropna(axis=0,how='all')
new_file.head()
new_file.to_csv("nema_gene_name_id.csv")
exit()
This will output a table of genes with 'seqid' and 'Name' in a .csv file: nema_gene_name_id.csv. Let's take a look at that file:
less nema_gene_name_id.csv
Notice there are multiple transcripts per gene model prediction. This .csv file can be used in tximport in downstream DE analysis.
- Putnam NH, Srivastava M, Hellsten U, Dirks B, Chapman J, Salamov A, Terry A, Shapiro H, Lindquist E, Kapitonov VV, Jurka J, Genikhovich G, Grigoriev IV, Lucas SM, Steele RE, Finnerty JR, Technau U, Martindale MQ, Rokhsar DS. (2007) Sea anemone genome reveals ancestral eumetazoan gene repertoire and genomic organization. Science. 317, 86-94.