Changes
Page history
Updated Practical 3 Processing 16S rRNA amplicon data (markdown)
authored
Feb 01, 2021
by
Ben Francis
Hide whitespace changes
Inline
Side-by-side
Practical-3-Processing-16S-rRNA-amplicon-data.md
View page @
98801c8c
...
@@ -135,40 +135,40 @@ library(reshape2)
...
@@ -135,40 +135,40 @@ library(reshape2)
library(RColorBrewer)
library(RColorBrewer)
plotTaxon <- function(seqtab, taxtable, taxon, min_abund, taxonomicLevel) {
plotTaxon <- function(seqtab, taxtable, taxon, min_abund, taxonomicLevel) {
seqtab <- data.frame(seqtab)
seqtab <- data.frame(seqtab)
colnames(seqtab) <- paste(taxtable[, 1], taxtable[, 2], taxtable[, 3], taxtable[, 4],
colnames(seqtab) <- paste(taxtable[, 1], taxtable[, 2], taxtable[, 3], taxtable[, 4],
taxtable[, 5], taxtable[, 6], as.character(1:length(taxtable[, 1])), sep="-")
taxtable[, 5], tax
tab
le
[,
6], as.character(1:length(taxtable[, 1])), sep="-")
seqtab <- seq
tab[,
grep("Chloroplast|Mitochondria", colnames(seqtab), invert=TRUE)]
seqtab <- seqtab
[, grep("Chloroplast|Mitochondria", colnames(seqtab), invert=TRUE)]
seqtab <- seqtab
/rowSums(seqtab)
seqtab <- seqtab
/rowSums(seqtab)
seqtab
1
<- seqtab
[, sapply(seqtab, function(x) max(x)) >= min_abund]
seqtab
1
<- seqtab[, sapply(seqtab, function(x) max(x))
>=
min_abund]
seqtab
2
<- seqtab[, sapply(seqtab, function(x) max(x))
<
min_abund]
seqtab2 <-
seqtab[, sapply(seqtab, function(x) max(x)) < min_abund]
colnames(
seqtab2
)
<-
paste0(colnames(seqtab2), "_Others")
colnames(
seqtab
2)
<-
paste0(colnames(
seqtab2
)
,
"_Other
s")
seqtab
3
<-
merge(seqtab1,
seqtab2,
by="row.name
s")
seqtab3 <- merge(seqtab1, seqtab2, by="
row.names
"
)
row.names(seqtab3) <-
row.names
(seqtab
)
row.names(
seqtab
3)
<-
row.names(
seqtab
)
seqtab <- seqtab
3
seqtab
<-
seqtab
3
seqtab
$sample <- rownames(
seqtab
)
seqtab
$sample <- rownames
(seqtab)
seqtab
.m <- melt
(seqtab)
seqtab.m <-
melt(seqtab)
seqtab.m <-
seqtab.m[grep(taxon, seqtab.m$variable),]
seqtab.m
<- seqtab.m[grep(taxon
, seqtab.m$variable)
,]
seqtab.m
$variable <- gsub("NA", "uncultured"
, seqtab.m$variable)
seqtab.m$variable <- gsub("
NA", "unculture
d", seqtab.m$variable)
seqtab.m$variable <- gsub("
.*Others", "Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abun
d", seqtab.m$variable)
seqtab.m$variable <-
gsub(".*Others", "Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund-Below_Min_Abund",
seqtab.m$variable
)
seqtab.m$variable <-
as.character(
seqtab.m$variable
seqtab.m$
variable <- as.character
(seqtab.m$variable
seqtab.m$
Domain <- sapply(strsplit
(seqtab.m$variable
, "-"), "[[", 1)
seqtab.m$
Domain
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
1
)
seqtab.m$
Phylum
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
2
)
seqtab.m$
Phylum
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
2
)
seqtab.m$
Class
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
3
)
seqtab.m$
Class
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
3
)
seqtab.m$
Order
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
4
)
seqtab.m$
Order
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
4
)
seqtab.m$
Family
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
5
)
seqtab.m$
Family
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
5
)
seqtab.m$
Genus
<- sapply(strsplit(seqtab.m$variable, "-"), "[[",
6
)
seqtab.m
$Genus <- sapply(strsplit(seqtab.m$variable, "-"), "[[", 6)
seqtab.m
<- cbind(seqtab.m, replicate(1,seqtab.m[colnames(seqtab.m) ==
seqtab.m <- cbind(seqtab.m, replicate(1,seqtab.m[colnames(seqtab.m) ==
taxonomicLevel]))
taxonomicLevel]))
colnames(seqtab.m) <- c(colnames(seqtab.m)[1:length(colnames(seqtab.m)) -1], "displayLevel")
colnames(seqtab.m) <- c(colnames(seqtab.m)[1:length(colnames(seqtab.m)) -1], "displayLevel")
nb.cols <- length(seqtab) - 1
nb.cols <- length(seqtab) - 1
mycolours <- rep(brewer.pal(12, "Paired"), ceiling(nb.cols/12))
mycolours <- rep(brewer.pal(12, "Paired"), ceiling(nb.cols/12))
ggplot(seqtab.m) +
ggplot(seqtab.m) +
geom_bar(aes(x=sample, y=value, fill=displayLevel), stat="identity", position="stack") +
geom_bar(aes(x=sample, y=value, fill=displayLevel), stat="identity", position="stack") +
scale_fill_manual(values=mycolours) +
scale_fill_manual(values=mycolours) +
labs(x="Sample", y="Proportion of Reads", fill=taxonomicLevel) +
labs(x="Sample", y="Proportion of Reads", fill=taxonomicLevel) +
theme_classic() +
theme_classic() +
theme(text=element_text(family="Serif", size=16)) +
theme(text=element_text(family="Serif", size=16)) +
theme(axis.text.x=element_text(angle=45, hjust=1))
theme(axis.text.x=element_text(angle=45, hjust=1))
}
}
```
```
...
...
...
...
