---
title: "PFAM domain plots"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/pfam_plots/',
  fig.width = 3.2, fig.height = 6.8, dpi = 350,
  fig.path='figure/4c-pfam_plots/',
  dev = c('png', 'pdf')
)
```

## Setup
```{r, cache=FALSE}
library(stringr)
library(tidyr)
library(dplyr)
library(glue)
library(ggplot2)
library(RColorBrewer)

theme_set(theme_classic())
theme_update(
  panel.background = element_rect(fill = "transparent", colour = NA),
  plot.background = element_rect(fill = "transparent", colour = NA),
  legend.background = element_rect(fill = "transparent", colour = NA)
)

setwd(paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/R_SCRIPT/"))
```

I loaded a tweaked version of the .txt file manually changed the header, including the hashtag column as it was freaking `read.table()` out.

check what all columns mean here:
http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf (pg 48 -- domain hits table)

relevant columns:

* 1: target name (gene/transcript), 3: target length (how much on the gene is matched)
* 3: target length
* 4, 5, 6: query/domain name, accession, qlen
* 7, 8, 9: e-value, score, bias
* 10, 11: hash, number of domains reported in the sequence (e.g. 1st out of 3 domains detected)
* 12, 13, 14, 15: c-e-value, i-e-value, score, bias
* 20, 21: to, from (alignment coordiante)
* 22: probability


The problem is that some genes are counted multiple times, because hmmsearch works at the transcript level. My rationale: Group by the gene ID (1), domain accession (5), find maximum: number of domains reported (11). This should remove duplicated hits due to a gene having multiple transcripts. Overlapping domains are acceptable.


## Load data
```{r}
astCal_raw <- read.table(
  file = paste0(Sys.getenv("HOME"), "/code/malawi_transposon/storage/cloud/other/pfam/AstCal_results_changedHeader.txt"),
  sep = "", header = TRUE, fill = TRUE, quote = NULL
)

mayZeb_raw <- read.table(
  file = paste0(Sys.getenv("HOME"), "/code/malawi_transposon/storage/cloud/other/pfam/MayZeb_results_changedHeader.txt"),
  sep = "", header = TRUE, fill = TRUE, quote = NULL
)
```

Private genes that are confirmed to have no ortholog in the other species through BioMart.
```{r}
df_compared_astCal <- readr::read_csv(
  paste0(Sys.getenv("HOME"), "/code/malawi_transposon/storage/local/script/gene_ontology/geneListFilter_markedlyAbsent-astCal.csv")
)
astCal_realPrivate_genes <- filter(df_compared_astCal, is.na(ortholog)) %>% .$gene %>% unique() 
astCal_Artifactual_genes <- filter(df_compared_astCal, !is.na(ortholog)) %>% .$gene %>% unique()

df_compared_mayZeb <- readr::read_csv(
  paste0(Sys.getenv("HOME"), "/code/malawi_transposon/storage/local/script/gene_ontology/geneListFilter_markedlyAbsent-mayZeb.csv")
)
mayZeb_realPrivate_genes <- filter(df_compared_mayZeb, is.na(ortholog)) %>% .$gene %>% unique() 
mayZeb_Artifactual_genes <- filter(df_compared_mayZeb, !is.na(ortholog)) %>% .$gene %>% unique()
```

## astCal
```{r}
pfam_astCal <- astCal_raw[,c(1,3:13,20:22)]

# de-duplication
pfam_astCal <- pfam_astCal %>% 
  mutate(gene = str_extract(target_name, "^[^|]+") ) %>% 
  rowwise() %>% 
  mutate(id = glue(target_name, '-', accession.1)) %>% 
  select(id, gene, target_name, hash, of, accession.1, query_name, qlen)
non_overlapping_ids_AstCal <- pfam_astCal %>% group_by(gene, accession.1) %>% mutate(id_max = id[which.max(of)])
pfam_astCal2 <- pfam_astCal %>% filter(id %in% unique(non_overlapping_ids_AstCal$id_max) )
```

Bar chart of top domains in truly private genes.
```{r pfam-astCal}
# reorder domains by count
freq_proteinDomAstCal <- count(filter(pfam_astCal2, gene %in% astCal_realPrivate_genes), query_name)
freq_proteinDomAstCal <- freq_proteinDomAstCal[order(-freq_proteinDomAstCal$n), ]
freq_proteinDomAstCal$query_name <- factor(
  freq_proteinDomAstCal$query_name, levels = freq_proteinDomAstCal$query_name[order(-freq_proteinDomAstCal$n)]
)

freq_proteinDomAstCal %>% 
  filter(n > 10) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#D65F5F', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```

Bar chart of top domains in artefacts.
```{r pfam-astCal-artifact}
# reorder domains by count
freq_proteinDomAstCal <- count(filter(pfam_astCal2, gene %in% astCal_Artifactual_genes), query_name)
freq_proteinDomAstCal <- freq_proteinDomAstCal[order(-freq_proteinDomAstCal$n), ]
freq_proteinDomAstCal$query_name <- factor(
  freq_proteinDomAstCal$query_name, levels = freq_proteinDomAstCal$query_name[order(-freq_proteinDomAstCal$n)]
)

freq_proteinDomAstCal %>% 
  filter(n > 5) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#D65F5F', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```

Plot everything together
```{r pfam-astCal-everything}
# reorder domains by count
freq_proteinDomAstCal <- count(pfam_astCal2, query_name)
freq_proteinDomAstCal <- freq_proteinDomAstCal[order(-freq_proteinDomAstCal$n), ]
freq_proteinDomAstCal$query_name <- factor(
  freq_proteinDomAstCal$query_name, levels = freq_proteinDomAstCal$query_name[order(-freq_proteinDomAstCal$n)]
)

freq_proteinDomAstCal %>% 
  filter(n > 11) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#D65F5F', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```

## mayZeb
```{r}
pfam_mayZeb <- mayZeb_raw[,c(1,3:13,20:22)]

# de-duplication
pfam_mayZeb <- pfam_mayZeb %>% 
  mutate(gene = str_extract(target_name, "^[^|]+") ) %>% 
  rowwise() %>% 
  mutate(id = glue(target_name, '-', accession.1)) %>% 
  select(id, gene, target_name, hash, of, accession.1, query_name, qlen)
non_overlapping_ids_AstCal <- pfam_mayZeb %>% group_by(gene, accession.1) %>% mutate(id_max = id[which.max(of)])
pfam_mayZeb2 <- pfam_mayZeb %>% filter(id %in% unique(non_overlapping_ids_AstCal$id_max) )
```

Bar chart of top domains in truly private genes.
```{r pfam-mayZeb}
# reorder domains by count
freq_proteinDomMayZeb <- count(filter(pfam_mayZeb2, gene %in% mayZeb_realPrivate_genes), query_name)
freq_proteinDomMayZeb <- freq_proteinDomMayZeb[order(-freq_proteinDomMayZeb$n), ]
freq_proteinDomMayZeb$query_name <- factor(
  freq_proteinDomMayZeb$query_name, levels = freq_proteinDomMayZeb$query_name[order(-freq_proteinDomMayZeb$n)]
)

freq_proteinDomMayZeb %>% 
  filter(n > 1) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#82C6E2', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```

Bar chart of top domains in artefacts.
```{r pfam-mayZeb-artifact}
# reorder domains by count
freq_proteinDomMayZeb <- count(filter(pfam_mayZeb2, gene %in% mayZeb_Artifactual_genes), query_name)
freq_proteinDomMayZeb <- freq_proteinDomMayZeb[order(-freq_proteinDomMayZeb$n), ]
freq_proteinDomMayZeb$query_name <- factor(
  freq_proteinDomMayZeb$query_name, levels = freq_proteinDomMayZeb$query_name[order(-freq_proteinDomMayZeb$n)]
)

freq_proteinDomMayZeb %>% 
  filter(n >= 5) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#82C6E2', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```

Plot everything together
```{r pfam-mayZeb-everything}
# reorder domains by count
freq_proteinDomMayZeb <- count(pfam_mayZeb2, query_name)
freq_proteinDomMayZeb <- freq_proteinDomMayZeb[order(-freq_proteinDomMayZeb$n), ]
freq_proteinDomMayZeb$query_name <- factor(
  freq_proteinDomMayZeb$query_name, levels = freq_proteinDomMayZeb$query_name[order(-freq_proteinDomMayZeb$n)]
)

freq_proteinDomMayZeb %>% 
  filter(n >= 5) %>% 
  ggplot(aes(x = query_name, y = n, width=0.70)) +
  geom_bar(stat = "identity", fill = '#82C6E2', color = 'black') + 
  labs(y = "Count", x = NULL) + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(expand = expansion(mult = c(.05, .25)) ) +
  geom_text(stat="identity", aes(label=n), hjust=-0.5, vjust=0.5) +
  coord_flip() +
  theme_bw()
```
