README.Rmd

---
# title: "gifrop"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
[![install with bioconda](https://img.shields.io/badge/install%20with-bioconda-brightgreen.svg?style=flat)](http://bioconda.github.io/recipes/gifrop/README.html)  

# gifrop  

**G**enomic **I**slands **f**rom **Ro**ary **P**angenomes.  

This program is supposed to identify 'genomic islands' from roary pangenomes. The intent was to try and identify units of horizontal gene transfer within very closely related strains using a pangenome framework.  

In general, this type of analysis works best for very closely related strains, for example, a collection of *Salmonella enterica serovar typhimurium*.  

In theory, this program will work on pangenomes of any size, though things get very messy with too many genomes or genomes that are too different from one another.

**Prerequisites:**  

1. A pangenome calculated by roary:  Any pangenome should work, even those that do not split paralogues (the -s option of roary).  

2. the *.gff files used by roary to create the pangenome.  I have only used those produced by prokka but others may work.  The increment between locus tags must be 1.  
  
If all you have are nucleotide fastas then you can use the included `pan_pipe` script to take care of the whole pipeline from annotation with prokka to pangenome calculation with roary to genomic island extraction with gifrop.  
  
More or less this is a wrapper for exploring roary results.  
  

## conda installation  
  
**create new conda environment**  
`conda create -n gifrop`  
  
**activate new environment**  
`conda activate gifrop`  
  
**install gifrop**  
`conda install -c conda-forge -c bioconda -c defaults gifrop`
  
or  

## manual installation  

Install the dependencies and make sure they are in your path:  
  
  1) [GNU parallel](https://www.gnu.org/software/parallel/)  
  2) [abricate](https://github.com/tseemann/abricate)  
    - You will also need these custom abricate databases  
      - https://github.com/Jtrachsel/megares_db_4_abricate  
      - https://github.com/Jtrachsel/viroseqs  
  3) [R 3.6.*](https://www.r-project.org/) and these R packages:  
    - 'dplyr','tidyr','readr','tibble','ggplot2','purrr','Biostrings','BSgenome', 'igraph', 'pheatmap'  
  4) [roary](https://sanger-pathogens.github.io/Roary/) (if using pan_pipe script)  
  5) [prokka](https://github.com/tseemann/prokka) (if using pan_pipe script)  
  6) make the scripts from this repo available in your path  
  

<br>  

## Running gifrop  

- Navigate to a directory containing prokka annotated gff files and the output of roary, specifically the `gene_presence_absence.csv` file.  
- run gifrop:  
  - `gifrop --find_islands --threads 8`  
- all outputs go to the `gifrop_out` directory.  

## pan_pipe script  

The included script `pan_pipe` is a complete pipeline for extracting genomic islands from a collection of nucleotide fastas.

These are the basic steps:  
1) Annotate with prokka  
2) Calculate pangenome with roary  
3) Extract genomic islands with gifrop

  
You are able to pass arguments to prokka, roary and gifrop using quoted strings of arguments, see `pan_pipe --help` for more information.  


<br>  
    

## How this program identifies pangenomic islands  

1) **remove the core genome from the pangenome.**  
    - Currently, this removes any gene that occurs in all isolates exactly 1 time.  This means that genes that only occur in some genomes and genes that occur twice or more in any genome are currently kept for consideration.   
2) **Identify strings of consecutive genes** (locus tag numbers)  
3) **Remove strings that contain fewer than the minimum number of genes.**  
    - These remaining strings of consecutive genes are 'genomic islands' 'pan-genomic islands?'  
4) **Classify pangenomic islands** by running abricate with the following databases:  
     1) [Megares2.0](https://megares.meglab.org/) database (bc it has metal tolerance)  
     2) [plasmid finder](https://cge.cbs.dtu.dk/services/PlasmidFinder/) (replicon genes)  
     3) [vfdb](http://www.mgc.ac.cn/VFs/main.htm) (virulence genes)  
     4) [ProphET](https://github.com/jaumlrc/ProphET) phage db  
5) **Cluster pangenomic islands**, this can get messy.  
    - make a graph where pangenomic islands are nodes and they are joined by edges when they share genes with other pangenomic islands.  The edges of the graph are weighted in proportion to the number of genes that are shared by the nodes
    - two levels of clustering:  
           - primary clustering: any islands that share any number of genes will be in the same primary cluster.  
           - secondary clustering: this applies louvain sp? clustering on the graph.  This should hopefully help to separate some of the more tangled islands.  
           - problems with clustering: in cases where there are two plasmids that are very different but have some shared genes, these will often be in the same cluster.  I have seen a case where the same transposon is inserted into different genomic islands and plasmids and then they get all tangled together.  ANyway, it gets messy so be careful what your inputs are.  
  6) **output pangenomic island info.**  
    All outputs go to the `gifrop_out` folder.  
      main outputs:  
      
           1. clustered_island_info.csv  
           2. pan_with_island_info.csv  
           3. pan_only_islands.csv  
           4. figures directory  
              - some basic exploratory figs  
              - heatmaps of the variable secondary clusters (buggy)  
           5. All the pangenomic islands  
              - both fasta format and gffs  
           6. logfiles  
  
  
<br>  
  
| file | description|  
|------|------------|  
|clustered_island_info.csv | This is a csv file that contains a detailed description of every genomic island detected in this pangenome. |  
|pan_with_island_info.csv | The roary 'gene_presence_absence.csv' file annotated with information about which genomic islands and clusters the genes belong to.|  
| pan_only_islands.csv | The same as the previous file but filtered to only include genes that occur on at least 1 genomic island|  

  
  

## graphical overview:  
 
 
```{r, echo=FALSE, message=FALSE}
set.seed(1)
library(tidyverse)


colours = c("core" = "#268bd2",
            "non-core" = "#b58900",
            "island" = "red",
            "non-island" = "grey")


### genome 1
tst <- tibble(locus_tag_num = seq(1,50),
              locus_tag = paste('genome1', locus_tag_num, sep = '_'),
       strt=seq(1,50),
       stp =seq(2,51),
       gene_type=sample(c('core', 'non-core'), 50, replace = TRUE),
       ymin=0,
       ymax=1,
       genome='genome1')
# tst

filt <- tst %>% filter(gene_type == 'non-core')

step4 <- split(filt$locus_tag_num, cumsum(c(1, diff(filt$locus_tag_num) != 1)))    # identifies streches of sequential loc.tags
#browser()
step5 <- step4[lapply(step4, length) >= 3]
islands <- unlist(step5)


tst$island <- ifelse(tst$locus_tag_num %in% islands, 'island', 'non-island')
genome1 <- tst
#### genome2 ###

tst <- tibble(locus_tag_num = seq(1,50),
              locus_tag = paste('genome2', locus_tag_num, sep = '_'),
              strt=seq(1,50),
              stp =seq(2,51),
              gene_type=sample(c('core', 'non-core'), 50, replace = TRUE),
              ymin=-2,
              ymax=-1,
              genome='genome2')
# tst

filt <- tst %>% filter(gene_type == 'non-core')

step4 <- split(filt$locus_tag_num, cumsum(c(1, diff(filt$locus_tag_num) != 1)))    # identifies streches of sequential loc.tags
#browser()
step5 <- step4[lapply(step4, length) >= 3]
islands <- unlist(step5)


tst$island <- ifelse(tst$locus_tag_num %in% islands, 'island', 'non-island')

genome2 <- tst

#### genome3 ###

tst <- tibble(locus_tag_num = seq(1,50),
              locus_tag = paste('genome3', locus_tag_num, sep = '_'),
              strt=seq(1,50),
              stp =seq(2,51),
              gene_type=sample(c('core', 'non-core'), 50, replace = TRUE),
              ymin=-4,
              ymax=-3,
              genome='genome3')
# tst

filt <- tst %>% filter(gene_type == 'non-core')

step4 <- split(filt$locus_tag_num, cumsum(c(1, diff(filt$locus_tag_num) != 1)))    # identifies streches of sequential loc.tags
#browser()
step5 <- step4[lapply(step4, length) >= 3]
islands <- unlist(step5)


tst$island <- ifelse(tst$locus_tag_num %in% islands, 'island', 'non-island')

genome3 <- tst



# end genomes

tst <- bind_rows(genome1, genome2, genome3)

# assign island clusters



# tst

clust1 <- c('genome1_16','genome1_17','genome1_18','genome1_19','genome2_19',
            'genome2_20','genome2_21','genome2_22','genome3_19','genome3_20',
            'genome3_21','genome3_22')

clust2 <- c('genome3_1', 'genome3_2', 'genome3_3',
            'genome3_37', 'genome3_38', 'genome3_39',
            'genome1_30', 'genome1_31', 'genome1_32',
            'genome2_30', 'genome2_31', 'genome2_32')

clust3 <- c('genome1_39', 'genome1_40', 'genome1_41', 'genome1_42',
            'genome2_12', 'genome2_13', 'genome2_14', 'genome2_15', 'genome2_16', 'genome2_17',
            'genome3_46', 'genome3_47', 'genome3_48')

clust4 <- c('genome1_44', 'genome1_45', 'genome1_46', 'genome1_47', 'genome1_48',
            'genome2_7', 'genome2_8', 'genome2_9',
            'genome3_14', 'genome3_15', 'genome3_16', 'genome3_17')


tet_resist <- c('genome1_39',
                'genome2_12',
                'genome3_46')

phage <- c('genome1_16','genome1_17','genome1_18','genome1_19',
           'genome2_19','genome2_20','genome2_21','genome2_22',
           'genome3_19','genome3_20','genome3_21','genome3_22')


tst <- tst %>% mutate(
  clust=case_when(
    locus_tag %in% clust1  ~ 'cluster1',
    locus_tag %in% clust2  ~ 'cluster2',
    locus_tag %in% clust3  ~ 'cluster3',
    locus_tag %in% clust4  ~ 'cluster4'),
  Gene_Type=case_when(
    locus_tag %in% tet_resist ~ 'tet_resistance',
    locus_tag %in% phage      ~ 'phage',
    FALSE                     ~ 'non-core')
  )

tst$gene_Type <- 'gene'

tst %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=gene_Type), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=2.5, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank()) +
  ggtitle('Step1: annotate with prokka') + 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')




tst %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=gene_type), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=2.5, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank()) +
  scale_fill_manual(values = colours) +
  ggtitle('Step2: calculate pangenome with roary')+ 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')

tst %>% filter(gene_type == 'non-core') %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=gene_type), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=3, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())+
  scale_fill_manual(values = colours) +
  ggtitle('Step3: remove core genome')+ 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')


tst %>% filter(island == 'island') %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=gene_type), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=3, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())+
  scale_fill_manual(values = colours) +
  ggtitle('Step4: remove small runs of non-core genes')+ 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')

tst %>% filter((locus_tag %in% c(clust1, clust2, clust3, clust4)) & island == 'island') %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=Gene_Type), color='black') +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=Gene_Type), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=3, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank()) +
  ggtitle('Step5: annotate islands with abricate')+ 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')



tst %>% filter((locus_tag %in% c(clust1, clust2, clust3, clust4)) & island == 'island') %>%
  ggplot(aes(x=strt, y=ymax)) +
  geom_rect(aes(xmin=strt, xmax=stp, ymin=ymin, ymax=ymax, fill=clust), color='black') +
  geom_text(aes(label=locus_tag), angle=90, size=3, hjust=1.1, vjust=1) +
  ylim(-4.5,1.5) + xlim(0,51) + theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank()) +
  ggtitle('Step6: graph based clustering')+ 
  annotate(geom='text', x=25, y=1.1, label='Genome1')+
  annotate(geom='text', x=25, y=-.9, label='Genome2')+
  annotate(geom='text', x=25, y=-2.9, label='Genome3')

```

## Known issues  

**Draft assemblies**  
gifrop can handle incomplete, draft assemblies, but the genomic islands it identifies **will never span contig breaks**.  That is, if there is a genomic island in your assembly that is divided across two contigs, it will be split into 2 different genomic islands.    
  
**Complete genomes**  
If you input complete genomes, (with closed circular chromosomes) and a genomic island spans the arbitrarily chosen start/end coordinates, it will be split into two genomic islands.  

**Combining of adjacent/nested genomic islands**  
If there is a situation where two distinct genomic islands are inserted adjacent to one another or one nested within another etc., they will be treated as one genomic island.  

**Error when identical contig names in different assemblies**  
I think running prokka with `--compliant` will avoid this.  



 
## TODO  

1) link abricate results to pangenome, use coords of hits with gffs to get loc_tags?  
2) Fix issue with viroseqs abricate db, puts 'phage' into resistance column.
3) Include example data?  


Suggestions welcome.