-
Notifications
You must be signed in to change notification settings - Fork 0
/
Lat_Normalisation_with_chloroplasts.Rmd
85 lines (76 loc) · 3.82 KB
/
Lat_Normalisation_with_chloroplasts.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
---
title: "Lat_Normalisation"
Date: 31.5.21
Correspondance: [email protected]
---
```{r Phylo object}
library(tidyverse)
library(dplyr)
setwd("~/Dropbox/UTS/PhD/Projects/Chapter1_Vibrio_latitude_study/Data/16S/Data")
sample_data <- read_csv('SMD.csv')
asv_cleaned_long <- read_csv("asv_cleaned_long_Chloroplasts.csv")
```
```{R Rarefaction with vegan - Filter to depth}
Read_Count <- asv_cleaned_long %>% group_by(SampleID) %>% summarise(Total = sum(Abundance))
asv_cleaned_long <- asv_cleaned_long %>% left_join(Read_Count)
Rareify.filt <- asv_cleaned_long %>% filter(Total > 30000)
Rareify.filt <- Rareify.filt %>% filter(!Abundance ==0)
Rareify.filt <- Rareify.filt %>% drop_na(SampleID)
```
```{r Rarefaction with vegan - Data Prep}
library(reshape2)
tableOFG = subset(Rareify.filt, select = c(ASV_name, ASV, SampleID, Abundance))
tableFGID = subset(tableOFG, select = c(ASV_name, SampleID, Abundance))
asvs.wide.OFG <- reshape2::dcast(tableFGID, SampleID~ASV_name, value.var = "Abundance", fill=0)
```
```{r Rareify}
rownames(asvs.wide.OFG)<-asvs.wide.OFG$SampleID
rarefy30k<-rrarefy(asvs.wide.OFG[,-c(1)], 30000)
rarefy30k.df<-as.data.frame(rarefy30k)
rarefy30k.df <- cbind(rownames(rarefy30k.df), data.frame(rarefy30k.df, row.names=NULL))
names(rarefy30k.df)[names(rarefy30k.df) == "rownames(rarefy30k.df)"] <- "ASV_name"
rarefy30k_long <- rarefy30k.df %>% pivot_longer(-ASV_name, names_to = "SampleID", values_to = "Abundance_Vegan_Rarefied")
rarefy30k_long <- rarefy30k_long %>% rename(SampleID. = ASV_name)
rarefy30k_long <- rarefy30k_long %>% rename(ASV_name = SampleID)
rarefy30k_long <- rarefy30k_long %>% rename(SampleID = SampleID.)
asv_rarefy30k_long <- asv_cleaned_long %>% left_join(rarefy30k_long)
Read_Count_rarefy30k <- asv_rarefy30k_long %>% group_by(SampleID) %>% summarise(Total = sum(Abundance_Vegan_Rarefied))
asv_rarefy30k_long <- asv_rarefy30k_long %>% group_by(SampleID) %>% mutate(RA = (Abundance_Vegan_Rarefied/30000)*100)
Read_Count_rarefy30k_RA <- asv_rarefy30k_long %>% group_by(SampleID) %>% summarise(Total = sum(RA))
asv_rarefy30k_long <- asv_rarefy30k_long %>% filter(!SampleID == "L88")
asv_rarefy30k_long <- asv_rarefy30k_long %>% filter(!SampleID == "L89")
asv_rarefy30k_long <- asv_rarefy30k_long %>% filter(!SampleID == "L90")
```
```{r write csv}
setwd("~/Dropbox/UTS/PhD/Projects/Chapter1_Vibrio_latitude_study/Data/16S/Data")
write_csv(asv_rarefy30k_long,'asv_rarefy30k_long_chloroplasts.csv')
```
```{r Rarefied30K NMDS Data preparation}
library(tidyverse)
library(vegan)
library(dplyr)
library(reshape2)
Rarefied30K_subset <- asv_rarefy30k_long %>% dplyr::select(ASV_name, Abundance_Vegan_Rarefied, SampleID)
Rarefied30K_subset <- Rarefied30K_subset %>% filter(!Abundance_Vegan_Rarefied ==0)
Rarefied30K.asvs.wide.OFG <- reshape2::dcast(Rarefied30K_subset , SampleID~ASV_name, value.var = "Abundance_Vegan_Rarefied", fill=0)
rownames(Rarefied30K.asvs.wide.OFG) <- Rarefied30K.asvs.wide.OFG$SampleID
Rarefied30K_nmds_df <- Rarefied30K.asvs.wide.OFG [,-c(1)]
```
```{r 30KRarefied NMDS}
Rarefied30K.nMDS <- metaMDS(Rarefied30K_nmds_df, distance = "bray", try =99, trymax=100, autotransform = F)
#stressplot(Rarefied30K.nMDS)
#plot(Rarefied30K.nMDS)
names(Rarefied30K.nMDS)
Rarefied30K.nMDS.points <- cbind(Rarefied30K.asvs.wide.OFG[,c(1)], as.data.frame(Rarefied30K.nMDS$points))
colnames(Rarefied30K.nMDS.points)[1] <- "SampleID"
Rarefied30K.nMDS.points <- Rarefied30K.nMDS.points %>% mutate(code=as.character(SampleID)) %>% left_join(sample_data)
Rarefied30K.nMDS.points$mycl450 <- as.character(Rarefied30K.nMDS.points$mycl450)
```
```{r 30KRarefied MDS}
library(tidyverse)
MDS.Rarefied30K <- ggplot(Rarefied30K.nMDS.points, aes(x=MDS1, y=MDS2, color = `mycl450`)) +
geom_point(aes(size=0.4, alpha=0.5)) +
theme_bw()
#scale_colour_manual(values=c("SubTropical"="#E66101", "Temperate"="#4DAC26" , "Tropical"="#0571B0"))
MDS.Rarefied30K
```