LDA_flu_test.Rmd

---
title: "LDA_flu_test"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## LDA 

https://tgmstat.wordpress.com/2014/01/15/computing-and-visualizing-lda-in-r/
https://gist.github.com/thigm85/8424654 

```{r}
library(MASS)
library(ggplot2)
library(scales)
```

```{r}
cpal <- c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#e6e600',
          '#a65628','#f781bf','#4d4d4d')
```


flu virus all test set 1  
```{r}
# label and seq data
uniflna_test1 <- 
  read.csv("/home/brian/Documents/data_mining/flu_project/uniflna_test1.csv", na.strings="")

# comparisons from sourmash
uniflna_test1_cmp <- 
  read.csv("/home/brian/Documents/data_mining/flu_project/uniflna_test1_cmp.csv")

# Label the rows
rownames(uniflna_test1_cmp) <- colnames(uniflna_test1_cmp)

# add gene column
uniflna_test1_wgs <- uniflna_test1_cmp
uniflna_test1_wgs$gene <- uniflna_test1$gene
uniflna_test1_wgs$segment <- factor(uniflna_test1$segment)

# Transform for plotting
uniflna_test1_cmp_mat <- as.matrix(uniflna_test1_cmp)
```



```{r}
lda <- lda(gene ~ ., 
           uniflna_test1_wgs[,-1002])

prop.lda = lda$svd^2/sum(lda$svd^2) 

plda <- predict(object = lda,
                newdata = uniflna_test1_wgs)
```



```{r}
dataset = data.frame(gene = uniflna_test1_wgs[,1001],
                     lda = plda$x)
```

```{r}
p1 <- ggplot(dataset) + geom_point(aes(lda.LD1, lda.LD2, colour = gene), size=2, alpha=0.2) + 
  #theme_minimal() +
  labs(x = paste("LD1 (", percent(prop.lda[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda[2]), ")", sep="")) +
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))
# https://stackoverflow.com/questions/5290003/how-to-set-legend-alpha-with-ggplot2

# ggsave("gene_LDA_plot.png", plot = p1, width = 7, height = 4, dpi = 200)

p1
```




```{r}
lda2 <- lda(segment ~ ., 
           uniflna_test1_wgs[,-1001])

prop.lda2 = lda2$svd^2/sum(lda2$svd^2) 

plda2 <- predict(object = lda2,
                newdata = uniflna_test1_wgs)
```

segment labels

```{r}
dataset2 = data.frame(segment = uniflna_test1_wgs[,1002],
                     lda2 = plda2$x)
```

```{r}
p2 <- ggplot(dataset2) + geom_point(aes(lda2.LD1, lda2.LD2, colour = segment), size=2, alpha=0.2) + 
  #theme_minimal() +
  labs(x = paste("LD1 (", percent(prop.lda2[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda2[2]), ")", sep="")) + 
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))
p2
```


segment labels not as reliable.

### Tsne for fun

```{r}
library(Rtsne)
```


```{r}
tsne_model <- Rtsne(uniflna_test1_cmp_mat, check_duplicates=FALSE, pca=TRUE, perplexity=50, theta=0.25, dims=2)
d_tsne = as.data.frame(tsne_model$Y) 
d_tsne$gene <- uniflna_test1_wgs[,1001]
#plot(d_tsne$V1, d_tsne$V2)
```

```{r}
ggplot(d_tsne, aes(V1, V2, colour = gene)) + geom_point(size=2, alpha=0.2) + 
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))
```


```{r}
itercosts <- data.frame(iter = 1:20, cost = tsne_model$itercosts)

ggplot(itercosts, aes(iter,cost)) + geom_line(color='blue')
```


second sample with 2500 sequences, k=7
uniflna_test2_cmp.csv

```{r}
# label and seq data
uniflna_test2 <- 
  read.csv("/home/brian/Documents/data_mining/flu_project/uniflna_test2.csv", na.strings="")

# comparisons from sourmash
uniflna_test2_cmp <- 
  read.csv("/home/brian/Documents/data_mining/flu_project/uniflna_test2_cmp.csv")

# Label the rows
rownames(uniflna_test2_cmp) <- colnames(uniflna_test2_cmp)

# add gene column
uniflna_test2_wgs <- uniflna_test2_cmp
uniflna_test2_wgs$gene <- uniflna_test2$gene
uniflna_test2_wgs$type <- uniflna_test2$type
uniflna_test2_wgs$HNtype <- uniflna_test2$HNtype
#uniflna_test2_wgs$segment <- factor(uniflna_test1$segment)

# Transform for plotting
uniflna_test2_cmp_mat <- as.matrix(uniflna_test2_cmp)
```



```{r}
lda3 <- lda(gene ~ ., 
           uniflna_test2_wgs)

prop.lda3 = lda3$svd^2/sum(lda3$svd^2) 

plda3 <- predict(object = lda3,
                newdata = uniflna_test2_wgs)
```



```{r}
dataset3 = data.frame(gene = uniflna_test2_wgs[,2501],
                     lda3 = plda3$x)
```

```{r}
p3 <- ggplot(dataset3) + geom_point(aes(lda3.LD1, lda3.LD2, colour = gene), size=2, alpha=0.2) + 
  #theme_minimal() +3
  labs(x = paste("LD1 (", percent(prop.lda3[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda3[2]), ")", sep="")) +
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))
# https://stackoverflow.com/questions/5290003/how-to-set-legend-alpha-with-ggplot2

# ggsave("gene2500_LDA_plot.png", plot = p3, width = 7, height = 4, dpi = 200)

p3
```


ABC type labels 

```{r}
lda4 <- lda(type ~ ., 
           uniflna_test2_wgs)

prop.lda4 = lda4$svd^2/sum(lda4$svd^2) 

plda4 <- predict(object = lda4,
                newdata = uniflna_test2_wgs)
```



```{r}
dataset4 = data.frame(type = uniflna_test2_wgs[,2502],
                     lda4 = plda4$x)
```

```{r}
p4 <- ggplot(dataset4) + geom_point(aes(lda4.LD1, lda4.LD2, colour = type), size=2, alpha=0.2) + 
  #theme_minimal() +3
  labs(x = paste("LD1 (", percent(prop.lda4[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda4[2]), ")", sep="")) +
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))
# https://stackoverflow.com/questions/5290003/how-to-set-legend-alpha-with-ggplot2

#ggsave("ABCtype_LDA_plot.png", plot = p4, width = 7, height = 4, dpi = 200)

p4
```


```{r}
table(uniflna_test2_wgs$type)
```


HN type labels 

```{r}
lda5 <- lda(HNtype ~ ., 
           uniflna_test2_wgs)

prop.lda5 = lda5$svd^2/sum(lda5$svd^2) 

plda5 <- predict(object = lda5,
                newdata = uniflna_test2_wgs)
```



```{r}
dataset5 = data.frame(HNtype = uniflna_test2_wgs[,2503],
                     lda5 = plda5$x)
```

```{r}
p5 <- ggplot(dataset5) + geom_point(aes(lda5.LD1, lda5.LD2, colour = HNtype), size=2, alpha=0.5) + 
  #theme_minimal() +3
  labs(x = paste("LD1 (", percent(prop.lda5[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda5[2]), ")", sep="")) + theme(legend.position="none")
  #scale_color_manual(values=cpal) + 
  #guides(colour = guide_legend(override.aes = list(alpha = 1)))


p5
```

```{r}
table(uniflna_test2_wgs$HNtype)
```


```{r fig.width=12}
ggplot(dataset5) + geom_point(aes(lda5.LD1, lda5.LD2, colour = HNtype), size=2, alpha=0.5) + 
  #theme_minimal() +3
  labs(x = paste("LD1 (", percent(prop.lda5[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda5[2]), ")", sep=""))
```



```{r}
library(dplyr)
```

```{r}
fds5 <- dataset5 %>% group_by(HNtype) %>% filter(n() > 3) %>% filter(HNtype != 'none')
```



```{r fig.width=12}
ggplot(fds5) + geom_point(aes(lda5.LD1, lda5.LD2, colour = HNtype), size=2, alpha=0.4) + 
  #theme_minimal() +
  labs(x = paste("LD1 (", percent(prop.lda5[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda5[2]), ")", sep=""))
```






```{r}
ggplot(fds5) + geom_point(aes(lda5.LD1, lda5.LD2, colour = HNtype), size=2, alpha=0.4) + 
  labs(x = paste("LD1 (", percent(prop.lda5[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda5[2]), ")", sep="")) + 
  theme(legend.position="none")
  #scale_color_manual(values=cpal) + 
  #guides(colour = guide_legend(override.aes = list(alpha = 1)))
# https://stackoverflow.com/questions/5290003/how-to-set-legend-alpha-with-ggplot2

#ggsave("HNType_LDA_plot.png", width = 7, height = 4, dpi = 200)
```


```{r}
fds5p20 <- dataset5 %>% group_by(HNtype) %>% filter(n() >= 20) %>% filter(HNtype != 'none')
```

```{r}
tempcolors = c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f',
                '#ff7f00','#cab2d6','#6a3d9a','#ffff99','#b15928','#4d4d4d')
ggplot(fds5p20) + geom_point(aes(lda5.LD1, lda5.LD2, colour = HNtype), size=2, alpha=0.4) + 
  labs(x = paste("LD1 (", percent(prop.lda5[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda5[2]), ")", sep="")) + 
  scale_color_manual(values=tempcolors) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

# ggsave("HNType_thirty_LDA_plot.png", width = 7, height = 4, dpi = 200)
```



### Make an MDS plot

```{r mds}
fit <- dist(uniflna_test2_cmp_mat)
fit <- cmdscale(fit)
x <- fit[, 1]
y <- fit[, 2]
plot(fit[ , 1], fit[ , 2], xlab = "Dimension 1", ylab = "Dimension 2")
```


```{r}
mds_data <- data.frame(x=x, y=y)
mds_data$gene <- uniflna_test2$gene
mds_data$type <- uniflna_test2$type
mds_data$HNtype <- uniflna_test2$HNtype
```


```{r}
ggplot(mds_data, aes(x,y)) + geom_point(aes(color=gene), size=2, alpha=0.2) +
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

```

```{r}
ggplot(mds_data, aes(x,y)) + geom_point(aes(color=type), size=2, alpha=0.2) +
  scale_color_manual(values=cpal) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

```