Notes

Kaiju Databases

  • MetaFunc relies on Kaiju’s nr_euk database

  • The nr_euk database is constructed from NCBI’s nr.gz. An entry line in nr.gz contains several accession numbers of proteins that are redundant (i.e. identical sequences).

  • As part of the database build, a .faa-file is made wherein ONLY the first protein of a fasta line from nr.gz is retained along with the sequences. A taxonomy ID of the protein is also included as a suffix after the accession number (e.g. WP_003131952.1_1357). If the proteins of the nr.gz come from more than one species or strain, the taxonomy ID of the lowest common ancestor of these organisms is appended as the suffix.

    • Example line in the .faa-file may look like:

    > WP_003131952.1_1357
MAQQRRGGFKRRKKVDFIAANKIEVVDYKDTELLKRFISERGKILPRRVTGTSAKNQRKVVNAIKRARVMALLPFVAEDQN

    where WP_003131952.1 is the first protein in a fasta line in nr.gz and 1357 is the lowest common ancestors of all redundant proteins in that fasta line.

Creating the NRGO Databases

Because only the first protein is in the output of Kaiju, the nrgo-pipeline first collects all redundant protein entries and builds a set of unique GO-terms for it.

However, GO annotations may not be suitable across different species, e.g. if proteins in a nr.gz line come from different species. Thus, in creating the NR –> GO database, only proteins identified to the species-level are considered.

Steps in Database Creation:

  1. Only species-level entries in the .faa files are obtained

  2. Only nr.gz lines whose first proteins match those in (1) are obtained. Note that the Kaiju output only is based on these proteins. There might be instances where the first protein is not found in the NCBI’s prot.accession2taxid file. In these lines the first accession with an entry in the prot.accession2taxid is used in Kaiju. For such lines in nr.gz, we compare the proteins in the entire line to (1) and use the protein that matches as the first protein.

  3. All proteins in a nr.gz line are converted to their uniprot counterparts using the idmapping_selected.tab.gz file from UniProt. Only Refseq, EMBL-CDS, and UniProtKB-AC accession numbers are considered. The corresponding GOs of the UniProt counterparts are obtained using the file goa_uniprot_all.gaf.gz from EBI.

  4. All GOs from (3) are then associated to the first protein in an SQLite database.

  5. The script prot2go.py in the pipeline takes the output protein accessions from Kaiju and looks up their GO annotations from the database in (4).

Description of the Process to Associate Function

The script that gathers gene ontology (GO) annotations for protein accession numbers does the following:

For each table in Directory: source/protein/, this script will get all protein accessions from reads that have been classified to the species level (and passed abundance cutoff) and that belong to a taxa that the user specified. “all” taxon is also an output but take note that proteins of organisms from different kingdoms could be annotated with the same GOs. For this pipeline, we use proportional read counts of the protein accessions (and will be referred to as ’reads’ henceforth).

Please see the following code for options:

python prot2go.py -h

If ‘grouping by condition’ is wanted, reads from samples that belong to a group/condition are averaged per accession number. The default of –-nogroup treats each sample separately. Each protein accession number is looked up in a nr(protein) –> GO sqlite database (see Databases for description of this database) to get a list of GO annotations for the specific protein accession number. We update this list of GO annotations by getting the ancestors that are related to these GO numbers by ‘is_a’ or ‘part_of’ terms (note that the entire list of GOs are considered in the update such that the GO terms and path/s to the top of the GO DAG are not doubled for each accession number). Each GO annotation of the accession number, including updated accessions, is then assigned its corresponding read and scaled read counts (summing up each time the same GO annotates another protein accession number). GOs are separated by their namespaces (Biological Process, Molecular Function, and Cellular Component) and percentage of reads covering a gene ontology term is calculated by dividing the scaled read count of that term by the total scaled read counts of the namespace and multiplying by 100% (for species-scaled, the same is done using species-scaled values). The percentage is obtained per namespace but the final output table contains all namespaces together, thus the percent columns would add up to more than 100%. If there are no GO annotating a protein accession number, it is binned into a ‘None’ file, which lists the protein accession numbers that do not have a GO annotation.

Note

The python script prot2go.py has –-all_ranks as an option allowing users to use proteins that matched to reads with taxon IDs at any taxon level. The default uses only proteins whose matches to reads identify taxa at species level. As such this option is still at beta stage and anyone using –-all_ranks should do so with caution.