forked from sydneyg/OTUvESV
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathTaxonomycorrelations_fungi_OTU_genusover1.R
More file actions
144 lines (107 loc) · 5.82 KB
/
Copy pathTaxonomycorrelations_fungi_OTU_genusover1.R
File metadata and controls
144 lines (107 loc) · 5.82 KB
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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
#June 22, 2018
#get dataframe for fungi OTU top genera over 1% relative abundance
#Reset R's Brain
rm(list=ls())
## Set working directory
#setwd("")
#set working directory
setwd("~/Dropbox/StatsandProgramming/16SElevationGradient/")
#load in desired libraries
library("ggplot2")
library("scales")
library("gridExtra")
library(plyr )
library(tidyverse)
library(stringr)
library(ggpubr)
library(vegan)
#source in functions
source('~/Dropbox/StatsandProgramming/source/gettaxondd.R', chdir = TRUE)
source('~/Dropbox/StatsandProgramming/source/getrowsums.R', chdir = TRUE)
################################################################################
############################## 1. DATA PROCESSING ##############################
################################################################################
# 1.1. Input the file of rarefied OTU-abundance-taxon table
#read in dataframes
#otu_abundance <- read.csv("data/fungi_ITS2/OTU_inoc_T1.csv", row.names=1)
otu_abundance_taxa <- read.csv("data/fungi_ITS2/otutab_inoculum_taxonomy_fungi.csv", row.names=1)
head(otu_abundance_taxa )
# 1.2. Extract otu-abundance table from OTU-abundance-taxon table
otu_abundance <- otu_abundance_taxa[, colnames(otu_abundance_taxa) != "Consensus.Lineage"]
head(otu_abundance)
#transform
otu_abundance_trans <- t(otu_abundance)
#rarefy
set.seed(10)
getrowsums(otu_abundance_trans)
library(vegan)
otu_abundance<- rrarefy(otu_abundance_trans,6665)
rowSums(otu_abundance)
#rarefied
otu_abundance <- t(otu_abundance)
colSums(otu_abundance)
#head(otu_taxa)
dd <- as.data.frame(otu_abundance_taxa$Consensus.Lineage)
dd$OTU <- row.names(otu_abundance_taxa)
names(dd) <- c("taxon","id")
# 1.3. Extract row of taxonomy from OTU-abundance-taxon table and modify the format
#otu_taxa <- otu_abundance_taxa["ConsensusLineage"] # Extract taxon table
library("stringr")
library("plyr")
otu_taxa1 <- ldply(str_split(string = dd$taxon, pattern=";"), rbind) # Divide a column using ";"and convert list to data frame
names(otu_taxa1) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
otu_taxa2 <- as.data.frame(lapply(otu_taxa1, gsub, pattern=" ", replacement=""))
otu_taxa3<- cbind(dd[,1:2 ],otu_taxa2)
#check column names of OTU table and taxonomy match
row.names(otu_abundance)==(otu_taxa3$id)
otu_abundance_taxa2 <- cbind(otu_abundance, otu_taxa3) #Combine otu ID and taxon table
head(otu_abundance_taxa2)
##########################################################################################
########################## 4. Community composition at phylum level ######################
##########################################################################################
# 4.1. Make dataframe for phylum composition
otu_abundance_taxa3 <- cbind(OTU_ID = rownames(otu_abundance_taxa2), otu_abundance_taxa2) #convert the rownames to a proper column of the data.frame
rownames(otu_abundance_taxa3) <- NULL
otu_abundance_taxa3$Phylum <- as.character(otu_abundance_taxa3$Phylum) #convert empty column to unidentified
otu_abundance_taxa3$Phylum[is.na(otu_abundance_taxa3$Phylum)]<-"p__unidentified"
otu_abundance_taxa3$Family <- as.character(otu_abundance_taxa3$Family) #convert empty column to unidentified
otu_abundance_taxa3$Family[is.na(otu_abundance_taxa3$Family)]<-"f__unidentified"
otu_abundance_taxa3$Family[otu_abundance_taxa3$Family=="f__"]<-"f__unidentified"
otu_abundance_taxa3$Genus <- as.character(otu_abundance_taxa3$Genus) #convert empty column to unidentified
otu_abundance_taxa3$Genus[is.na(otu_abundance_taxa3$Genus)]<-"g__unidentified"
otu_abundance_taxa3$Genus[otu_abundance_taxa3$Genus=="g__"]<-"g__unidentified"
otu_abundance_taxa3$Kingdom_Phylum <- str_c(otu_abundance_taxa3$Phylum,otu_abundance_taxa3$Family, otu_abundance_taxa3$Genus, sep=";") #combine kingdom and phylun
otu_abundance_taxa4 <- otu_abundance_taxa3[, !(colnames(otu_abundance_taxa3) %in% c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"))]
otu_abundance_taxa4 <- otu_abundance_taxa4[, -c(22,23)]
head(otu_abundance_taxa4)
# 4.2. Caluculate mean value for each phylum on each date
library("reshape2")
otu_abundance_taxa5 <- melt(otu_abundance_taxa4, id.vars= c("OTU_ID","Kingdom_Phylum"), variable.name="Sample", value.name="Abundance")
otu_abundance_taxa5 <- subset(otu_abundance_taxa5, Abundance>0)
#get total number of sequences per sample
totalseqs <- colSums(otu_abundance_taxa4[ ,2:21])
avgnumseqs <- sum(totalseqs)/20
otu_abundance_taxa5$Abundance <- (otu_abundance_taxa5$Abundance/avgnumseqs)*100 # Check your sequence number and convert to relative abundance (%)
head(otu_abundance_taxa5)
phylum_abundance <- ddply(otu_abundance_taxa5, c("Kingdom_Phylum","Sample"), summarise,
sum = sum(Abundance, na.rm=TRUE)
)
head(phylum_abundance)
phylum_abundance$Site <- str_sub(phylum_abundance$Sample,2,2) #get a column of site names
head(phylum_abundance)
phylum_abundance2 <- ddply(phylum_abundance, c("Kingdom_Phylum","Site"), summarise, # Caluculate mean value on each date
mean = mean(sum, na.rm=TRUE),
sd = sd(sum, na.rm=TRUE),
n = sum(!is.na(sum)),
se = sd/sqrt(n),
max = max(sum, na.rm=TRUE)
)
head(phylum_abundance2)
# 4.3. Classify phylum with mean abundance < XX% (1% here) into "other phylum".
# Twelve phylum and others are recommended to use palette="Paired".
# "maximum abundance" rather than "mean abundance" will be better when compositional change is large.
phylum_over1 <- subset(phylum_abundance2, mean>1) # Subset the phylum with >1%
#if one site has more than 1% of that phylum then it stays
phylum_over1_vector <- as.vector(unique(phylum_over1$Kingdom_Phylum)) # Subset the phylum with >1%
length(phylum_over1_vector) #34
write.csv(phylum_over1, "Figures/otuvzotu_taxonomy/fungi_OTU_genusover1%relativeabundance.csv")