/
LjAt_NC_load_data.R
executable file
·61 lines (39 loc) · 1.52 KB
/
LjAt_NC_load_data.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
# originally by Ruben Garrido-Oter
# garridoo@mpipz.mpg.de
# load libraries
library("ggplot2")
library("scales")
library("grid")
library("vegan")
# load plotting functions
source("plotting_functions.R")
source("plotting_parameters.R")
source("cpcoa.func.R")
source("paths.R")
# files
design.file <- paste(data.dir, "LjAt_NC_design.txt", sep="")
otu_table.file <- paste(results.dir, "LjAt_NC_ASV_table.txt", sep="")
taxonomy.file <- paste(results.dir, "LjAt_NC_taxonomy.txt", sep="")
# load data
design <- read.table(design.file, header=T, sep="\t")
otu_table <- read.table(otu_table.file, sep="\t", header=T, check.names=F)
taxonomy <- read.table(taxonomy.file, sep="\t", header=T, check.names=F)
# treat the order Burkholderiales (formerly Proteobacteriales) as Proteobacteria (Class)
taxonomy$Class <- as.character(taxonomy$Class)
idx <- taxonomy$Order=="Burkholderiales"
taxonomy$Class[idx] <- "Betaproteobacteria"
# rarefy ASV table
otu_table_raref <- rrarefy(otu_table, sample=1000)
# re-order data matrices
idx <- design$Original.SampleID %in% colnames(otu_table)
design <- design[idx, ]
idx <- match(design$Original.SampleID, colnames(otu_table))
otu_table <- otu_table[, idx]
otu_table_raref <- otu_table_raref[, idx]
colnames(otu_table) <- design$SampleID
# remove chloroplast reads
chloro <- taxonomy$ASV[taxonomy$Phylum=="Cyanobacteria/Chloroplast"]
otu_table <- otu_table[!rownames(otu_table) %in% chloro, ]
# normalize otu table
design$depth <- colSums(otu_table)
otu_table <- apply(otu_table, 2, function(x) x/sum(x))