-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathRelative Abundance (1-Level, 1-Category).r
62 lines (54 loc) · 2.38 KB
/
Relative Abundance (1-Level, 1-Category).r
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
library(ggplot2)
library(stringr)
library(stringi)
# Import files
setwd("~/R/Analysis/1_Test")
DESIGN <- read.csv("experimental_design.csv",header=T)
setwd("~/R/Analysis/1_Test/ITS")
ASV.table <- read.table(file="rarefied_ASV_table.txt",header=T)
# % table
ASV <- ASV.table [,1:(ncol(ASV.table)-7)]
taxonomy <- ASV.table [,(ncol(ASV.table)-6):ncol(ASV.table)]
percent <- ASV / mean(colSums(ASV)) *100
# Remove "p__" before phylum name
taxonomy <- data.frame(lapply(taxonomy, function(x){gsub(pattern="p__", replacement = "", x)}),stringsAsFactors = FALSE) # Change p__ --> the first alphabet of the level you analyze
# Aggregate
agrregated <- aggregate(percent, by=list(taxonomy$Phylum),FUN = sum,na.rm=F) # Change Phylum --> the level you analyze
row.names(agrregated)<-agrregated[,1]
agrregated <- agrregated[,-1]
agrregated <- data.frame(agrregated)
rowMeans <- rowMeans(agrregated)
agrregated <- cbind(agrregated,rowMeans)
# Main + <1% abund
major <- agrregated[agrregated[,"rowMeans"] > 1,]
majors <- majors[order(majors$rowMeans,decreasing = T),]
minors <- agrregated[agrregated[,"rowMeans"] < 1,]
Others <- colSums (minors)
selected <- rbind (majors, others)
rownames (selected)[nrow(selected)] <- "Others"
selected <- selected[,-ncol(selected)]
# Make dataset
selected.t <- t (selected)
write.csv(selected.t, "aggregated.family.table.csv") #Change
bind <- cbind (selected.t, DESIGN)
table <- aggregate(selected.t, by=list(bind$Group),FUN = mean) # Change Group --> the category you compare
Group <- as.character(table[,1])
table <- data.frame(table)
data <- data.frame()
for (i in 2:(ncol(table))){
Abundance <-table[,i]
Example <- colnames (table)[i]
bind <- cbind(Group, Abundance, Example)
data<-rbind (data,bind)}
data$Abundance <- as.numeric(as.character(data$Abundance))
rownames<-colnames(table)[ncol(table):2]
data$Example <- factor(data$Example, levels = rownames)
# ggplot
ggplot (data, mapping=aes(x=Group, y=Abundance, fill=Example))+ # Change Group --> the category you compare
geom_bar(aes(), stat="identity", position="stack",color="black")+
# scale_fill_manual(values = c("gray","#C77CFF","#7CAE00","#00BFC4","#F8766D"))+ # if you want to change the colors
theme_classic()+
theme(text=element_text(size=14,color="black"),axis.text=element_text(size=12,color="black"))+
labs (y="Abundance (%)")
# Save
ggsave(file = "Rel.Abund.png")