kml/runVariation.R at main · Brodinlab/kml · GitHub

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
145
146
147
148
149
150
151
152
153
#!/usr/bin/env Rscript


########################
# SETTING UP LIBRARIES #
########################

library("argparse")
library("dplyr")
library("phyloseq")
library("vegan")
library("logr")


#######################
# SETTING UP LOGGING #
######################
timestamp <- format(Sys.time(), "%Y%m%d-%H%M%S")
logfile <- paste0("logfile_", timestamp, ".log")
lf = log_open(logfile)

options("logr.compact" = TRUE)

#######################
# SETTING UP ARGPARSER #
########################

sep("parsing cli args")

parser <- ArgumentParser(description = "script to run a permanova for trajectory assignments")
parser$add_argument("--input", type = "character", required = TRUE,
                    help = "Path to CLR tx rds.")
parser$add_argument("--metadata", type = "character", required = TRUE,
                    help = "Path metadata")
parser$add_argument("--subset", type = "character", required = TRUE,
                    help = "Naming variable only. Is this data a particular subset?")
parser$add_argument("--transformation", type = "character", required = TRUE,
                    help = "Naming variable only. Is this data transformed?")
parser$add_argument("--level", type = "character", required = TRUE,
                    help = "Naming variable only. What taxonomic level does the data represent?")

# Parse the command line arguments
args <- parser$parse_args()

# Access the input argument
input_arg <- args$input
metadata_arg <- args$metadata
tax_arg <- args$tax
subset_arg <- args$subset
transformation_arg <- args$transformation
level_arg <- args$level

options("logr.notes" = FALSE)
put(paste("user provided input: ", input_arg, sep = ""))
put(paste("user provided metadata: ", metadata_arg, sep = ""))
put(paste("user provided taxtable: ", tax_arg, sep = ""))
put(paste("user provided subset: ", subset_arg, sep = ""))
put(paste("user provided transformation: ", transformation_arg, sep = ""))
put(paste("user provided level: ", level_arg, sep = ""))
options("logr.notes" = TRUE)

#############################
# SETTING UP BASE FUNCTIONS #
#############################

check_dirs <- function(subdir_path) {

        if (!dir.exists(subdir_path)) {
        dir.create(subdir_path, recursive = TRUE)
        put(paste("Directory", subdir_path, "created."))
    } else {
        put(paste("Directory", subdir_path, "already exists."))
    }
}


#######################
# SETTING UP PATHINGS #
#######################

sep("preparations to iterate adonis2")

output_path = paste("output", "/", subset_arg, "/", transformation_arg, "/", level_arg, sep = "")
check_dirs(output_path)

#######################
# SETTING UP METADATA #
#######################

path_to_subject_assignments = paste(output_path, "/kmlSeedStatistics.rds", sep = "")
subject_assignments = readRDS(path_to_subject_assignments)[["taxaSubjectAssignments"]]

meged_subjects_assignments <- subject_assignments %>%
    lapply(function(df) {dplyr::select(df, -count, -prop)}) %>% # Remove columns
    purrr::reduce(dplyr::full_join, by = "subject") %>%
    distinct(subject, .keep_all = TRUE) # Merge on 'subject'


metadata = read.csv(metadata_arg, row.names = "X")
adonis2_md = merge(metadata, meged_subjects_assignments, by = "subject")
rownames(adonis2_md) = metadata$Sample


############################
# SETTING UP ADONIS2 INPUT #
############################

adonis2_input = readRDS(input_arg)

# running adonis2
sep("running adoni2")

adonis2_output = list()

for (level in c("species", "genus", "family", "order", "class")) {

    put(paste("working on level: ", level, sep = ""))

    adonis2_output[[level]] = list()

    for (taxa in names(subject_assignments)) {

        put(paste("working on taxa: ", taxa, sep = ""))

        adonis2_output[[level]][[taxa]] = list()

        taxa_columns = paste(taxa, "_clusters", sep = "")

        adonis2_taxa_md = adonis2_md %>% dplyr::rename(taxa_var = paste(taxa_columns))
        rownames(adonis2_taxa_md) = adonis2_taxa_md$Sample
        adonis2_taxa_md = adonis2_taxa_md[rownames(adonis2_input[[level]]), ]

        if (!identical(rownames(adonis2_taxa_md), rownames(adonis2_input[[level]]))) {

            put(taxa)
            put("DATAFRAMES ARE SORTED INACCURATELY")
            put("RESULTS ARE NOT REPRESENTATIVE")

        } else {

            adonis2_model = adonis2(adonis2_input[[level]] ~ taxa_var, data = adonis2_taxa_md, permutations = 999, method = "euclidean", by = "terms", strata = adonis2_taxa_md$subject)
            adonis2_model[["taxa"]] = taxa
            adonis2_model[["level"]] = level
            adonis2_model[["terms"]] = rownames(adonis2_model)
            adonis2_output[[level]][[taxa]] = adonis2_model

        }

    }

}

saveRDS(adonis2_output, paste(output_path, "/adonis2_output.rds", sep = ""))