---
title: "RefPlantNLR appendix S1"
author: "J Kourelis"
output: 
  html_document:
    toc: TRUE
---

```{r libraries, echo=FALSE, warning=FALSE, message=FALSE}
library(tidyverse)
library(readxl)
library(knitr)
library(splitstackshape)
library(RColorBrewer)
```

# Description

R script used to generate the gff3 output for the RefPlantNLRs as well as the figures.

# The script to generate the NLR annotations

## InterProScan

* Run InterProScan (used version v5.42-78.0) with the following code:

```text
interproscan.sh -i RefNLRs_AA.fasta -appl Pfam,Gene3D,SUPERFAMILY,PRINTS,Coils,SMART,CDD,ProSiteProfiles -f gff3 -cpu 10 -dp
```

* Optional: Run the following R script to compress the size of the GFF files
* Remove second overlapping SSF52540 signature from DM2h^Ler going from position 268-649

```{r, echo=FALSE, warning=FALSE, message=FALSE}

sample <- read.delim("Compress/RefNLRs_AA.fasta.gff3", header=FALSE, sep="\t", stringsAsFactors = FALSE)  %>%
  filter(!grepl("#", V1),
         V3 != "polypeptide",
         !is.na(V4))

write.table(sample, "RefNLRs/InterProScan.gff3", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = FALSE)

rm(sample)
```

## Read in the proteome annotation
  
* Read in the proteome annotations (stored in the folder gff/RefNLRs)

```{r, echo=FALSE, warning=FALSE, message=FALSE}

RefNLRs <- read.delim("RefNLRs/InterProScan.gff3", header=FALSE, sep="\t", stringsAsFactors = FALSE) %>%
  rename(seqname = V1,  # Name of the chromosome or scaffold
         source = V2,  # Name of the program that generated this feature, or the data source (database or project name)
         feature = V3,  # Feature type name, e.g. Gene, Variation, Similarity
         start = V4,  # Start position of the feature, with sequence numbering starting at 1
         end = V5,  #  End position of the feature, with sequence numbering starting at 1
         score = V6,  # A floating point value
         strand = V7,  # Defined as + (forward) or - (reverse)
         frame = V8,  # One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on
         attribute = V9) %>%  # A semicolon-separated list of tag-value pairs, providing additional information about each feature
  filter(!grepl("#", seqname),
         feature != "polypeptide",
         !is.na(start)) %>%
  mutate(seqname = str_replace(seqname, "Bph19", "Bph1/9"))

```

## Make pretty GFF file with annotation

* Requirements (all stored in the folder Description as .tsv files):
    + [CATH-Gene3D](https://www.ebi.ac.uk/interpro/beta/entry/unintegrated/#table)
    + [Pfam](https://www.ebi.ac.uk/interpro/beta/entry/unintegrated/Pfam/#table)    
    + [SUPERFAMILY](https://www.ebi.ac.uk/interpro/beta/entry/unintegrated/ssf/#table)
    + [InterPro](ftp://ftp.ebi.ac.uk/pub/databases/interpro/entry.list)

```{r, echo=FALSE, warning=FALSE, message=FALSE}

# GFF annotation of NLRs

## InterProScan description

InterProScan <- read.delim("Description/20200423_InterPro_79.0_InterPro.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) 

## SMART annotation

### SMART gff

NLR_SMART <- RefNLRs %>%
  filter(source == "SMART")

### SMART description

InterPro_SMART <- read.delim("Description/20200423_InterPro_79.0_SMART.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,")) %>%
  mutate_all(na_if,"")

### SMART gff annotated

NLR_SMART <- NLR_SMART %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_4 = str_replace(attribute_4, "Name=", ""),
         attribute_6 = str_replace(attribute_6, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_SMART, by = c("attribute_4" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_4)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";SMART=", attribute_4, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## CDD annotation

### CDD gff

NLR_CDD <- RefNLRs %>%
  filter(source == "CDD")

### CDD description

InterPro_CDD <- read.delim("Description/20200423_InterPro_79.0_CDD.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,")) %>%
  mutate_all(na_if,"")

### CDD gff annotated

NLR_CDD <- NLR_CDD %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_5 = str_replace(attribute_5, "Name=", ""),
         attribute_7 = str_replace(attribute_7, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_CDD, by = c("attribute_5" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_5)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";CDD=", attribute_5, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## ProSiteProfiles annotation

### ProSiteProfiles gff

NLR_ProSiteProfiles <- RefNLRs %>%
  filter(source == "ProSiteProfiles")

### ProSiteProfiles description

InterPro_ProSiteProfiles <- read.delim("Description/20200423_InterPro_79.0_PrositeProfiles.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,")) %>%
  mutate_all(na_if,"")

### ProSiteProfiles gff annotated

NLR_ProSiteProfiles <- NLR_ProSiteProfiles %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_5 = str_replace(attribute_5, "Name=", ""),
         attribute_7 = str_replace(attribute_7, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_ProSiteProfiles, by = c("attribute_5" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_5)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";ProSiteProfiles=", attribute_5, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## Coils annotation

### Coils gff

NLR_Coils <- RefNLRs %>%
  filter(source == "Coils") %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_4 = str_replace(attribute_4, "Name=", ""),
         attribute = paste0("Name=", attribute_4, ";Coils=", attribute_4, ";note=", attribute_4, ";InterProScan=None")) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## Gene3D annotation

### Gene3D gff

NLR_Gene3D <- RefNLRs %>%
  filter(source == "Gene3D")

### Gene3D description

InterPro_Gene3D <- read.delim("Description/20200423_InterPro_79.0_Gene3D.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,"),
         ENTRY_NAME = if_else(Accession == "G3DSA:3.30.70.100", "HMA", ENTRY_NAME)) %>%
  mutate_all(na_if,"")

### Gene3D gff annotated

NLR_Gene3D <- NLR_Gene3D %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_4 = str_replace(attribute_4, "Name=", ""),
         attribute_6 = str_replace(attribute_6, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_Gene3D, by = c("attribute_4" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_4)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";Gene3D=", attribute_4, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## Pfam annotation

### Pfam gff

NLR_Pfam <- RefNLRs %>%
  filter(source == "Pfam") 

### Pfam description

InterPro_Pfam <- read.delim("Description/20200423_InterPro_79.0_Pfam.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,")) %>%
  mutate_all(na_if,"")

### Pfam gff annotated

NLR_Pfam <- NLR_Pfam %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_5 = str_replace(attribute_5, "Name=", ""),
         attribute_7 = str_replace(attribute_7, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_Pfam, by = c("attribute_5" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_5)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";Pfam=", attribute_5, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## SUPERFAMILY annotation

### SUPERFAMILY gff

NLR_SUPERFAMILY <- RefNLRs %>%
  filter(source == "SUPERFAMILY") 

### SUPERFAMILY annotations

InterPro_SUPERFAMILY <- read.delim("Description/20200423_InterPro_79.0_SUPERFAMILY.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,")) %>%
  mutate_all(na_if,"")

### SUPERFAMILY gff annotated

NLR_SUPERFAMILY <- NLR_SUPERFAMILY %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_4 = str_replace(attribute_4, "Name=", ""),
         attribute_6 = str_replace(attribute_6, "Dbxref=InterPro:", "")) %>%
  left_join(InterPro_SUPERFAMILY, by = c("attribute_4" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_4)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";SUPERFAMILY=", attribute_4, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

## PRINTS annotation

### PRINTS gff

NLR_PRINTS <- RefNLRs %>%
  filter(source == "PRINTS") 

### PRINTS description

InterPro_PRINTS <- read.delim("Description/20200423_InterPro_79.0_PRINTS.tsv", header=TRUE, sep="\t", stringsAsFactors = FALSE) %>%
  left_join(InterProScan, by = c("Integrated.Into" = "ENTRY_AC")) %>%
  select(Accession, Name, Integrated.Into, ENTRY_NAME) %>%
  mutate(Name = str_replace(Name, "\\;", "\\,"),
         ENTRY_NAME = if_else(Accession == "PR00364", "DISEASERSIST", ENTRY_NAME)) %>%
  mutate_all(na_if,"")

### PRINTS gff annotated

NLR_PRINTS <- NLR_PRINTS %>%
  cSplit("attribute", sep = ";", direction = "wide") %>%
  mutate(feature = source) %>%
  mutate(attribute_5 = str_replace(attribute_5, "Name=", "")) %>%
  left_join(InterPro_PRINTS, by = c("attribute_5" = "Accession")) %>%
  mutate(Name = if_else(!is.na(ENTRY_NAME), ENTRY_NAME,
                        if_else(!is.na(Name), Name, attribute_5)),
         InterProScan = if_else(!is.na(Integrated.Into), Integrated.Into, "None"),
         attribute = paste0("Name=", Name, ";PRINTS=", attribute_5, ";note=", Name, ";InterProScan=", InterProScan)) %>%
  select(seqname, source, feature, start, end, score, strand, frame, attribute)

##  all together

NLRs_combined <- rbind(NLR_Gene3D, NLR_Pfam, NLR_SUPERFAMILY, NLR_PRINTS, NLR_SMART, NLR_CDD, NLR_Coils, NLR_ProSiteProfiles)

NLRs_combined %>%
  write.table("RefNLRs/20200528_RefNLRs_AA.gff3", sep="\t", row.names = FALSE, col.names = FALSE, na = "", quote=FALSE)

NLRs_combined %>%
  mutate(start = (start*3)-2,
         end = end*3) %>%
  filter(!grepl("AhRRS5|Le4|Pid3-W3|Pid3-W4|Pid3-I1|Pi64|R9a|Rpi-mcq1.1|Rpi-sto2|NbNRC4a|NbNRC4b|Ag-0|RPS4-Ws|RRS1-Ws|Pik-e", seqname)) %>% # no nucleotide sequences
  write.table("RefNLRs/20200528_RefNLRs_CDS.gff3", sep="\t", row.names = FALSE, col.names = FALSE, na = "", quote=FALSE)

rm(InterPro_Gene3D, InterPro_SUPERFAMILY, InterPro_Pfam, InterPro_PRINTS, InterPro_SMART, InterPro_CDD, InterPro_ProSiteProfiles, InterProScan, NLR_Gene3D, NLR_Pfam, NLR_SUPERFAMILY, NLR_PRINTS, NLR_SMART, NLR_CDD, NLR_Coils, NLR_ProSiteProfiles, NLRs, NLRs_combined, NACHT)
```

# Generate the NB-ARC phylogeny

* extract the SSF52540 using Geneious Prime, export to RefNLRs/SSF52540.csv
* Deduplicate identical SSF52540 domains:

```{r, echo=FALSE, warning=FALSE, message=FALSE}

SSF52540 <- read.csv("RefNLRs/SSF52540.csv") %>%
  mutate(Name = str_replace(Name, " - .*", "")) %>%
  group_by(Name) %>%
  mutate(ID = row_number()) %>% # Identify NLRs with multiple NB-ARC domains
  ungroup() %>%
  mutate(Name = if_else(ID != "1", paste0(Name, "_", ID), as.character(Name))) %>% # Give unique identifier to NB-ARC from NLRs with multiple NB-ARC domains
  select(-ID) %>%
  group_by(Sequence) %>%
  arrange(Name) %>%
  summarise(Name = paste(Name, collapse = "//")) %>% # Collapse identical NB-ARC domains with a // separator
  ungroup() %>%
  select(Name, Sequence)

write.csv(SSF52540, "RefNLRs/RefNLRs_SSF52540_DEDUP.csv", row.names = FALSE, quote = FALSE)

```

* Align using ClustalOmega
* Use QKphylogeny to remove data missing in >= 95% of positions

```text
python QKphylogeny_alignment_analysis.py -a RefNLRs_SSF52540_DEDUP_ClustalOmega.phy -d 0.05 -t protein -o RefNLRs_SSF52540_DEDUP_ClustalOmega.phy_missing5.phy
```

* Use RAxML to generate the NB-ARC phylogeny

```text
raxmlHPC-PTHREADS-AVX -T 6 -s RefNLRs_SSF52540_DEDUP_ClustalOmega_missing5.phy -n 20200528_RefNLR -m PROTGAMMAAUTO -f a -# 100 -x 8153044963028367 -p 644124967711489 
```

# Generate the reduced redundancy set

* Use Geneious Prime to extract the NB-ARC domain (SSF52540)
* Use CD-HIT to generate reduced redundancy set based on the extracted NB-ARC domains

```text
cdhit -i RefNLRs_SSF52540.fasta -o RefNLRs_SSF52540_90 -c 0.90 -n 5 -M 16000 -T 10 -d 0
```

* Use CD-HIT output to assign representative sequences per genus

```{r, echo=FALSE, warning=FALSE, message=FALSE}

# Get all the clusters from the CD-HIT output
Cluster <- read.delim("RefNLRs_SSF52540_90.clstr", header=FALSE, sep="\t", stringsAsFactors = FALSE) %>%
  cSplit("V2", sep = ", >", direction = "wide") %>%
  cSplit("V2_2", sep = "... ", direction = "wide") %>%
  rename(Cluster = V1,  # query (e.g., gene) sequence id
         Length_AA = V2_1,  # subject (e.g., reference genome) sequence id
         seqname = V2_2_1,  # percentage of identical matches
         pident = V2_2_2) %>%  # alignment length
  mutate(Length_AA = str_replace(Length_AA, "aa", ""),
         seqname = str_replace(seqname, ">", ""),
         seqname = str_replace(seqname, "\\.\\.", ""),
         pident = str_replace(pident, "at ", ""),
         pident = str_replace(pident, "%", ""),
         Representative = case_when(pident == "*" ~ seqname, # Name the representative of each cluster
                                    TRUE ~ NA_character_),
         Cluster = case_when(str_detect(Cluster, ">") ~ Cluster,
                                    TRUE ~ NA_character_),
         Cluster = as.numeric(str_replace(Cluster, ">Cluster ", "")),
         pident = str_replace(pident, "\\*", "100.00")) %>% # Set 100% identity for the cluster representative
  fill(Cluster, .direction = "down") %>%
  group_by(Cluster) %>%
  arrange(Representative) %>%
  fill(Representative, .direction = "down") %>%
  ungroup() %>%
  filter(!is.na(pident)) %>%
  mutate(Cluster = paste("Cluster", formatC(format(Cluster, scientific = FALSE), width = 5, flag = "0"), sep = "_"),
         Cluster = str_replace_all(Cluster, " ", "0")) %>% # Consistently name the clusters
  group_by(seqname) %>% # If NLR has multiple NB-ARC select largest and cluster accordingly
  arrange(desc(Length_AA)) %>%
  slice(1) 

# Get out the representative
Representative_Part1 <- read_xlsx("RefNLRs/20200528_Table_KnownNLRs.xlsx", na = "NA") %>%
  filter(ID != "MLA38-1") %>%
  mutate(Genus = str_remove(Host, "\\s.*")) %>%
  select(seqname = ID, Genus) %>%
  left_join(Cluster, by = "seqname") %>%
  filter(seqname == Representative) # Select the representative for each cluster

# Retain sequences as representative if multuple genera were clustered together in a single cluster
Representative_Part2 <- read_xlsx("RefNLRs/20200528_Table_KnownNLRs.xlsx", na = "NA") %>%
  filter(ID != "MLA38-1") %>%
  mutate(Genus = str_remove(Host, "\\s.*")) %>%
  select(seqname = ID, Genus) %>%
  left_join(Cluster, by = "seqname") %>%
  filter(!seqname %in% Representative_Part1$seqname) %>% # Remove all the representative sequences
  group_by(Cluster, Genus) %>% # select the largest sequence per Cluster/Genus
  arrange(desc(Length_AA)) %>%
  slice(1) %>%
  anti_join(Representative_Part1, by = c("Cluster", "Genus")) %>% # Remove the sequences in each cluster belonging to the same Genus as the representative sequence
  mutate(Representative = seqname)

# Final set
## Select all the representative sequences, label them
Representative <- bind_rows(Representative_Part1, Representative_Part2) %>%
  select(ID = Representative, Cluster) %>%
  mutate(Representative = "YES")

## Add the non-representative sequences, label them
Representative <- Cluster %>%
  select(ID = seqname, Cluster) %>%
  filter(!ID %in% Representative$ID) %>%
  mutate(Representative = "NO") %>%
  bind_rows(Representative)

rm(Cluster, Representative_Part1, Representative_Part2)

```

# Generate figures

```{r, echo=FALSE, warning=FALSE, message=FALSE}

RefNLRs_AA <- read.csv("RefNLRs/RefPlantNLR_v20200528_416_AA.csv")

RefNLRs_SSF52540 <- read.csv("RefNLRs/RefPlantNLR_SSF52540_v20200528_416_AA.csv")

RefNLRs_INFO <- read_xlsx("RefNLRs/20200528_Table_KnownNLRs.xlsx", na = "NA") %>%
  mutate(Genus = str_remove(Host, "\\s.*"),
         Class = case_when(grepl("CNL", Subclass) ~ "CC-NLR",
                           grepl("GT", Subclass) ~ "TIR-NLR",
                           grepl("GR", Subclass) ~ "CCR-NLR",
                           grepl("G10", Subclass) ~ "CCG10-NLR",
                           TRUE ~ "Other"),
         Structure = case_when(Class == "CCG10-NLR" ~  paste0(Structure, " (CCG10-NLR)"),
                               Class == "CCR-NLR" & Structure == "NL" ~ "NL (CCR-NLR)",
                               Class == "CC-NLR" & Structure == "NL" ~ "NL (CC-NLR)",
                               Class == "CC-NLR" & Structure == "NLO" ~ "NLO (CC-NLR)",
                               Class == "CC-NLR" & Structure == "ONL" ~ "ONL (CC-NLR)",
                               TRUE ~ Structure)) %>%
  left_join(Representative, by = "ID")

RefNLRs_INFO %>% 
  left_join(RefNLRs_AA, by = c("ID" = "Name")) %>%
  distinct(Sequence, .keep_all = TRUE) %>%
  group_by(Class) %>%
  summarise(n = n(),
            median = median(Length_AA),
            mean = mean(Length_AA),
            min = min(Length_AA),
            max = max(Length_AA),
            stdev = sd(Length_AA)) %>%
  mutate(SE = stdev / sqrt(n),
         lower.ci = mean - qt(1 - (0.05 / 2), n - 1) * SE,
         upper.ci = mean + qt(1 - (0.05 / 2), n - 1) * SE)

RefNLRs_INFO %>% 
  left_join(RefNLRs_SSF52540, by = c("ID" = "Name")) %>%
  distinct(Sequence, .keep_all = TRUE) %>%
  mutate(Length_SSF52540 = str_length(Sequence)) %>%
  summarise(n = n(),
            median = median(Length_SSF52540),
            mean = mean(Length_SSF52540),
            min = min(Length_SSF52540),
            max = max(Length_SSF52540),
            stdev = sd(Length_SSF52540)) %>%
  mutate(SE = stdev / sqrt(n),
         lower.ci = mean - qt(1 - (0.05 / 2), n - 1) * SE,
         upper.ci = mean + qt(1 - (0.05 / 2), n - 1) * SE)

## Figure 1A
RefNLRs_stack <- RefNLRs_INFO %>%
  group_by(Genus, Class) %>%
  summarise(n = n())
 
ggplot(data = RefNLRs_stack, aes(x = Genus, y = n, fill = Class, label = n)) + 
  geom_bar(stat = "identity") +
  geom_text(size = 3, position = position_stack(vjust = 0.5)) +
  scale_fill_brewer(palette="Set3") +
  guides(fill = guide_legend(title = NULL)) + 
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_discrete(limits = c("Hordeum", "Secale", "Triticum", "Aegilops", "Dasypyrum", "Oryza", "Sorghum", "Zea", "Musa", "Beta", "Spinacia", "Helianthus", "Lactuca", "Nicotiana", "Capsicum", "Solanum", "Vitis", "Linum", "Cucumis", "Prunus", "Malus", "Rosa", "Arachis", "Medicago", "Cajanus", "Glycine", "Phaseolus", "Gossypium", "Brassica", "Arabidopsis")) + 
  coord_flip() +
  ylab("# Functionally validated NLRs") + 
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.y = element_line(colour = NA),
        axis.line.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.y = element_text(angle = 0, hjust = 1, size = 8),
        axis.title.y=element_blank())

## Figure 1B
RefNLRs_stack_dedup <- RefNLRs_INFO %>%
  filter(Representative == "YES") %>%
  group_by(Genus, Class) %>%
  summarise(n = n())
 
ggplot(data = RefNLRs_stack_dedup, aes(x = Genus, y = n, fill = Class, label = n)) + 
  geom_bar(stat = "identity") +
  geom_text(size = 3, position = position_stack(vjust = 0.5)) +
  scale_fill_brewer(palette="Set3") +
  guides(fill = guide_legend(title = NULL)) + 
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_discrete(limits = c("Hordeum", "Secale", "Triticum", "Aegilops", "Dasypyrum", "Oryza", "Sorghum", "Zea", "Musa", "Beta", "Spinacia", "Helianthus", "Lactuca", "Nicotiana", "Capsicum", "Solanum", "Vitis", "Linum", "Cucumis", "Prunus", "Malus", "Rosa", "Arachis", "Medicago", "Cajanus", "Glycine", "Phaseolus", "Gossypium", "Brassica", "Arabidopsis")) + 
  coord_flip() +
  ylab("# Functionally validated NLRs") + 
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.y = element_line(colour = NA),
        axis.line.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.y = element_text(angle = 0, hjust = 1, size = 8),
        axis.title.y=element_blank())

## Figure 2A and 3A
RefNLRs_Unique <- RefNLRs_INFO %>%
  left_join(RefNLRs_AA, by = c("ID" = "Name")) %>%
  distinct(Sequence, .keep_all = TRUE)

RefNLRs_Unique %>%
  group_by(Class) %>%
  summarise(n = n())

ggplot(data = RefNLRs_Unique, aes(x = Length_AA, fill = Class)) +
  geom_histogram(binwidth = 50) +
  ggtitle("Unique NLR length") + 
  scale_fill_brewer(palette="Set3") +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_continuous(breaks = c(0, 500, 1000, 1500, 2000, 2500)) +
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12),
        axis.title.x=element_blank())

ggplot(data = RefNLRs_Unique, aes(x = Structure)) + 
  geom_bar(stat = "count") +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -1) + 
  ggtitle("Unique NLRs / domain architecture") + 
  scale_x_discrete(limits = c("CNL", "NL (CC-NLR)", "ONL (CC-NLR)", "NLO (CC-NLR)", "CNLO", "OCNLO", "CONL", "BCCNL", "BCNL", "CNNL", "NL (CCG10-NLR)", "ONL (CCG10-NLR)", "RNL", "NL (CCR-NLR)", "RNOO", "TNL", "TNLN", "TN", "TNOOO", "TNLO", "TNLOO", "OOOTNLO", "TNOTNL")) + 
  scale_fill_brewer(palette="Set3") +
  guides(fill = guide_legend(title = NULL)) + 
  ylab("count") + 
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 60, hjust = 1, size = 12),
        axis.title.x=element_blank())

## Figure 2C and 3B
RefNLRs_Dedup <- RefNLRs_INFO %>%
  filter(Representative == "YES")

ggplot(data = RefNLRs_Dedup, aes(x = Length_AA, fill = Class)) +
  geom_histogram(binwidth = 50) +
  ggtitle("Unique NLR length") + 
  scale_fill_brewer(palette="Set3") +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_continuous(breaks = c(0, 500, 1000, 1500, 2000, 2500)) +
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12),
        axis.title.x=element_blank())

ggplot(data = RefNLRs_Dedup, aes(x = Structure)) + 
  geom_bar(stat = "count") +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -1) + 
  ggtitle("Unique NLRs / domain architecture") + 
  scale_x_discrete(limits = c("CNL", "NL (CC-NLR)", "ONL (CC-NLR)", "NLO (CC-NLR)", "CNLO", "OCNLO", "CONL", "BCCNL", "BCNL", "CNNL", "NL (CCG10-NLR)", "ONL (CCG10-NLR)", "RNL", "NL (CCR-NLR)", "RNOO", "TNL", "TNLN", "TN", "TNOOO", "TNLO", "TNLOO", "OOOTNLO", "TNOTNL")) + 
  scale_fill_brewer(palette="Set3") +
  guides(fill = guide_legend(title = NULL)) + 
  ylab("count") + 
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 60, hjust = 1, size = 12),
        axis.title.x=element_blank())

## Figure 2B
RefNLRs_SSF52540_Unique <- RefNLRs_INFO %>%
  left_join(RefNLRs_SSF52540, by = c("ID" = "Name")) %>%
  mutate(Length_SSF52540 = str_length(Sequence)) %>%
  distinct(Sequence, .keep_all = TRUE)

RefNLRs_SSF52540_Unique %>%
  group_by(Class) %>%
  summarise(n = n())

ggplot(data = RefNLRs_SSF52540_Unique, aes(x = Length_SSF52540, fill = Class)) +
  geom_histogram(binwidth = 5) +
  ggtitle("Unique NB-ARC length") + 
  scale_fill_brewer(palette="Set3") +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_continuous(breaks = c(150, 200, 250, 300, 350)) +
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12),
        axis.title.x=element_blank())

## Figure 2D
RefNLRs_SSF52540_Dedup <- RefNLRs_INFO %>%
  left_join(RefNLRs_SSF52540, by = c("ID" = "Name")) %>%
  mutate(Length_SSF52540 = str_length(Sequence)) %>%
  filter(Representative == "YES")

ggplot(data = RefNLRs_SSF52540_Dedup, aes(x = Length_SSF52540, fill = Class)) +
  geom_histogram(binwidth = 5) +
  ggtitle("Unique NB-ARC length") + 
  scale_fill_brewer(palette="Set3") +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)) +
  scale_x_continuous(breaks = c(150, 200, 250, 300, 350)) +
  theme(panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "#BFBFBF", size = 0.6),
        panel.grid.major.x = element_line(colour = NA),
        axis.line.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.y = element_line(color = "#BFBFBF", size = 0.6),
        axis.ticks.x = element_line(color = "#BFBFBF", size = 0.6),
        axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12),
        axis.title.x=element_blank())

```

# NLRs per taxa

* Download [Appendix S1](https://bsapubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2Fajb2.1019&file=ajb21019-sup-0001-AppendixS1.docx) from [Smith and Brown, 2018](https://bsapubs.onlinelibrary.wiley.com/doi/full/10.1002/ajb2.1019)
* Convert to .xlsx

```{r, echo=FALSE, warning=FALSE, message=FALSE}

# The RefPlantNLRs
RefNLRs_INFO <- read_xlsx("RefNLRs/20200528_Table_KnownNLRs.xlsx", na = "NA") %>%
  mutate(Genus = str_remove(Host, "\\s.*"),
         Class = case_when(grepl("CNL", Subclass) ~ "CC-NLR",
                           grepl("GT", Subclass) ~ "TIR-NLR",
                           grepl("GR", Subclass) ~ "CCR-NLR",
                           grepl("G10", Subclass) ~ "CCG10-NLR",
                           TRUE ~ "Other"),
         Structure = case_when(Class == "CCG10-NLR" ~  paste0(Structure, " (CCG10-NLR)"),
                               Class == "CCR-NLR" & Structure == "NL" ~ "NL (CCR-NLR)",
                               Class == "CC-NLR" & Structure == "NL" ~ "NL (CC-NLR)",
                               Class == "CC-NLR" & Structure == "NLO" ~ "NLO (CC-NLR)",
                               Class == "CC-NLR" & Structure == "ONL" ~ "ONL (CC-NLR)",
                               TRUE ~ Structure))

# Smith and Brown major taxa
Smith_Brown <- read_xlsx("RefNLRs/Smith_Brown_2018_TableS1.xlsx") %>%
  left_join(RefNLRs_INFO, by = c("Clade" = "Order")) %>%
  select(Clade, Class, ID) %>%
  filter(!is.na(ID)) %>%
  group_by(Clade, Class) %>%
  summarise(n = n())

Smith_Brown <- read_xlsx("RefNLRs/Smith_Brown_2018_TableS1.xlsx") %>%
  select(Clade) %>%
  left_join(Smith_Brown, by = "Clade") %>%
  spread(Class, n)

# Write to file
write.csv(Smith_Brown, "RefNLRs/BioRxiv_v1/Table_S2.csv", row.names = FALSE, quote = FALSE, na = "0")

```

