---
title: "Genome annotation notebook"
author: "Nina Dombrowski"
affiliation: "NIOZ"
date: "`r Sys.Date()`"
knit: (function(input_file, encoding) {out_dir <- 'docs';rmarkdown::render(input_file,encoding=encoding,output_file=file.path(dirname(input_file), out_dir, 'Annotations_tutorial.html'))})
output:
  rmdformats::readthedown:
    highlight: kate
    code_folding: hide
editor_options: 
  
  chunk_output_type: console
---

```{r, echo=FALSE}
library(htmltools)
htmltools::img(src = knitr::image_uri("/Users/ninadombrowski/Desktop/WorkingDir/General/Pictures_general/text4236.png"),                alt = 'logo', 
               style = 'position:absolute; top:50px; right:100px; padding:10px;height:100px;width:100px')

```



```{css, echo=FALSE}
.spoiler {
  visibility: hidden;
}

.spoiler::before {
  visibility: visible;
  content: "Spoiler alert! Hover over me to see the answer.";
  font-weight: bold
}

.spoiler:hover {
  visibility: visible;
}

.spoiler:hover::before {
  display: none;
}
```


```{r knitr setup, include=FALSE,  eval=TRUE, echo=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
options(knitr.table.format = "html")
knitr::opts_chunk$set(eval=TRUE, cache=FALSE, message=FALSE, warning=FALSE, 
                      comment = "", results="markup")
#https://bookdown.org/yihui/rmarkdown/html-document.html
#install.packages('knitr', ependencies = TRUE)
#install.packages("devtools", lib="~/R/lib")
#library(DT)
#devtools::session_info()
```


###################################################################################
###################################################################################
# Introduction
###################################################################################
###################################################################################

The goal of this tutorial is to annotate the proteins from two genomes. Therefore, we will:

- download two genomes from NCBI  
- annotate them using several databases such as the COGs, arCOGS, KOs, PFAM, ...
- combine the results into one comprehensive table.

If you already have proteins for your genome of interest, you do not need to do the first steps. You can either start with the mapping file part (which runs prokka for gene calling and includes some bash steps to add summary statistics, such as the bin id, contig name, contig length to the protein files) or if you already have protein files you can start at the step 'do the annotations'. The only thing you might want to generate is a file that lists all the proteins you are working with.

**General comments**

- A lot of programs in this workflow use for loops or GNU parallel to speed things up. While this might not make so much sense for only two genomes, if you work with thousands of files this becomes very useful. 
- For most annotation steps we can use of several CPUs. For this tutorial we only use few CPUs but this should be adjusted according to the available resources you have available. More on how to check resources can be found in the [Unix tutorial](https://github.com/ndombrowski/Unix_tutorial)
- This script uses a lot of AWK for parsing. If you want to have some more insights into how to use AWK we also have a [AWK tutorial](https://github.com/ndombrowski/AWK_tutorial) available.
- Once you have the annotation table, we provide an [R tutorial here](https://github.com/ndombrowski/R_tutorial) that gives an example how to summarize the data.
- The individual database searches might take a bit (esp. the interproscan and the diamond search against the NCBI database). 
- We will provide all intermediate data for you in a repository so you can have a look at what kind of data is generated.

**General info on this tutorial**

- `r text_spec("grey box.", background = "#BEBEBE", color = "black", bold = F)` : The code that we want to execute
- `r text_spec("red box.", background = "#D05A6E", color = "white", bold = T)` : If text is highlighted by a red box,  this highlights additional code that is useful but code that we WILL NOT use during this tutorial
- `r text_spec("Hidden code", background = "white", color = "black", bold = F)` : In some cases the code might be hidden and you only see a little box at the right hand side. Hidden codes means that from the previous sessions you should have enough knowledge to do this step or it is is a test of your skills. If you run into issues you can open the code and check how to run your command by clicking on the `Code` button to reveal the code
- `r text_spec("hyperlink", background = "white", color = "red", bold = F)` : If the text is highlighted like this during a text section, then you can click on it and get forwarded to websites with more information, i.e. a publication or an online tutorial.
- **Exercises/Questions/Homework:** : This will appear if there are some questions for you do check if you understood everything but this is not required to finish the tutorial. Some of these are more time-consuming and might be best to run on your own after you have finished the complete workflow.
- Sometimes we also just asked questions to think about the output of certain programs. The text is then hidden and you have to hover over the **"spoiler"** and or click a **Show/hide button**to see the answer to the question.




###################################################################################
###################################################################################
# General
###################################################################################
###################################################################################

Below you find a list of the used programs. If you want to run this pipeline, please install the listed Custom scripts and databases can be found in a [Zenodo repository](https://zenodo.org/deposit/6489670).

**General notice for people interested in eukaryotic genomes:**

This tutorial is designed specifically with prokaryotic genomes in mind but can be used for for eukaryotic proteins as well. This especially is easy if you already have proteins for your genomes of interest. If you do not have such files beware that the tool we use for protein-calling, prokka, is specific for prokaryotic genomes. Tools for identifying eukaryotic proteins are, among others, [Augustus](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160219/) for single genomes or [MetaEuk](https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-020-00808-x) for metagenomes.

Something to keep in mind with the used databases is that these were often generated using only archaeal and/or bacterial genomes. Only of a few of the databases used here, such as the KOhmms, also including eukaryotic genomes. A good alternative to add are the [KOGs](https://pubmed.ncbi.nlm.nih.gov/14759257/) that are specifically designed for eukaryotes. The most recent profiles provided by EggNOG can be found [here](http://eggnog5.embl.de/#/app/downloads). 


## Version programs
##################################

The versions might differ from what you have available on your system and should not be problematic but for consistency we list the version numbers that were used to generate this workflow as well as links with install instructions.

1. Python 2.7.5
2. perl v5.16.3
3. prokka 1.14-dev ; install instructions can be found [here](https://github.com/tseemann/prokka)
4. HMMER 3.1b2; install instructions can be found [here](http://hmmer.org/documentation.html)
5. blastp v2.7.1+; install instructions can be found [here](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download)
6. diamond 0.9.22; install instructions can be found [here](https://github.com/bbuchfink/diamond/wiki)
7. interproscan: InterProScan-5.29-68.0; install instructions can be found [here](https://interproscan-docs.readthedocs.io/en/latest/)
8. GNU parallel 20181222; install instructions can be found [here](https://www.gnu.org/software/parallel/)



## Version databases
##################################

Databases (and mapping files) used in the current version can be downloaded from [here] (https://zenodo.org/deposit/6489670).

1. KOhmms: downloaded 90419 from [here](https://www.genome.jp/tools/kofamkoala/)
2. PGAP (former TIGRs): downloaded Nov 2021 from [here](https://ftp.ncbi.nlm.nih.gov/hmm/current/)
3. Pfams: RELEASE 34.0, downloaded Nov 2021 from [here](http://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/)
4. CAZymes: dbCAN-HMMdb-V9 downloaded on 21 April 2020 from [here](https://bcb.unl.edu/dbCAN2/download/)
5. Merops: downloaded from Merops server in November 2018 from [here](https://www.ebi.ac.uk/merops/download_list.shtml)
6. HydDB: downloaded form HydDB webserver in July 2020 from [here](https://services.birc.au.dk/hyddb/browser/)
7. COGs: downloaded from NCBI Octover 2020 from  [here](ftp://ftp.ncbi.nlm.nih.gov/pub/COG/COG2020/)
8. TransporterDB: downloaded from TCDB April 2021 [here](http://www.tcdb.org)
9. ncbi_nr (maintained and downloaded by Alejandro Abdala Asbun). Last update was done on August 2021
10. ArCOGs: downloaded in 2019 and available [here](https://ftp.ncbi.nih.gov/pub/wolf/COGs/arCOG/) (update date: 2018)



## External scripts needed
##################################

Custom scripts used in the current version can be downloaded from [here] (https://zenodo.org/deposit/6489670).

**perl**

- length+GC.pl

**python**

- Parse_prokka_for_MAGs_from_gbk-file.py
- Split_Multifasta.py
- parse_IPRdomains_vs2_GO_2_ts_sigP.py
- parse_diamond_blast_results_id_taxid.py



###################################################################################
###################################################################################
# Setup working environment
###################################################################################
###################################################################################


##  Define variables
##################################

Notice: If you have done the assembly and binning tutorials, you can see that we do things a bit differently here. Since for the annotations we work with a lot of different databases, there are two ways how to deal with a lot of files in different locations: 

(a) we can add the path for each database into the code itself  
(b) we can define variables and store the file paths inside these variables

Especially if you work with collaborators it often is easier to use the second option, as you just have to adjust the paths in one spot.

<u>Instructions: </u>

- Change the first variable, the wdir, to the path were you want to work from and change the paths for the databases according to where you downloaded the databases.
- The two last lines of code are needed for the software prokka and are specific for the system the script was written on. Adjust according to how prokka is setup in your system.


```{bash class.source = 'fold-show', eval=FALSE} 
#set the path to our working directory
wdir="/export/lv1/user/ndombrowski/Mai_testing"

#nr of cpus we generally want to use
cpus=5

#location of some custom python and perl scripts
python_scripts="/export/lv1/user/spang_team/Scripts/Others/"
perl_scripts="/export/lv1/user/spang_team/Scripts/Others/"

#locations of the indiv. databases
cog_db="/export/data01/databases/cogs/hmms/NCBI_COGs_Oct2020.hmm"
arcog_db="/export/data01/databases/arCOG/arCOGs2019/All_Arcogs_2018.hmm"
ko_db="/export/data01/databases/ko/combined_hmms/All_KOs.hmm"
pfam_db="/export/data01/databases/pfam/Pfam-A.hmm"
tigr_db="/export/data01/databases/PGAP_ProtFam/hmm_PGAP.lib"
cazy_db="/export/data01/databases/cazy_spt/dbCAN-HMMdb-V7.txt"
transporter_db="/export/data01/databases/tcdb/tcdb"
hyd_db="/export/data01/databases/HydDB/HydDB_uniq"
ncbi_nr_db="/export/data01/databases/ncbi_nr/diamond/nr.dmnd"

#locations of the mapping files for the indiv. databases (links the database ID to a more useful description)
cog_mapping="/export/data01/databases/cogs/hmms/cog_mapping.txt"
arcog_mapping="/export/lv1/user/spang_team/Databases/arCOGs2019/ar14.arCOGdef18.final.tab"
ko_mapping="/export/data01/databases/ko/ko_list_for_mapping_clean"
pfam_mapping="/export/data01/databases/pfam/Pfam-A.clans.cleaned.tsv"
tigr_mapping="/export/lv1/user/spang_team/Databases/PGAP/hmm_PGAP_edited.tsv"
cazy_mapping="/export/data01/databases/cazy_spt/CAZY_mapping_2.txt"
transporter_mapping="/export/lv1/user/spang_team/Databases/TransporterDB/tcdb_definition_2021.tsv"
hyd_mapping="/export/data01/databases/HydDB/HydDB_mapping"
taxon_db="/export/data01/databases/taxmapfiles/ncbi_nr/prot.accession2taxid.gz"

#set environment to work for prokka
export PERL5LIB=/export/lv1/public_perl/perl5/lib/perl5

#add libraries so that tbl2asn works
export PATH=$PATH:/opt/biolinux/prokka/binaries/linux/
```




##  Set our working directory
##################################

Let's first go to our working directory and make a new folder for the Annotation tutorial.

```{bash class.source = 'fold-show', eval=FALSE}
#check that we set the variable to our working directory correctly
echo $wdir

#change the directory to your home directory
cd $wdir

#make a folder for the Annotations and move into this directory
mkdir Annotations
cd Annotations
```


###################################################################################
###################################################################################
# Prepare genomes of interest
###################################################################################
###################################################################################

In this workflow we will work with 2 genomes that we download from NCBI but the steps can easily be adjusted for your own genomes. You could also work with the bins we have assembled and binned in the previous workflows, if you feel comfortable with the command line feel free to use those instead.

First, we will download two archaeal genomes from NCBI and rename them a bit to avoid issues with the workflow. The bins from the previous workflows (inside ../Binning/4_FinalBins/renamed) already have a good header and should be ready to use in case you want to repeat the workflow with those genomes.

After we have downloaded the 2 genomes from NCBI the file name and file header should look something like this:

```{r fig1, echo=FALSE, fig.cap="", out.width = '70%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations1b.png")
```

We can see that both the file name and file header are long and has a lot of symbols UNIX might have trouble with (brackets, spaces, dots...) and generally, something to watch out for when naming files is:

- Have short but descriptive file names and file headers. I.e. if we look at the header of the NCBI genomes you will see the header is rather long, so it is useful to condense it. <br/><br/>

- Avoid extra symbols (space, comma, semicolon, quotes, brackets) as these can easily result in unwanted behavior. I.e. a very simple naming could only include two symbols:

  - underscores, `_` , for everything where we normally would use a space
  - minus symbol  `-`  for key information that we might want to separate later (i.e. genome name and protein ID)<br/><br/>

Another useful thing is to add information into the headers that allow us to concatenate multiple genomes and afterwards still know from what genome a sequence came from. Therefore, we recommend to:

  - Add the genome name into the fasta header and use a unique delimiter to separate it from the contig id. That way you can concatenate your genomes and afterwards still know from what genome a contig came from
- have a consistent way to separate the genome and protein ID across different analyses. In phylogenies its very common to concatenate several genes/proteins from individual genomes and for this it is useful to remove the proteinID from the fasta header.

```{bash class.source = 'fold-show', eval=FALSE}
#make folders in which we store the genomes
mkdir fna
mkdir fna/renamed

#download genomes from NCBI using wget
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/008/085/GCA_000008085.1_ASM808v1/GCA_000008085.1_ASM808v1_genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/017/945/GCA_000017945.1_ASM1794v1/GCA_000017945.1_ASM1794v1_genomic.fna.gz

#unzip files and check the file names
gzip -d *gz
ls -l *.fna

#shorten the file name by cutting everything away after the dot
for f in *.fna; do mv $f "${f/.*}".fna; done

#check the file name again
ls -l *.fna

#check how the header looks --> we have a long name and a lot of symbols we do not want, lets remove these
head GCA_000008085.fna 

#shorten the fasta headers by removing everything after the space
for f in *.fna; do cut -d " " -f1  $f > fna/${f}; done
head fna/GCA_000008085.fna 

#add the genome name into the header (this is useful if we want to concatenate a lot of genomes into one file)
cd fna
for i in *fna; do awk '/>/{sub(">","&"FILENAME"-");sub(/\.fna/,x)}1' $i > renamed/$i; done
head renamed/GCA_000008085.fna 

#remove intermediate files
rm *fna
cd ..
rm *fna
```

The cleaned up files and file names should look something like this afterwards:

```{r fig2, echo=FALSE, fig.cap="", out.width = '70%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations2b.png")
```

<br/><br/>

### The steps above work fine but what if we have hundreds of genomes from NCBI?

<div class="alert-danger">
For the tutorial we will NOT run the code below but continue at the step: **Make a filelist for all genomes of interest**. </br></br>
However, if you deal a lot with genomes from NCBI and you want to have a file that links a ncbi tax id to a tax string, this still might be useful to know. I.e. if you work with a genome of interest, such as E.coli, you might want to include several closely related genomes as reference in your analysis and NCBI is an excellent resource providing you with a multitude of genomes as well as metadata.
</div>


Notice: 

- The file that links the taxonomy ID with the taxa string (lineages-2020-12-10.csv) was generated with ncbitax2lin (instructions added below)
- this part of the workflow was set up to run at servers at the NIOZ and the paths need to be adjusted to work on your system

The first part of this workflow is to get a list of all available archaeal (or bacterial) genomes from ncbi and add a taxa string into the file we downloaded. 


```{bash class.source = 'fold-show', eval=FALSE}
#prepare some folders
mkdir ncbi_genomes
cd ncbi_genomes

#get a list of all avail. archaeal genomes
wget https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt
```

This file includes a list of all (archaeal in this case) genomes available at NCBI and relevant metadata such as the project ID, taxonomy ID, species name, download link for individual genome files etc.

Next, we will clean the file a bit and since the species name might not be informative enough (i.e. we might want to know to what phylum a species belongs), we want to add the full taxonomy string in.

```{bash class.source = 'fold-show', eval=FALSE}
#clean up the file by removing the first row and any `#` and move the tax id to the front
sed 1d assembly_summary_genbank.txt | sed 's/# //g' | awk -F"\t" -v OFS="\t" '{print $6, $0} ' | sed 's/ /_/g' > assembly_summary_edit.txt

#replace empty rows with a "none" and remove the header
awk 'BEGIN { FS = OFS = "\t" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i="none"}; 1' assembly_summary_edit.txt | sed '1d' > temp1

#add in the tax string (to be able to better select genomes)
#the code for how to generate the csv needed here can be found below
#merge assembly info with tax string
LC_ALL=C join -a1  -j1 -e'-' -o 0,1.2,1.3,1.4,1.5,1.6,1.8,1.9,1.10,1.11,1.13,1.16,1.17,1.18,1.19,1.21,2.2 <(LC_ALL=C sort temp1) <(LC_ALL=C sort  /export/lv1/user/spang_team/Databases/NCBI_taxonomy/Dec2020/lineages-2020-12-10.csv) -t $'\t' > temp2

#add in headers
echo -e "taxid\taccession\tbioproject\tbiosample\twgs\trefseq_category\tspeciesID\torganism\tintraspecific_name\tisolate\tassembly\trel_date\tasm\tsubmitter\tgcf_name\tftp_path\ttax" | cat -  temp2 >  assembly_summary_archaea_Dec2020_tax.txt

cd ..
```      

In this file you can just look for the genomes you are interested in. You can subset for lineages of interest by selection specific groups via their taxonomy and you might run CheckM as well to select them based on their genome stats.

Once you have your genomes selected that you are interested in you then can copy the ftp_links into a document called `get_genomes.txt` via nano and then do some small changes to the file with sed. Afterwards you can download all the genomes from NCBI.

```{bash class.source = 'fold-show', eval=FALSE}
# extend path to be able to download the genomes
sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/)(GCA/)([0-9]{3}/)([0-9]{3}/)([0-9]{3}/)(GCA_.+)|\1\2\3\4\5\6/\6_genomic.fna.gz|' get_genomes.txt  > get_genomes_v2.txt

# download genomes
for next in $(cat get_genomes_v2.txt); do wget  "$next"; done
```

Once you have the genomes, you can clean the names as suggested above and continue with the steps below.

In case you want to generate the csv file that links the ncbi taxID to taxstring, which includes the phylum, class, order, etc information, this is the workflow (you will need to download [ncbitax2lin](https://github.com/zyxue/ncbitax2lin)):

```{bash class.source = 'fold-show', eval=FALSE}
#change to the directory in which we installed ncbitax2lin
cd /export/lv1/user/spang_team/Scripts/NCBI_taxdump/ncbitax2lin

#start bashrc that allows to run conda
source ~/.bashrc.conda2

#setup virtual environment (see more instructions here: https://github.com/zyxue/ncbitax2lin)
#virtualenv venv/

#re-do taxonomy mapping file
make

#unzip file
gzip -d  lineages-2020-12-10.csv.gz

#close virtual environment
#source deactivate venv/

#make a new directory for the most current tax
mkdir /export/lv1/user/spang_team/Databases/NCBI_taxonomy/Dec2020

#change dir
cd /export/lv1/user/spang_team/Databases/NCBI_taxonomy/Dec2020

#cp the taxa file into our new dir and replace spaces with underscore
sed 's/ /_/g'  ~/../spang_team/Scripts/NCBI_taxdump/ncbitax2lin/lineages-2020-12-10.csv  > temp1

#only print the levels until species level
awk -F',' -v OFS="\t" '{print $1,$2,$3,$4,$5,$6,$7,$8}' temp1 > temp2

#add in “none” whenever a tax level is emtpy
awk 'BEGIN { FS = OFS = "\t" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i="none"}; 1' temp2 > temp3

#merge columns 2-8, thus we have 1 column with the tax number and one column with the tax string
awk ' BEGIN { FS = OFS = "\t" } {print $1,$2","$3","$4","$5","$6","$7","$8}' temp3 > temp4

#remove header
sed '1d' temp4 > lineages-2020-12-10.csv

#clean up
rm temp*


```


## Make a filelist for all genomes of interest

Make a file that lists all the genomes we are working with. Such a file list is useful when working with for loops.

```{bash class.source = 'fold-show', eval=FALSE}
#make a folder for all the secondary files we need
mkdir FileLists

#list our genomes of interest
ls fna/* | sed 's/\.fna//g' | sed 's/fna\///g' > FileLists/GenomeList

#check content of file
head FileLists/GenomeList
```

The last command should print something like this:

```{r, echo=FALSE, fig.cap="", out.width = '50%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations7.png")
```




## Call proteins with prokka and clean files


Now that we have the genomes organized, lets identify the coding regions and call the proteins with Prokka. Another common tool you might encounter for this step is prodigal. Prokka basically uses prodigal but is doing some extra annotations, which we also want to look at.

Notice:

- For Prokka we have some specific setups that are needed for Prokka to run on the servers this workflow was written for. If you are on your own system adjust accordingly.(the export functions are part of the start of the script, were we define all the additional environmental variables).
- Our genomes are archaeal, if you have other genomes change the **--kingdom Archaea** part below. This setting has no major differences for calling proteins but the Prokka annotations might slightly change.
- Prokka is not designed to work with eukaryotic genomes, an alternative to look into for eukaryotes is for example Augustus.
- Prokka does very quick and dirty annotations, these are good for a first glimpse into what you have but in our experience there can be several issues, so use them with caution.


### Run prokka

[Prokka](https://github.com/tseemann/prokka): a tool for rapid prokaryotic genome annotation [(Seeman et al., 2014)](https://pubmed.ncbi.nlm.nih.gov/24642063/).

```{bash class.source = 'fold-show', eval=FALSE}
#prepare new directories to store our data
mkdir faa
mkdir faa/Prokka

#run prokka
for sample in `cat FileLists/GenomeList`; do prokka fna/renamed/${sample}.fna --outdir faa/Prokka/${sample} --prefix $sample --kingdom Archaea --addgenes --force --increment 10 --compliant --centre UU --cpus 2 --norrna --notrna ; done

#control number of Genomes
ll faa/Prokka/*/*faa | wc -l

#check how the genome header looks like
head faa/Prokka/GCA_000008085/GCA_000008085.faa
```

We have 2 genomes and the header of the first file should look something like this:

```{r, echo=FALSE, fig.cap="", out.width = '60%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations8.png")
```

We can see that prokka gives a unique protein ID and adds an annotation to the file header.

Additionally, Prokka generates several output files, these include:

<p align="left">
  <img width=450, height=450, src="/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/prokka.png">
</p>

Additonally, prokka does not only call proteins but also identifies tRNAs, does protein anntations, etc.

<p align="left">
  <img width=600, height=600, src="/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/prokka2.png">
</p>




### Concatenate genomes

If we look at the fasta headers of the newly generated files, we see that the header is not very descriptive. I.e. we have something like this:   
``>CALJACFG_00010 ATP-dependent DNA helicase Hel308``.<br/><br/>
This header is very long and uses a lot of extra symbols (a space, and a minus symbol, which we already use for other purposes). Therefore, we first want to clean things up a bit.

Additionally, since for the following workflow we want to concatenate the two genomes, i.e. merge the two files into one, it is useful to add the genome name into the fasta header before merging the files. 

This is optional, but it is useful to use a special delimiter to later separate the genome from the protein name. Therefore we add the file name and separate it with the old sequence header with a `-`.


```{bash class.source = 'fold-show', eval=FALSE}
#go into the prokka folder
cd faa/Prokka/

#organize files
mkdir renamed
cp GCA*/*faa .

#add the filename into the fasta header
for i in *faa; do awk '/>/{sub(">","&"FILENAME"-");sub(/\.faa/,x)}1' $i > renamed/$i; done

#cleanup temp files
rm *faa

#concatenate files
cat renamed/*faa > All_Genomes.faa

#clean the header/shorten it by removing everything after the space
cut -f1 -d " "  All_Genomes.faa > All_Genomes_clean.faa

#go back to the wdir
cd ../..

#check the fasta header
head faa/Prokka/All_Genomes_clean.faa 

#make a blast db out of concat genomes (we need this whenever we want to run blast)
makeblastdb -in faa/Prokka/All_Genomes_clean.faa  -out faa/All_Genomes_clean   -dbtype prot -parse_seqids
```

Now our file header should look something like this: 

```{r, echo=FALSE, fig.cap="", out.width = '50%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations9.png")
```



###################################################################################
###################################################################################
# If you work with your own genomes
###################################################################################
###################################################################################

- Ideally have not only the protein ID but also the genome ID sequence header
- make sure the file header is concise and does not have ANY spaces and that ideally uses a '-' to separate the binID from the protein ID
- If you have a concise header + the bin ID in the header, first concatenate the protein sequences of your genomes in one file.


More specifically you want to:

(a) for the fna files
  - generate a folder named `fna/renamed`
  - add the individual genomes in there


(b) for the faa files:
  - generate the file in the folder `faa/Prokka`
  - concatenate the proteins from all your genomes of interest and name the file `All_Genomes_clean.faa`
  - continue with the steps below



###################################################################################
###################################################################################
# Make a mapping file with all relevant genome info
###################################################################################
###################################################################################


Now that our files are ready we first want to make a table with the basic stats (i.e. contig length, contig gc, protein length, ...).



## Get contig stats
##################################

Prokka has the habit to rename the fasta header. I.e. the original contig name is renamed.
Therefore, it is hard to link from what contig the proteins originally came from. The workflow below generates a mapping file that links the original contig name with the prokka contig name.

The goal is to generate a file with this information:

- old contig name, binID, contig GC, contig length

**Notice**

- This part depends a lot on the files having short naming schemes and no extra symbols. So make sure that the final file, has the 4 columns we would expect.

<div class="alert-danger">
**Comment if working with your own genomes:**<br/><br/>

If working with your own genomes, esp. step 4 can be problematic. This step will ONLY work if your fasta header includes the binID and the contig ID that are separated with a `-`. If you do not have both pieces of info, have extra symbols symbol or if that symbol is somewhere else in the filename/header this part of the workflow will very likely not work (without needing adjustments). 
<br/><br/>
In this case  it will be easiest skip this step and do the following:<br/><br/>
(A) generate the file FileLists/1_Bins_to_protein_list.txt in the section **Prepare list that links proteins with binIDs** <br/><br/>
(B) continue with the step **Do the annotations**.  <br/><br/>
(C) During the first annotation step with the COG database, ignore the step #combine with previous data frame and do `cp NCBI_COGs/A_NCBI_COGs.tsv merge/temp_A.tsv` to generate the first table.
</div>


```{bash class.source = 'fold-show', eval=FALSE}
#make folder
mkdir contig_mapping

#for each contig list the gc and contig length (we use this later to control whether the merging of information works ok)
for sample in `cat FileLists/GenomeList`; do perl $perl_scripts/length+GC.pl fna/renamed/${sample}.fna > contig_mapping/${sample}_temp1; done

#add number of contigs as separate column
for sample in `cat FileLists/GenomeList`; do awk '$1=(FNR FS $1 FILENAME)' contig_mapping/${sample}_temp1 > contig_mapping/${sample}_temp2; done

#cleanup file: add in binIDs, so we have: contig id, binID, contig nr, contig gc and contig length
for sample in `cat FileLists/GenomeList`; do awk 'BEGIN{OFS="\t"}{split($2,a,"-"); print a[2],$1,$3,$4}' contig_mapping/${sample}_temp2  > contig_mapping/${sample}_temp3; done

for sample in `cat FileLists/GenomeList`; do awk 'BEGIN{OFS="\t"}{split($1,a,"/"); print a[1],a[2], $2,$3,$4}' contig_mapping/${sample}_temp3 | sed 's/contig_mapping//g' | sed 's/_temp1//g'  > contig_mapping/${sample}_temp4; done

#combine results
cat contig_mapping/*temp4 > contig_mapping/Contig_Old_mapping.txt

#check file. heder should have --> contigID, binID, contig NR, GC, contig length
head contig_mapping/Contig_Old_mapping.txt

#cleanup temp files
rm contig_mapping/*temp*

#create file for merging with prokka contigIDs. Here, we add the binID with the contig number together to recreate the prokka contig name
awk 'BEGIN{OFS="\t"}{print $2"_contig_"$3,$0}' contig_mapping/Contig_Old_mapping.txt > contig_mapping/Contig_Old_mapping_for_merging.txt

#check file: prokka contig name, prokka binID, binID, contig NR, GC, contig length
head contig_mapping/Contig_Old_mapping_for_merging.txt

```

In the file we see the contig name, contig GC and contig length. Since the two genomes are complete genomes we only work with two contigs in total. If you work with MAGs this table will get much longer.

```{r, echo=FALSE, fig.cap="", out.width = '60%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations10.png")
```





## Prepare list that links proteins with binIDs
##################################

Now, we want to generate a file that lists all proteins we are working and a second column with the bin ID. 
So we want two columns: proteinID (accession), binID. 

```{bash class.source = 'fold-show', eval=FALSE}
#get list of protein accession numbers
grep "^>" faa/Prokka/All_Genomes_clean.faa  > temp

#remove the prokka annotations
cut -f1 -d " " temp > temp2

#remove the ``>``
sed 's/>//g' temp2 > FileLists/AllProteins_list.txt

#check file; we want one column with the proteinIDs
head FileLists/AllProteins_list.txt 

#remove temp files
rm temp*

#Modify protein list to add in a column with binID
awk -F'\t' -v OFS='\t' '{split($1,a,"-"); print $1, a[1]}' FileLists/AllProteins_list.txt | LC_ALL=C sort > FileLists/1_Bins_to_protein_list.txt

#check file content: we now have the accession and the binID
head FileLists/1_Bins_to_protein_list.txt

#check with how many proteins we work
wc -l FileLists/1_Bins_to_protein_list.txt
```

We now have a list with 2,069 proteins and the respective genome each protein belongs to. The first few lines of the file should look something like this:

```{r, echo=FALSE, fig.cap="", out.width = '70%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations11.png")
```

Throughout the workflow (or whenever generating your own pipeline), check that the number of protein stays consistent.

Additionally, while running things for the first time it is always good to check the files we generate to get a better feeling on:

- what is happening at specific points in the script
- Find potential bugs



## Prepare a file with prokka annotations & Lengths & other genome info
##################################

Now, we want to generate a file that links the contig info and all the protein info.

```{bash class.source = 'fold-show', eval=FALSE}
#compress gbk files (generated by prokka)
gzip faa/Prokka/*/*gbk

#make list for looping
for f in faa/Prokka/*/*gbk.gz; do echo $f >> FileLists/list_of_gzips; done

#control that the nr of files we have is correct
wc -l FileLists/list_of_gzips

'''
We deal with 2 genomes
'''

#create a summary file based on the gbk file
#lists the binname, contig name, contig length, proteinID, prokka annotation
python $python_scripts/Parse_prokka_for_MAGs_from_gbk-file_v2.py -i FileLists/list_of_gzips -t FileLists/Contig_Protein_list.txt

#check output of the script: Prokka binID, Prokka contig ID, contig length, gene start, gene stop, strand, accession, prokka annotation
head FileLists/Contig_Protein_list.txt
```

The output should look something like this

```{r, echo=FALSE, fig.cap="", out.width = '90%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations12.png")
```

Next, we want to have a list that assigns the new ID prokka gives to a genome to the old genome name we got from NCBI:

```{bash class.source = 'fold-show', eval=FALSE}
#make prokka to binId mapping file
awk 'BEGIN{OFS="\t"}{split($1,a,"-"); print a[2],$2}' FileLists/1_Bins_to_protein_list.txt | sed 's/_[0-9]*\t/\t/g' | sort | uniq > FileLists/Prokka_to_BinID.txt

#check file: Prokka bin ID, old bin ID
head FileLists/Prokka_to_BinID.txt
```

```{r, echo=FALSE, fig.cap="", out.width = '50%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations13.png")
```

Finally, we can merge everything together.

```{bash class.source = 'fold-show', eval=FALSE}
#link old to new binIDs
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' FileLists/Prokka_to_BinID.txt FileLists/Contig_Protein_list.txt | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$10,$2,$2,$3,$4,$5,$6,$7,$8}'  > temp_Bin_Contig_Protein_list.txt

#add in extra column for contig nr
awk -F'\t' -v OFS='\t'   'gsub("[A-Za-z0-9]*_","",$4)' temp_Bin_Contig_Protein_list.txt | awk -F'\t' -v OFS='\t' '{print $2"_contig_"$4,$1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' > FileLists/Bin_Contig_Protein_list.txt

#check file
head FileLists/Bin_Contig_Protein_list.txt

#merge with old contig IDs
#headers: accession, BinID, newContigID, oldContigID, mergeContigID,ContigLengthNew, LengthContigOld, GC,ProteinID, protein start, protein end, protein strand, prokka annotation
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' contig_mapping/Contig_Old_mapping_for_merging.txt FileLists/Bin_Contig_Protein_list.txt | awk 'BEGIN{FS="\t";OFS="\t"}{print $3"-"$7,$3,$4,$13,$1,$6,$17,$16,$7,$8,$9,$10,$11 }' >  FileLists/Bin_Contig_Protein_list_merged.txt

#check file
head FileLists/Bin_Contig_Protein_list_merged.txt

#prepare a file with protein length and gc to add protein length into our file
perl $perl_scripts/length+GC.pl faa/Prokka/All_Genomes_clean.faa > temp

#merge with contig/protein info
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' temp FileLists/Bin_Contig_Protein_list_merged.txt | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,$4,$6,$7,$8,$9,$10,$11,$12,$16,$13}' > temp2

#add a header
echo -e "accession\tBinID\\tNewContigID\tOldContigId\tContigNewLength\tContigOldLength\tGC\tProteinID\tProteinStart\tProteinEnd\tStrand\tProteinLength\tProkka" | cat - temp2 > FileLists/B_GenomeInfo.txt

#check file
less -S FileLists/B_GenomeInfo.txt

#control that all is ok
wc -l FileLists/B_GenomeInfo.txt
#2070

#remove temp files
rm temp*

```

The file stored in `FileLists/B_GenomeInfo.txt` should look something like this:

<p>
  <img width=900, height=900, src="/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/annotations14.png">
</p>



###################################################################################
###################################################################################
# Do the annotations
###################################################################################
###################################################################################

All searches work have a similar workflow:

- Run the search (with hmmer, diamond, blast, etc.)
- For each protein, make sure we have just one unique hit (usually we select the best hit for each protein using the best e-value or bitscore). I.e. a single protein might be similar to several dehydrogenases and we only keep the best hit.
- Clean up the output from the search.
- Merge the output with our protein list (that way we make sure we always have the same number of rows in all the files we generate).
- add gene descriptions to our table, i.e. when we do the search against the COGs we only have the COG IDs, but we want to add a gene description here.
- add an informative header.

For all the searches, we keep the original input table and the final output table and the other intermediate files we delete.

Once we have run our searches against several different databases, we merge all the results in a step-wise manner (using the protein IDs) and then in the final step clean the final data frame by removing redundant information.

We will add basic information on the different databases and search algorithms in the last section of this tutorial, so if you want more details check there.


## COG search
##################################

The COGs (Clusters of Orthologous Genes) is a database of orthologous gene clusters generated from a database from 1,187 bacteria and 122 archaea (2020 version) and therefore is useful if you work with organisms from these two groups.

Possible settings to change:

- number of CPUs 
- e-value threshold 

```{bash class.source = 'fold-show', eval=FALSE} 
#1. Setup directories
mkdir NCBI_COGs
mkdir merge

#2. run hmmsearch against all COGs 
hmmsearch --tblout NCBI_COGs/sequence_results.txt -o NCBI_COGs/results_all.txt --domtblout NCBI_COGs/domain_results.txt --notextw --cpu $cpus $cog_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select hits above a certain e-value
sed 's/ \+ /\t/g' NCBI_COGs/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $3, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-3'  > NCBI_COGs/sequence_results_red_e_cutoff.txt

#get best hit/protein based on bit score, and e-value
sort -t$'\t' -k3,3gr -k4,4g NCBI_COGs/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  NCBI_COGs/All_NCBI_COGs_hmm.txt

#merge with contig list
LC_ALL=C join -a1  -j1 -e'-' -t $'\t' -o 0,2.2,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort NCBI_COGs/All_NCBI_COGs_hmm.txt) | LC_ALL=C sort  > NCBI_COGs/temp4

#merge with COG mapping file 
LC_ALL=C join -a1 -1 2 -2 1 -e'-' -t $'\t' -o1.1,0,2.3,2.2,2.5,1.3 <(LC_ALL=C sort -k2 NCBI_COGs/temp4) <(LC_ALL=C sort -k1 $cog_mapping) | LC_ALL=C  sort > NCBI_COGs/temp5

#add headers
echo -e "accession\tNCBI_COG\tNCBI_COG_Description\tPathwayID\tPathway\tNCBI_COG_evalue" | cat - NCBI_COGs/temp5 > NCBI_COGs/A_NCBI_COGs.tsv

#combine with previous data frame
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' NCBI_COGs/A_NCBI_COGs.tsv  FileLists/B_GenomeInfo.txt > merge/temp_A.tsv

#check file
less -S merge/temp_A.tsv

#clean-up
rm NCBI_COGs/temp*
````


In step 2, hmmsearch will generate multiple outputs (for details on the individual data found in these tables [see chaper 5, page 45](http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf)) of the HMMER manual:

1. The --tblout output option produces the target hits table. The target hits table consists of one line for each different query/target comparison that met the reporting thresholds, ranked by decreasing statistical significance (increasing E-value).
2. The --domtblout option produces the domain hits table. There is one line for each domain. There may be more than one domain per sequence. T
3. A table with the raw output, including alignments (-o option). Notice: This file can get rather large and if you have space issues, its not essential and can be deleted.

After parsing our table that is stored in `NCBI_COGs/A_NCBI_COGs.tsv` contains the protein accession, COG ID, COG description and e-value. If other relevant data is available, in the case of the COGs individual COGs are sorted into individual pathways, this information is added as well:

```{r, echo=FALSE, fig.cap="", out.width = '100%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations15.png")
```



**Exercise**

Compare the files sequence_results.txt and sequence_results_red_e_cutoff.txt. What is the difference between them?

<div id="spoilerA" style="display:none">

The commands to check this are:

```{bash class.source = 'fold-show', eval=FALSE}
head NCBI_COGs/sequence_results.txt<br/>
head NCBI_COGs/sequence_results_red_e_cutoff.txt <br/><br/>
```

We can see that we removed some columns from the second table and that we excluded any rows starting with a #

</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoilerA') .style.display=='none')
              {document.getElementById('spoilerA') .style.display=''}
            else{document.getElementById('spoilerA') .style.display='none'}">
  Show/hide
</button>
<br/><br/>



## ArCOGs search
##################################

The arCOGs is a database of orthologous groups based on archaeal genomes and therefore is especially useful when working with archaea.

Possible settings to change:

- number of CPUs (default = 2)
- e-value threshold (default = 1e-3)


```{bash class.source = 'fold-show', eval=FALSE}
#prepare folders
mkdir arcogs

#run hmmsearch
hmmsearch --tblout arcogs/sequence_results.txt -o arcogs/results_all.txt --domtblout arcogs/domain_results.txt --notextw --cpu $cpus $arcog_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select sequences above a certain e-value
sed 's/ \+ /\t/g' arcogs/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $3, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-3'  > arcogs/sequence_results_red_e_cutoff.txt

#get best hit based on bit score, and then e-value
sort -t$'\t' -k3,3gr -k4,4g arcogs/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  arcogs/All_arcogs_hmm.txt

#separate header with arcog and ID
awk -F'\t' -v OFS='\t' '{split($2,a,"."); print $1, a[1], $3,$4}' arcogs/All_arcogs_hmm.txt | LC_ALL=C sort > arcogs/temp3

#merge with contig list
LC_ALL=C join -a1  -j1 -e'-' -t $'\t' -o 0,2.2,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort arcogs/temp3) | LC_ALL=C sort  > arcogs/temp4

#merge with arcog names
LC_ALL=C join -a1 -1 2 -2 1 -e'-' -t $'\t' -o1.1,0,2.3,2.4,2.2,1.3 <(LC_ALL=C sort -k2 arcogs/temp4) <(LC_ALL=C sort -k1 $arcog_mapping) | LC_ALL=C  sort > arcogs/temp5

#add in headers
echo -e "accession\tarcogs\tarcogs_geneID\tarcogs_Description\tPathway\tarcogs_evalue" | cat - arcogs/temp5 > arcogs/C_arcogs.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' arcogs/C_arcogs.tsv merge/temp_A.tsv > merge/temp_BC.tsv

#check what happened
less -S merge/temp_BC.tsv 

#clean-up
rm arcogs/temp*
```


Now, we have in our full table not only the information for the COGs but also the arCOGs. This way we can go to every individual protein and check whether they hit to the same protein.

```{r, echo=FALSE, fig.cap="", out.width = '100%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations16.png")
```

In general it is useful to use more then one database to analyse your genome since no database is perfect and having several databases to compare the results allows us to judge how confident we are with a certain annotation.

Unfortunately, they is no consistent naming convention. I.e. in the table above you can see for CALJACFG_00080 that the COG link this protein to a tRNA_A37_methylthiotransferase and the arCOG links this protein to a 2-methylthioadenine_synthetase but despite the different names both names actually are for the same enzyme.

**Question**

We now have 3 annotations for each protein from prokka, the arcogs and the COGs.  
Look at the gene with ID 00260 for the genome GCA_000008085 and check if all 3 databases give consistent annotations. If not, which annotation do you trust most?   
Since the prokka generates different BinIDs for each run, I can not give you the exact accession, but using grep to search for the genome and protein ID should be enough for you to find the protein.

[CALJACFG_00260 hits to sulfide-quinone reductase with prokka and to a NADH_dehydrogenase with the arcogs and the COGs. Generally, the e-values are pretty good with the last two databases, while prokka by default does not give us any score to judge the confidence of this hit. We therefore should be skeptical about the annotation as sulfide quinone reductase.]{.spoiler}






## KO search using Hmm's
##################################

The KO (KEGG Orthology) database is a database of molecular functions represented in terms of functional orthologs and was generated using archaeal, bacterial and eukaryotic genomes.

Possible settings to change:

- number of CPUs 
- e-value threshold 


```{bash class.source = 'fold-show', eval=FALSE}
#setup directories
mkdir KO_hmm

#run KO search 
hmmsearch --tblout KO_hmm/sequence_results.txt -o KO_hmm/results_all.txt --domtblout KO_hmm/domain_results.txt --notextw --cpu $cpus $ko_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select sequences above a certain e-value
sed 's/ \+ /\t/g' KO_hmm/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $3, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-3'  > KO_hmm/sequence_results_red_e_cutoff.txt

#get best hit based on bit score, and then e-value
sort -t$'\t' -k3,3gr -k4,4g KO_hmm/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  KO_hmm/All_KO_hmm.txt

#merge with protein list
LC_ALL=C join -a1 -t $'\t' -j1 -e'-' -o 0,2.2,2.3,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort KO_hmm/All_KO_hmm.txt) | sort -u -k1 > KO_hmm/temp

#merge with KO_hmm names
LC_ALL=C join -a1 -1 2 -2 1 -e'-' -t $'\t' -o 1.1,1.2,1.4,1.3,2.2,2.12  <(LC_ALL=C sort -k2 KO_hmm/temp) <(LC_ALL=C sort -k1 $ko_mapping) | LC_ALL=C  sort > KO_hmm/temp3

#add in an extra column that lists whether hits have a high confidence score
awk  -v OFS='\t' '{ if ($4 > $5){ $7="high_score" }else{ $7="-" } print } ' KO_hmm/temp3 > KO_hmm/temp4

#add header
echo -e "accession\tKO_hmm\te_value\tbit_score\tbit_score_cutoff\tDefinition\tconfidence" | cat - KO_hmm/temp4 > KO_hmm/M_KO_hmm.tsv

#combine with previous data frame
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' KO_hmm/M_KO_hmm.tsv merge/temp_BC.tsv > merge/temp_BCM.tsv

#control lines
wc -l merge/*

#clean up
rm KO_hmm/temp*
```







## Pfam search
##################################

The Pfam database is a large collection of protein families that are generated from [UniProt Reference Proteomes](https://www.uniprot.org/help/reference_proteome) and includes archaea, bacteria and eukaryotes.

Possible settings to change:

- number of CPUs
- e-value threshold


```{bash class.source = 'fold-show', eval=FALSE}
#1. Setup directories
mkdir pfam

#2. run hmmsearch against all pfams
hmmsearch --tblout pfam/sequence_results.txt -o pfam/results_all.txt --domtblout pfam/domain_results.txt --notextw --cpu $cpus $pfam_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select sequences above a certain e-value
sed 's/ \+ /\t/g' pfam/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $4, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-3'  > pfam/sequence_results_red_e_cutoff.txt

#get best hit based on bit score, and then e-value
sort -t$'\t' -k3,3gr -k4,4g pfam/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  pfam/All_pfam.txt

#clean pfam id
awk -F'\t' -v OFS='\t' '{split($2,a,"."); print  $1,a[1], $3, $4}' pfam/All_pfam.txt > pfam/temp

#merge with protein list
LC_ALL=C join -e'-' -a1 -t $'\t' -j1 -o 0,2.2,2.3,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort pfam/temp) | sort -u -k1 > pfam/temp3

#merge with pfam names
LC_ALL=C join -a1 -1 2 -2 1 -t $'\t' -e'-' -o1.1,0,2.5,1.4 <(LC_ALL=C sort -k2 pfam/temp3) <(LC_ALL=C sort -k1 $pfam_mapping) | LC_ALL=C sort > pfam/temp5

#add in header
echo -e "accession\tPFAM_hmm\tPFAM_description\tPfam_Evalue" | cat - pfam/temp5 > pfam/E_Pfam.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' pfam/E_Pfam.tsv merge/temp_BCM.tsv > merge/temp_BCME.tsv

#control lines
wc -l merge/*

#clean up
rm pfam/temp*
```






## TIGR search
##################################

The original TIGRFAMs database was a research project of The Institute for Genomic Research (TIGR) and its successor, the J. Craig Venter Institute (JCVI) and are now maintained by NCBI. TIGRFAMs is a collection of manually curated protein families focusing primarily on prokaryotic sequences.

Possible settings to change:

- number of CPUs
- e-value threshold


```{bash class.source = 'fold-show', eval=FALSE}
#1. Setup directories
mkdir TIGRs

#2. run hmmsearch against all tigrs
hmmsearch --tblout TIGRs/sequence_results.txt -o TIGRs/results_all.txt --domtblout TIGRs/domain_results.txt --notextw --cpu $cpus $tigr_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select sequences above a certain e-value
sed 's/ \+ /\t/g' TIGRs/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $4, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-3'  > TIGRs/sequence_results_red_e_cutoff.txt

#get best hit based on bit score, and then e-value
sort -t$'\t' -k3,3gr -k4,4g TIGRs/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  TIGRs/All_TIGR.txt

#merge with protein list
LC_ALL=C join -a1 -e'-' -t $'\t' -j1 -o 0,1.2,2.2,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort TIGRs/All_TIGR.txt) | sort -u -k1 > TIGRs/temp

#merge with tigr names
LC_ALL=C join -a1 -1 3 -2 1 -e'-' -t $'\t' -o1.1,0,2.11,2.13,1.4  <(LC_ALL=C sort -k3 TIGRs/temp) <(LC_ALL=C sort -k1 $tigr_mapping)| LC_ALL=C  sort > TIGRs/temp3

#add header
echo -e "accession\tTIRGR\tTIGR_description\tEC\tTIGR_Evalue" | cat - TIGRs/temp3 > TIGRs/F_TIGR.tsv

#combine with previous data frame
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' TIGRs/F_TIGR.tsv merge/temp_BCME.tsv > merge/temp_BCMEF.tsv

#control lines
wc -l merge/*

#clean up
rm TIGRs/temp*
```






## CazyDB search
##################################

The CAZy (Carbohydrate-Active enZYmes) database describes the families of structurally-related catalytic and carbohydrate-binding modules (or functional domains) of enzymes that degrade, modify, or create glycosidic bonds. As such, this database does not consist of clusters of orthologous groups but individual protein sequences that can be linked to CAZys.

Possible settings to change:

- number of CPUs
- e-value threshold (default = 1e-10). We are a bit stricter here, because the database is smaller --> so the chance of false-positives is higher.


```{bash class.source = 'fold-show', eval=FALSE}
#make folder
mkdir CAZYmes

#2. run hmmsearch against the CAZy database
hmmsearch --tblout CAZYmes/sequence_results.txt -o CAZYmes/results_all.txt --domtblout CAZYmes/domain_results.txt --notextw --cpu $cpus $cazy_db faa/Prokka/All_Genomes_clean.faa

#format the full table and only select sequences above a certain e-value
sed 's/ \+ /\t/g' CAZYmes/sequence_results.txt | sed '/^#/d'| sed 's/ /\t/g'| awk -F'\t' -v OFS='\t' '{print $1, $3, $6, $5}' | awk -F'\t' -v OFS='\t' '($4 + 0) <= 1E-10'  > CAZYmes/sequence_results_red_e_cutoff.txt

#get best hit based on bit score, and then e-value
sort -t$'\t' -k3,3gr -k4,4g CAZYmes/sequence_results_red_e_cutoff.txt | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k3,3gr -k4,4g >  CAZYmes/All_vs_CAZYmes.txt

#clean things and remove .hmm from the file
sed s'/\.hmm//g'  CAZYmes/All_vs_CAZYmes.txt > CAZYmes/temp

#merge with contig list
LC_ALL=C join -a1  -j1 -e'-' -t $'\t' -o 0,2.2,2.4 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort CAZYmes/temp)  | LC_ALL=C sort  > CAZYmes/temp3

#merge with mapping file
LC_ALL=C join -a1  -1 2 -2 1 -e'-' -t $'\t' -o 1.1,1.2,1.3,2.2 <(LC_ALL=C sort -k2 CAZYmes/temp3) <(LC_ALL=C sort $cazy_mapping) | LC_ALL=C sort  > CAZYmes/temp4

#add in headers
echo -e "accession\tCAZy\tCAZy_evalue\tDescription" | cat - CAZYmes/temp4 > CAZYmes/G_CAZy.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' CAZYmes/G_CAZy.tsv merge/temp_BCMEF.tsv > merge/temp_BCMEFG.tsv

#control lines
wc -l merge/*

#clean up
rm CAZYmes/temp*
```




## Transporter DB search
##################################

The transporter classification database (also tcdb)  details a comprehensive classification system for membrane transport proteins known as the Transporter Classification (TC) system. The database consists of individual protein sequences linked to this classification system.

Possible settings to change:

- number of CPUs (default setting for -num_threads = 2)
- e-value threshold (default setting for -evalue = 1e-10)  We are a bit stricter here, because the database is smaller --> so the chance of false-positives is higher.


```{bash class.source = 'fold-show', eval=FALSE}
#make folder
mkdir TransporterDB

#2. run protein blast against all TransporterDB
blastp -num_threads $cpus -outfmt 6 -query faa/Prokka/All_Genomes_clean.faa -db $transporter_db -out TransporterDB/All_vs_TPDB.tsv -evalue 1e-10

#get best hit based on bit score, and then e-value
sort -t$'\t' -k12,12gr -k11,11g TransporterDB/All_vs_TPDB.tsv | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k12,12gr -k11,11g >  TransporterDB/temp

#merge with contig list
LC_ALL=C join -a1  -j1 -e'-' -t $'\t' -o 0,2.2,2.11 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort TransporterDB/temp)  | LC_ALL=C sort  > TransporterDB/temp3

#merge with mapping file
LC_ALL=C join -a1  -1 2 -2 1 -e'-' -t $'\t' -o 1.1,1.2,1.3,2.3 <(LC_ALL=C sort -k2 TransporterDB/temp3) <(LC_ALL=C sort -k1 $transporter_mapping) | LC_ALL=C sort  > TransporterDB/temp4

#add in headers
echo -e "accession\tTCBD_ID\tTCBD_evalue\tTCBD_description" | cat - TransporterDB/temp4 > TransporterDB/I_TPDB.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' TransporterDB/I_TPDB.tsv merge/temp_BCMEFG.tsv > merge/temp_BCMEFGHI.tsv

#control lines
wc -l merge/*

#clean up
rm TransporterDB/temp*
```





## HydDB search 
##################################

The hydrogenase database (HydDB) contains individual protein sequences from hydrogenases and is a curated database of hydrogenases by known type.

Possible settings to change:

- number of CPUs (default setting for -num_threads = 2)
- e-value threshold (default setting for -evalue = 1e-10)  We are a bit stricter here, because the database is smaller --> so the chance of false-positives is higher.


```{bash class.source = 'fold-show', eval=FALSE}
#make folder
mkdir HydDB

#run search
blastp -num_threads $cpus -outfmt 6 -query faa/Prokka/All_Genomes_clean.faa -db $hyd_db -out HydDB/All_vs_HydDB.tsv -evalue 1e-10

#get best hit based on bit score, and then e-value
sort -t$'\t' -k12,12gr -k11,11g HydDB/All_vs_HydDB.tsv | sort -t$'\t' --stable -u -k1,1  | sort -t$'\t' -k12,12gr -k11,11g >  HydDB/temp

#merge with contig list; control: grep for 05160
LC_ALL=C join -a1  -j1 -e'-'  -t $'\t' -o 0,2.2,2.3,2.11,2.12 <(LC_ALL=C sort FileLists/1_Bins_to_protein_list.txt) <(LC_ALL=C sort HydDB/temp) | LC_ALL=C sort  > HydDB/temp2

#merge with HydDB names
LC_ALL=C join -a1 -1 2 -2 1 -e'-' -t $'\t' -o1.1,0,2.2,1.4  <(LC_ALL=C sort -k2  HydDB/temp2) <(LC_ALL=C sort -k1 $hyd_mapping)  | LC_ALL=C  sort > HydDB/temp3

#add in headers
echo -e "accession\tHydDB\tDescription\tHydDB_evalue" | cat - HydDB/temp3 > HydDB/J_HydDB.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' HydDB/J_HydDB.tsv merge/temp_BCMEFGHI.tsv > merge/temp_BCMEFGHIJ.tsv

#control lines
wc -l merge/*

#clean up
rm HydDB/temp*
```



## IPRscan --> done
##################################

InterPro provides functional analysis of proteins by classifying them into families and predicting domains and important sites, such as transmembrane motives (TMs) or signal proteins.

<div class="alert-danger">
Both IPRscan and the Diamonds search against NCBI_nr will take quite some time. While you are waiting for this: <br/><br/>

Continue with the step **Modify final dataframe** and instead of using merge/temp_BCMEFGHIJKL.tsv  use the file merge/temp_BCMEFGHIJ.tsv in the first step. Follow through this section, download the table and have a look at what you got. 

</div>

### genome cleanup

IPRscan does not like to deal with ``*`` in the protein files, which is used for stop codons. If we work with programs that don't recognize ``*``, we should remove this symbol first.

```{bash class.source = 'fold-show', eval=FALSE} 
#interproscan does not like if we have stars in our protein files, so we need to remove them
cd faa/Prokka
sed 's/*//g' All_Genomes_clean.faa >  All_Genomes_clean_interproscan.faa 
cd ../..
```

Before running IPRscan, we can split our big concatenated file into several smaller files.
Then we can use parallel GNU to speed things up. 

We set the number of processes to run with ``-j``. Always make sure that you adjust this depending on the avail resources. 

### run iprscan

Possible settings to change:

- Number of files we split the concatenated fasta file into via Split_Multifasta.py (default of -n = each split file will contain 1000 sequences)
- number of parallel jobs to run via parallel (default setting for -j= 4)
- e-value threshold (default setting for -evalue = 1e-20)


```{bash class.source = 'fold-show', eval=FALSE}

mkdir IPRscan

#split files --> 2069 proteins were split into 3 files
python $python_scripts/Split_Multifasta.py -m faa/Prokka/All_Genomes_clean_interproscan.faa -n 1000

#move files into a better location
mkdir split_faa
mv File*.faa split_faa

#iprscan: (interproscan-5.31-70.0)
#"cite:Tange (2011): GNU Parallel - The Command-Line Power Tool, ;login: The USENIX Magazine, February 2011:42-47."

parallel -j3 'i={}; interproscan.sh -i $i -d IPRscan/ -appl TIGRFAM,SFLD,SUPERFAMILY,Gene3D,Hamap,Coils,ProSiteProfiles,SMART,CDD,PRINTS,ProSitePatterns,Pfam,ProDom,MobiDBLite,PIRSF,TMHMM -T IPRscan/temp --iprlookup --goterms -cpu 5' ::: split_faa/File*.faa

#if parallel is not avail on your system, use a for-loop
#for sample in split_faa/File*.faa; do  interproscan.sh -i $sample -d IPRscan/ -appl TIGRFAM,SFLD,SUPERFAMILY,Gene3D,Hamap,Coils,ProSiteProfiles,SMART,CDD,PRINTS,ProSitePatterns,Pfam,ProDom,MobiDBLite,PIRSF,TMHMM -T IPRscan/temp --iprlookup --goterms -cpu 5; done

#Concat result files
cd IPRscan
mkdir single_files
mv File* single_files
mkdir Concat_results/

cat single_files/File*.faa.xml > Concat_results/All_bins_iprscan-results.xml
cat single_files/File*.faa.gff3 > Concat_results/All_bins_bins_iprscan-results.gff3
cat single_files/File*.faa.tsv > Concat_results/All_bins_bins_iprscan-results.tsv

#cleanup
rm single_files/File*

cd ..

#clean up fasta header so that it is exactly the same as the accession ID in the interproscan results
python $python_scripts/parse_IPRdomains_vs2_GO_2_ts_sigP.py -s faa/Prokka/All_Genomes_clean.faa -i IPRscan/Concat_results/All_bins_bins_iprscan-results.tsv -o IPRscan/All_bins_bins_iprscan-results_parsed.tsv

#remove issue with spaces
sed 's/ /_/g' IPRscan/All_bins_bins_iprscan-results_parsed.tsv | LC_ALL=C  sort > IPRscan/temp.txt

#print only columns of interest
awk -F'\t' -v OFS="\t"  '{print $1, $2,$3,$4, $5, $8, $9}' IPRscan/temp.txt > IPRscan/temp2

#add header
echo -e "accession\tIPR\tIPRdescription\tIPR_PFAM\tIPR_PFAMdescription\tTMHMM\tSignalP" | cat - IPRscan/temp2 > IPRscan/K_IPR.tsv

#combine with previous data frame
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' IPRscan/K_IPR.tsv merge/temp_BCMEFGHIJ.tsv > merge/temp_BCMEFGHIJK.tsv

#control lines
wc -l merge/*

#clean up
rm IPRscan/temp*

```




## Diamond against NCBI NR
##################################

Possible settings to change:

- number of CPUs (--threads)
- e-value threshold (--evalue)
- Sequence numbers to display (--seq)

Notice:

- This step can take quite some time.
- the NCBI nr database is constantly updated (and very large), which is why we are not providing the database with the other files. For information on how to download ncbi_nr check out the information [here](https://www.ncbi.nlm.nih.gov/books/NBK569856/).


```{bash class.source = 'fold-show', eval=FALSE}
#make folder
mkdir NCBI_NR

#run diamond against NR database
diamond blastp -q faa/Prokka/All_Genomes_clean.faa --more-sensitive --evalue 1e-3 --threads $cpus --seq 50 --db $ncbi_nr_db --taxonmap $taxon_db --outfmt 6 qseqid qtitle qlen sseqid salltitles slen qstart qend sstart send evalue bitscore length pident staxids -o NCBI_NR/All_NCBInr.tsv

#Select columns of interest in diamond output file
awk -F'\t' -v OFS="\t" '{ print $1, $5, $6, $11, $12, $14, $15 }'  NCBI_NR/All_NCBInr.tsv | sed 's/ /_/g' > NCBI_NR/temp

#Parse Diamnond Results
python $python_scripts/parse_diamond_blast_results_id_taxid.py -i FileLists/AllProteins_list.txt -d NCBI_NR/temp -o NCBI_NR/temp2

#rm header
sed 1d NCBI_NR/temp2 > NCBI_NR/temp3

#add an '-' into empty columns
awk -F"\t" '{for(i=2;i<=NF;i+=2)gsub(/[[:blank:]]/,"_",$i)}1' OFS="\t" NCBI_NR/temp3 > NCBI_NR/temp4

#the python script above sometimes leaves an empty 7th column, this gets rid of that issue
awk -F'\t' -v OFS="\t"  '{if (!$7) {print $1,$2, $4 , $6, "-"} else {print $1, $2, $4, $6, $7}}' NCBI_NR/temp4 | LC_ALL=C sort > NCBI_NR/temp5

#split columns with two tax ids
awk -F'\t' -v OFS='\t' '{split($5,a,";"); print $1, $2, $3, $4, a[1]}' NCBI_NR/temp5 > NCBI_NR/temp6

#in column 2 remove everything after < (otherwise the name can get too long)
awk -F'\t' -v OFS='\t' '{split($2,a,"<"); print $1, a[1], $3, $4, $5}' NCBI_NR/temp6 > NCBI_NR/temp6a

#merge with tax names
LC_ALL=C join -a1 -1 5 -2 1 -e'-' -t $'\t'  -o1.1,1.2,1.3,1.4,1.5,2.2  <(LC_ALL=C sort -k5  NCBI_NR/temp6a) <(LC_ALL=C sort -k1 ~/../spang_team/Databases/NCBI_taxonomy/Jul2019/taxonomy5.txt ) | LC_ALL=C  sort > NCBI_NR/temp7

#add in header
echo -e "accession\tTopHit\tE_value\tPecID\tTaxID\tTaxString" | cat - NCBI_NR/temp7 > NCBI_NR/L_Diamond.tsv

#combine with previous dataframe
awk 'BEGIN{FS="\t";OFS="\t"}FNR==NR{a[$1]=$0;next}{print $0,a[$1]}' NCBI_NR/L_Diamond.tsv merge/temp_BCMEFGHIJK.tsv > merge/temp_BCMEFGHIJKL.tsv

#control lines
wc -l merge/*

#clean up
rm NCBI_NR/temp*
```




###################################################################################
# Sidenote
###################################################################################

The step-by-step addition of the individual dataframes can get rather tedious esp. if you need to run things over again.
We will not do this in this tutorial but there are other means to merge several files with a common column, which you might want to explore and for which we could write a small python script later on.


```python

from functools import reduce
reduce(lambda x,y: pd.merge(x,y, on='Col1', how='outer'), [df1, df2, df3])

```



###################################################################################
###################################################################################
# Modify final dataframe
###################################################################################
###################################################################################

Now, lets cleanup the final dataframe and then we made it!

Basically, we do several cosmetic steps. I.e. for each database search we get a column with the accessionID, in the second step, we remove the duplicated columns to slim the data down.


```{bash class.source = 'fold-show', eval=FALSE}
#duplicate first column and add new header for this column
awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$0}' merge/temp_BCMEFGHIJKL.tsv | sed '1s/accession/for_merge/' > merge/temp_0.tsv

#rm redundant headers and reconstruct accession ID
awk 'BEGIN{FS="\t";OFS="\t"}NR==1{for(x=1;x<=NF;x++)if($x!="accession")l[x]++;}{for(i=1;i<=NF;i++)if(i in l)printf (i==NF)?$i"":$i"\t";printf "\n"}' merge/temp_0.tsv > merge/temp_1.tsv

#change column name
awk 'BEGIN{FS="\t";FS="\t"; OFS="\t"}{if(NR==1) $1="accession"} {print $0 }' merge/temp_1.tsv > merge/Annotations.txt

#control first 50 columns
head -50 merge/Annotations.txt > merge/temp_2.tsv

#clean up
rm merge/temp_[0-9].tsv
gzip merge/*tsv
```



**Question/Homework --> Let's figure out what Igniccous is doing**

- Download the table `Annotations.txt` to your personal computer using scp.
- Compare the different annotations, have a look at how they differ.
- Make one separate table for each of the two genomes
- For Ingicoccus (GCA_000017945), copy the accession numbers and KO ids into a new file and upload or copy and paste them to the [KEGG mapper](https://www.genome.jp/kegg/mapper.html). Use this tool to get some first insights into the metabolic capacities of that organism. 

Based on this try to answer:

  - **Does Igni use glycolysis or gluconeogenesis (or both)?** <br/>
  This question allows us to figure out whether our organism might use carbon sources like sugars (via glycolysis, if we have this we might look further into whether we have a fermenting organism) or uses acetyl-CoA to make their own sugars, which could be then used in other pathways (gluconeogenesis). <br/>

  - **Is Igni likely an anaerobic or aerobic organism?** <br/>
  This question can help us to decide what growth conditions we need and whether O2 is required for our organism to generate energy.

  - **Is Ingi and autotroph or heterotroph?** <br/>
  Again, this question is important if we want to know how to grow our organism, do we need to feed it sugars so that it grows, or can it grow on CO2. If it is involved in CO2 and/or methane cycling it might also have an impact in global carbon cycles.

  - **Can Igni reduce sulfur?**  <br/>
  If Igni grows on sulfur, we could add sulfur to a potential growth medium. Plus, for the environmental side of things we want to look into lithotrophy and how organism in our environment are involved in nitrogen, sulfur, metal or hydrogen cycling.

  - **If Igni reduces sulfur, what does it oxidize?** <br/>
  Same as above.

The Input for the website should look something like this:
```{r fig3, echo=FALSE, fig.cap="", out.width = '80%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations3.png")
```



Answer for **Does Igni use glycolysis or gluconeogenesis?**

<div id="spoiler1" style="display:none">

Ignicoccus has a gluconegenic pathway. The key enzyme to distinguish glycolysis from gluconeogenesis is [EC:3.1.3.11](https://biocyc.org/gene?orgid=META&id=F16B-MONOMER) which allows this organism to go from acetate to glucose.
To figure this out, next to adding a literature search, there are two routes to take:<br/>
Option A:
```{r fig4A2, echo=FALSE, fig.cap="", out.width = '100%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/KeegMapper1.png")
```
Option B:
```{r fig4b2, echo=FALSE, fig.cap="", out.width = '100%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/KeegMapper2.png")
```

</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoiler1') .style.display=='none')
              {document.getElementById('spoiler1') .style.display=''}
            else{document.getElementById('spoiler1') .style.display='none'}">
  Show/hide
</button>
<br/><br/>


Answer for **Is Igni likely an anaerobic or aerobic organism?**

<div id="spoiler2" style="display:none">

There are several things to check, but oxidative phosphorylation is a good start. We can see that complex IV, the terminal oxidase of the respiratory chain that transfers electrons from cyt c to O2, is absent. So this organism is likely an anaerobe
```{r fig5, echo=FALSE, fig.cap="", out.width = '80%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations5.png")
```

</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoiler2') .style.display=='none')
              {document.getElementById('spoiler2') .style.display=''}
            else{document.getElementById('spoiler2') .style.display='none'}">
  Show/hide
</button>
<br/><br/>

Answer for **Is Ingi and autotroph or heterotroph?**


<div id="spoiler3" style="display:none">

Igni is an autotroph that encodes the dicarboxylate-hydroxybutyrate cycle (big red circle highlights the genes of the pathway.<br/><br/>
This one is a bit more tricky, since the malate dehydrogenase, 1.1.1.37, is absent if we look at the map (small red box).<br/>
However, if we dig into the annotation tables, we can find the blast hit to WP_012123403.1 (likely gene ID 13020) hits to a malate dehydrogenase and we can see that most databases view the gene as a malate/lactate dehydrogenase (i.e. arCOG00246). Malate dehydrogenases are very easy to mix-up with lactate dehydrogenases and that happened here with the KO annotations. 
<br/><br/>

As a general tip, if only a few genes of a pathway are lacking --> have a second look into the table and search the other annotations. If we have genes where its hard to distinguish between substrate specificity sometimes a tree might help or looking at the catalytic site of a protein, depending on how well studied our protein of interest is
```{r fig6, echo=FALSE, fig.cap="", out.width = '80%', fig.align="left"}
knitr::include_graphics("/Users/ninadombrowski/Desktop/WorkingDir/Notebooks/Pictures/Annotations6.png")
```

</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoiler3') .style.display=='none')
              {document.getElementById('spoiler3') .style.display=''}
            else{document.getElementById('spoiler3') .style.display='none'}">
  Show/hide
</button>
<br/><br/>

Answer for **Can Igni reduce sulfur? **

<div id="spoiler4" style="display:none">

Yes, Igni can reduce sulfur using K18367/arCOG01069, which encodes CoA-dependent NAD(P)H_sulfur_oxidoreductase [EC:1.8.1.18]. <br/><br/>
This one is a bit tricky and something you do not find this hit looking at the pathway maps alone. 
<br/>To find the issue with the maps: go to the metabolic pathway overview, click on show matched objects and search for K18367. You will see we will not find a hit here, but if we search for this KO in our table we can find it, so what is the issue? <br/><br/>If we search for K18367 across the KEGG database, we find some info [here](https://www.genome.jp/dbget-bin/www_bget?ko:K18367). We first see that this enzyme is not linked to a pathway (and a problem of the pathway map tool) and also the explanation that this enzyme is involved in NAD(P)H-linked sulfur reduction. We can confirm this by showing that (a) this gene is similar to a characterized gene in Pyrococcus but also since we find sulfite transporters (i.e. arCOG02050/WP_011998105.1) which might be involved in getting the produced H2S outside of the cell.
</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoiler4') .style.display=='none')
              {document.getElementById('spoiler4') .style.display=''}
            else{document.getElementById('spoiler4') .style.display='none'}">
  Show/hide
</button>
<br/><br/>


Answer for **If Igni reduces sulfur, what does it oxidize? **

<div id="spoiler5" style="display:none">

Igni gains energy by reduction of elemental sulfur to hydrogen sulfide using molecular hydrogen as the electron donor. <br/><br/>Again, this is not something the Kegg mapper will tell us. To figure this out we can again use different approaches: (a) search the literature and find out what related organisms do, i.e. Pyrococcus and see whether our organism has a similar gene set, (b) search our Annotation table for potential electron donors that might make sense and are also used by other S reducers. 
<br/><br/>The second approach (and some background reading in most cases) will lead us to K06281 and K06282. These hit to a [NiFe]_Group_1g hydrogenase and we can check [here](https://services.birc.au.dk/hyddb/browser/class/nife-group-1g/) what exactly this hydrogenase typically is doing. We need to check this, because for hydrogenases it is sometimes tricky to say, whether they are involved energy reactions. If we check this specific example we find that it is involved in hydrogenotrophic respiration using sulfur as terminal electron acceptor, so perfect!
<br/><br/>A word of warning, the HydDB for hydrogenases will never be perfect. So also here, for more uncertain hits it might make sense to build a tree for hydrogenases to see if we really have a 1g hydrogenase.
</div>

<button title="Click to show answer" type="button"
   onclick="if(document.getElementById('spoiler5') .style.display=='none')
              {document.getElementById('spoiler5') .style.display=''}
            else{document.getElementById('spoiler5') .style.display='none'}">
  Show/hide
</button>
<br/><br/>



###################################################################################
###################################################################################
# What next?
###################################################################################
###################################################################################

Depending on the number of genomes we have the final table will get rather large. Therefore, it makes a lot of sense to condense the information. The following ideas will not be part of this tutorial but will give you some pointers on how do to that:

## R magic

We provide a R tutorial [here](https://github.com/ndombrowski/R_tutorial). This tutorial starts with basics in R but ends with an example to deal with the annotation tables here. More specifically, it will show you how to:

- Make count tables for each genome, i.e. for MAG1 count how often we find arCOG0001. This calculation is done each genome that is part of your analysis and thus provides you with more condensed table
-  Condense the data on the same on taxonomic or custom levels. I.e. the example in the R tutorial is for ~30 genomes of the same phylum. However, the genomes come from different environments. Using the tutorial you will be able to ask, how many of the aquifer genomes have arCOG0001? How many of the marine genomes have this gene?

## Python magic

This section will be extended once the author of this tutorial has time (o.O). However, once you have to deal with 5000 genomes python becomes much more attractive to do the parsing esp. compared to RStudio. 





###################################################################################
###################################################################################
# Background info on the used tools and databases
###################################################################################
###################################################################################

Here, you will some background info on the used approaches, databases and tools.
The idea is to provide links to the manuals and some other information that might (or might not) be useful.

Note, this is an on-going part of this document, so it might change over time.


## Gene predictions

There are several tools to identify genes in our genomes

- [Prodigal](https://github.com/hyattpd/Prodigal): a prokaryotic gene recognition and translation initiation site identification tool. [Hyatt et al., 2010](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2848648/)
- [Prokka](https://github.com/tseemann/prokka): A rapid prokaryotic genome annotation tool. [Seeman 2014](https://pubmed.ncbi.nlm.nih.gov/24642063/) This tool is based on prodigal but does some basic genome annotations.

Something important to keep in mind is that these gene prediction tools generally used for genomes from archaea and bacteria. 
If you are interested in eukaryotic genomes, things like splicing and alternative reading frames make gene prediction more tricky and there are special tools out there that are more suited for Eukaryotes.


## Information on the different search tools
##################################

### BLAST

BLAST (Basic Local Alignment Search Tool) finds regions of local similarity between protein or nucleotide sequences [(Altschul 1990)](https://pubmed.ncbi.nlm.nih.gov/2231712/).

NCBI provides a [quick tutorial](https://www.ncbi.nlm.nih.gov/books/NBK1734/) explaining its general usage and you can find some more details in the algorithm and usage of the website [here](https://www.ccg.unam.mx/~vinuesa/tlem/pdfs/Bioinformatics_explained_BLAST.pdf).

There are different versions of blast that you can use depending on what comparisons you want to make.

- blastn = nucleotide blast searches with a nucleotide “query” against a database of nucleotide “subject” sequences.
- blastp = protein blast searches with a protein “query” against a database of protein “subject” sequences.
- translated blasts searches use a genetic code to translate either the “query,” database “subjects,” or both, into protein sequences, which are then aligned as in “blastp.”
  - tblastn = a protein sequence query is compared to the six-frame translations of the sequences in a nucleotide database.
  - blastx = a nucleotide sequence query is translated in six reading frames, and the resulting six-protein sequences are compared, in turn, to those in a protein sequence database.
  - tblastx = both the “query” and database “subject” nucleotide sequences are translated in six reading frames, after which 36 (6 × 6) protein “blastp” comparisons are made.

Additionally, we can also work with alignments using psi-blast = protein sequence profile search method that builds off the alignments generated by a run of the BLASTp program. Some more info into this algorithm can be found [here](https://www.ncbi.nlm.nih.gov/books/NBK2590/#:~:text=PSI-BLAST%20(Position-Specific,threshold%20using%20protein–protein%20BLAST.)



### Diamond

Diamond is an alternative to BLAST for protein and translated DNA searches [(Buchfink et al., 2015)](https://www.nature.com/articles/nmeth.3176) If you check the diamond [website](http://www.diamondsearch.org/index.php) you will find that some advantages over blast are listed:

- Pairwise alignment of proteins and translated DNA at 100x-20,000x speed of BLAST.
- Frameshift alignments for long read analysis.
- Low resource requirements and suitable for running on standard desktops or laptops.
- Various output formats, including BLAST pairwise, tabular and XML, as well as taxonomic classification.

In principle diamond works very similar compared to blast and uses the same input files.


### HMMER

We also use HMMER for searching databases and generating alignments (but generally we prefer to use mafft-linsi for alignments ourselves).

The approach taken by HMMER is quite different from blast and diamond as HMMER implements methods using probabilistic models called profile hidden Markov models (profile HMMs).

If you look at some HMM databases used in this tutorial with ie. `less arcog.hmm`, you will see that the files look very different compared to what we use as input for blast and diamond. One of the biggest differences is that the profiles used in HMMER searches are build from multiple sequence alignments.

For more information you can check the [HMMER website](http://hmmer.org) or the [hmmer paper](https://pubmed.ncbi.nlm.nih.gov/20180275/)


## Information on the different databases
##################################

## COGs

COG = Clusters of Orthologous Groups.

The COG protein database was generated by comparing predicted and known proteins in all completely sequenced microbial genomes to infer sets of orthologs. Each COG consists of a group of proteins found to be orthologous across at least three lineages.

The different COGs are characterized into different pathways and the pathway description can be found [here](http://clovr.org/docs/clusters-of-orthologous-groups-cogs/)

Some papers you might want to check if you want to know more:

- [The NCBI website](https://www.ncbi.nlm.nih.gov/research/cog-project/)
- [The COG approach](https://pubmed.ncbi.nlm.nih.gov/28968633/)
- [The Original COG paper ](https://pubmed.ncbi.nlm.nih.gov/9381173/)


## arCOGs

arCOG = Archaeal Clusters of Orthologous Genes.
arCOGs are generated in a similar manner as the COGs, just with a focus on archaeal genomes.

Some papers you might want to check if you want to know more:

- [The NCBI website](https://www.ncbi.nlm.nih.gov/research/cog-project/)
- [The first description of the arCOGs](https://pubmed.ncbi.nlm.nih.gov/25764277/)


## KOs and KEGG

KO = KEGG ORTHOLOGY

The KO (KEGG Orthology) database is a database of molecular functions represented in terms of functional orthologs. By using the HMMs we assign K numbers to our proteins by HMMER/HMMSEARCH against the KOfam database.

Details on how the HMMs and cutoffs were generated can be found [here](https://academic.oup.com/bioinformatics/article/36/7/2251/5631907)

Something useful to know is that all KOs, ie. K00161, are linked to metabolic maps that have quite good descriptions. For example, if you would google K00161 you will come across this [website](https://www.genome.jp/dbget-bin/www_bget?ko:K00161) that gives:

- a detailed description of what gene we have
- if applicable the enzyme number (EC number)
- The pathways this gene can be involved in. If you click on the pathway you see its exact position in a pathway
- A paper that first describes the gene

If you click on the enzyme number (EC 1.2.4.1) you will get other useful information

- alternative names
- the "chemical" class it belongs to
- the chemical reaction 
- useful comments about the reaction and whether it is part of an enzyme complex
- if it is an enzyme complex, we get the KOs for the other subunits


## Biocyc: Alternative to KEGG to see how a gene is linked to a pathway

If you want to find out more about a gene, you can check also metacyc.
I.e. if we search for ``"pyruvate dehydrogenase" biocyc``in goggle we should land [here](https://biocyc.org/META/NEW-IMAGE?type=ENZYME&object=CPLX-2022)

Similar as with KEGG, we get some basic information about this enzyme, such as:

- the different subunits
- a basic description of what that enzyme does
- the reaction, incl. the direction
- Functions of the other subunits
- some paper references


## PFAM

[PFAM](https://pfam.xfam.org) is a database is a large collection of protein families, each represented by multiple sequence alignments and hidden Markov models (HMMs).


## TIGRFAMs

The [TIGRFAMs](http://tigrfams.jcvi.org/cgi-bin/index.cgi) database is a resource consisting of curated multiple sequence alignments, Hidden Markov Models (HMMs) for protein sequence classification, and associated information designed to support automated annotation of (mostly prokaryotic) proteins.

## CAZy

CAZy = Carbohydrate-Active enZYmes

A  database dedicated to the display and analysis of genomic, structural and biochemical information on Carbohydrate-Active Enzymes (CAZymes).

Resources:

- [Here](http://www.cazy.org) we have a database with more information on indiv. CAZymes and potential substrates
- [dbCAN](dbCAN2 meta server) is a webserver for CAZYme annotation. We also get the sequences we use for our workflow from this website.
- The [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2686590/) describing the CAZY database


## Transporter DB

[This](http://www.tcdb.org) is the website that gives more information about potential transporters found in our search.


## HydDB

HydDB = Hydrogenase database

Resources:

- [The original paper](https://www.nature.com/articles/srep34212)
- The [website](https://services.birc.au.dk/hyddb/) we originally got the sequences from. Also provides more detailed information on the different hydrogenase subcategories.


## IPRscan

Short for InterProScan

Interpro allows for the classification of protein families and predicting domains and important sites.

If we see an ID like this ``IPR000013``, this is a so called interpro domain that might be found on your protein. Notice, a protein can have more than one domain. 

Resources:

- The [website](https://www.ebi.ac.uk/interpro/search/sequence/) from which to download all information
- [The original paper](https://academic.oup.com/nar/article/47/D1/D351/5162469), which provides also an overview of the different databases integrated with InterPro.


## NCBI-nr

As we have seen above, NCBI (the National Center for Biotechnology Information) provides a lot of tools, such as blast, and databases, such as the COGs and arCOGs. 
Additionally, it is one of the largest repositories of sequences from metagenomes, amplicon libraries and genomes. 

The nr database (here: ncbi_nr) is compiled by the NCBI as a protein database for Blast searches.It contains non-identical sequences from GenBank CDS translations, PDB, Swiss-Prot, PIR, and PRF and is frequently updated. However, it is a quite large database making searches more time consuming.



## Uniprot - checking your results

Another good database, next to KEGG and metacyc, to check proteins is [Uniprot](https://www.uniprot.org).
The mission of UniProt is to provide the scientific community with a comprehensive, high-quality and freely accessible resource of protein sequence and functional information.

We can search for:

- protein names, i.e. "pyruvate dehydrogenase"
- proteins from certain taxa using the taxonomy option, i.e. ``taxonomy:bacteria pyruvate dehydrogenase``
- KOs, arCOGs, etc ...

I.e. if we would search for `` taxonomy:bacteria "pyruvate dehydrogenase e1" `` we should get:

- A list of all proteins classified as E1 subunit of the pyruvate dehydrogenase in bacteria
- some of these sequences are reviewed (i.e. might have some more supporting info) and others are unreviewed (annotation only based on sequence similarity)
- a list of taxa that have this protein (under view by click on taxonomy) and then we can see in exactly what bacteria these enzyme is found

If we click on an individual example, i.e. from E.coli, we get:

- Details on the reaction
- Publications about this specific enzyme in E.coli
- Co-factors
- Information about potential functions and biological processes
- a reference to MetaCyc
- The links to several databases (i.e. protein databases, KEGG, IPR domains, etc.)