|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("adegenet"))
|
|
|
4 suppressPackageStartupMessages(library("ape"))
|
|
|
5 suppressPackageStartupMessages(library("data.table"))
|
|
|
6 suppressPackageStartupMessages(library("dbplyr"))
|
|
|
7 suppressPackageStartupMessages(library("dplyr"))
|
|
|
8 suppressPackageStartupMessages(library("ggplot2"))
|
|
|
9 suppressPackageStartupMessages(library("knitr"))
|
|
|
10 suppressPackageStartupMessages(library("maps"))
|
|
|
11 suppressPackageStartupMessages(library("mapproj"))
|
|
|
12 suppressPackageStartupMessages(library("optparse"))
|
|
|
13 suppressPackageStartupMessages(library("poppr"))
|
|
|
14 suppressPackageStartupMessages(library("RColorBrewer"))
|
|
|
15 suppressPackageStartupMessages(library("RPostgres"))
|
|
|
16 suppressPackageStartupMessages(library("SNPRelate"))
|
|
|
17 suppressPackageStartupMessages(library("tidyr"))
|
|
|
18 suppressPackageStartupMessages(library("vcfR"))
|
|
|
19 suppressPackageStartupMessages(library("vegan"))
|
|
|
20 suppressPackageStartupMessages(library("yarrr"))
|
|
|
21 theme_set(theme_bw())
|
|
|
22
|
|
|
23 DEFAULT_MISSING_NUMERIC_VALUE <- -9.000000;
|
|
|
24
|
|
|
25 option_list <- list(
|
|
|
26 make_option(c("--database_connection_string"), action="store", dest="database_connection_string", help="Corals (stag) database connection string"),
|
|
|
27 make_option(c("--input_affy_metadata"), action="store", dest="input_affy_metadata", help="Affymetrix 96 well plate input file"),
|
|
|
28 make_option(c("--input_pop_info"), action="store", dest="input_pop_info", help="Population information input file"),
|
|
|
29 make_option(c("--input_vcf"), action="store", dest="input_vcf", help="VCF input file"),
|
|
|
30 make_option(c("--output_nj_phylogeny_tree"), action="store", dest="output_nj_phylogeny_tree", default=NULL, help="Flag to plot neighbor-joining phylogeny tree"),
|
|
|
31 make_option(c("--output_stag_db_report"), action="store", dest="output_stag_db_report", help="Flag to output stag db report file")
|
|
|
32 )
|
|
|
33
|
|
|
34 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
|
|
|
35 args <- parse_args(parser, positional_arguments=TRUE);
|
|
|
36 opt <- args$options;
|
|
|
37
|
|
|
38 get_file_path = function(dir, file_name) {
|
|
|
39 file_path = paste(dir, file_name, sep="/");
|
|
|
40 return(file_path);
|
|
|
41 }
|
|
|
42
|
|
|
43 get_database_connection <- function(db_conn_string) {
|
|
|
44 # Instantiate database connection.
|
|
|
45 # The connection string has this format:
|
|
|
46 # postgresql://user:password@host/dbname
|
|
|
47 conn_items <- strsplit(db_conn_string, "://")[[1]];
|
|
|
48 string_needed <- conn_items[2];
|
|
|
49 items_needed <- strsplit(string_needed, "@")[[1]];
|
|
|
50 user_pass_string <- items_needed[1];
|
|
|
51 host_dbname_string <- items_needed[2];
|
|
|
52 user_pass_items <- strsplit(user_pass_string, ":")[[1]];
|
|
|
53 host_dbname_items <- strsplit(host_dbname_string, "/")[[1]];
|
|
|
54 user <- user_pass_items[1];
|
|
|
55 pass <- user_pass_items[2];
|
|
|
56 host <- host_dbname_items[1];
|
|
|
57 dbname <- host_dbname_items[2];
|
|
|
58 conn <- DBI::dbConnect(RPostgres::Postgres(), host=host, port="5432", dbname=dbname, user=user, password=pass);
|
|
|
59 return (conn);
|
|
|
60 }
|
|
|
61
|
|
4
|
62 log_data_frame <- function(name, df) {
|
|
|
63 cat("\n", name, ":\n");
|
|
|
64 show(df);
|
|
|
65 cat("\n\n");
|
|
|
66 }
|
|
|
67
|
|
0
|
68 time_elapsed <- function(start_time) {
|
|
|
69 cat("Elapsed time: ", proc.time() - start_time, "\n\n");
|
|
|
70 }
|
|
|
71
|
|
|
72 time_start <- function(msg) {
|
|
|
73 start_time <- proc.time();
|
|
|
74 cat(msg, "...\n");
|
|
|
75 return(start_time);
|
|
|
76 }
|
|
|
77
|
|
|
78 write_data_frame <- function(dir, file_name, data_frame) {
|
|
|
79 cat("\nWriting file: ", file_name, "\n");
|
|
|
80 file_path <- get_file_path(dir, file_name);
|
|
|
81 write.table(data_frame, file=file_path, quote=FALSE, row.names=FALSE, sep="\t");
|
|
|
82 }
|
|
|
83
|
|
|
84 # Prepare for processing.
|
|
|
85 output_data_dir = "output_data_dir";
|
|
|
86 output_plots_dir = "output_plots_dir";
|
|
|
87 # Read in VCF input file.
|
|
|
88 start_time <- time_start("Reading VCF input");
|
|
|
89 vcf <- read.vcfR(opt$input_vcf);
|
|
|
90 time_elapsed(start_time);
|
|
|
91
|
|
|
92 # Convert VCF file into a genind for the Poppr package.
|
|
|
93 start_time <- time_start("Converting VCF data to a genind object");
|
|
|
94 genind_obj <- vcfR2genind(vcf);
|
|
4
|
95 cat("\ngenind_obj:\n");
|
|
|
96 genind_obj
|
|
|
97 cat("\n\n");
|
|
0
|
98 time_elapsed(start_time);
|
|
|
99
|
|
|
100 # Add population information to the genind object.
|
|
|
101 population_info_data_table <- read.table(opt$input_pop_info, check.names=FALSE, header=F, na.strings=c("", "NA"), stringsAsFactors=FALSE, sep="\t", quote="");
|
|
|
102 colnames(population_info_data_table) <- c("row_id", "affy_id", "user_specimen_id", "region");
|
|
4
|
103 cat("\npopulation_info_data_table:\n");
|
|
|
104 population_info_data_table
|
|
|
105 cat("\n\n");
|
|
0
|
106 #write_data_frame(output_data_dir, "population_info_data_table", population_info_data_table);
|
|
|
107 genind_obj@pop <- as.factor(population_info_data_table$region);
|
|
|
108 strata(genind_obj) <- data.frame(pop(genind_obj));
|
|
|
109
|
|
|
110 # Convert genind object to a genclone object.
|
|
|
111 start_time <- time_start("Converting the genind object to a genclone object");
|
|
|
112 genind_clone <- as.genclone(genind_obj);
|
|
4
|
113 cat("\ngenind_clone:\n");
|
|
|
114 genind_clone
|
|
|
115 cat("\n\n");
|
|
0
|
116 time_elapsed(start_time);
|
|
7
|
117 # Remove genind object from memory.
|
|
|
118 rm(genind_obj);
|
|
0
|
119
|
|
|
120 # Calculate the bitwise distance between individuals.
|
|
|
121 start_time <- time_start("Calculating the bitwise distance between individuals");
|
|
|
122 bitwise_distance <- bitwise.dist(genind_clone);
|
|
|
123 time_elapsed(start_time);
|
|
|
124
|
|
|
125 # Multilocus genotypes (threshold of 3.2%).
|
|
4
|
126 cat("\nFiltering multilocus genotypes with threshold of 3.2%...\n\n");
|
|
0
|
127 mlg.filter(genind_clone, distance=bitwise_distance) <- 0.032;
|
|
|
128
|
|
|
129 # Create list of MLGs.
|
|
4
|
130 cat("\nCreating list of mlg_ids...\n\n");
|
|
0
|
131 mlg_ids <- mlg.id(genind_clone);
|
|
8
|
132 cat("\nCreated list of mlg_ids...\n\n");
|
|
0
|
133
|
|
|
134 # Read user's Affymetrix 96 well plate tabular file.
|
|
8
|
135 cat("\nCreating affy_metadata_data_frame...\n\n");
|
|
0
|
136 affy_metadata_data_frame <- read.table(opt$input_affy_metadata, header=FALSE, stringsAsFactors=FALSE, sep="\t", na.strings=c("", "NA"), quote="");
|
|
8
|
137 cat("\nCreated affy_metadata_data_frame...\n\n");
|
|
0
|
138 colnames(affy_metadata_data_frame) <- c("user_specimen_id", "field_call", "bcoral_genet_id", "bsym_genet_id", "reef",
|
|
|
139 "region", "latitude", "longitude", "geographic_origin", "colony_location",
|
|
|
140 "depth", "disease_resist", "bleach_resist", "mortality","tle",
|
|
|
141 "spawning", "collector_last_name", "collector_first_name", "organization", "collection_date",
|
|
|
142 "email", "seq_facility", "array_version", "public", "public_after_date",
|
|
|
143 "sperm_motility", "healing_time", "dna_extraction_method", "dna_concentration", "registry_id",
|
|
|
144 "result_folder_name", "plate_barcode");
|
|
|
145 affy_metadata_data_frame$user_specimen_id <- as.character(affy_metadata_data_frame$user_specimen_id);
|
|
4
|
146 log_data_frame("affy_metadata_data_frame", affy_metadata_data_frame);
|
|
0
|
147 user_specimen_ids <- as.character(affy_metadata_data_frame$user_specimen_id);
|
|
4
|
148 cat("\nuser_specimen_ids:\n", toString(user_specimen_ids), "\n\n");
|
|
0
|
149 # The specimen_id_field_call_data_table looks like this:
|
|
|
150 # user_specimen_ids V2
|
|
|
151 # 1090 prolifera
|
|
|
152 # 1091 prolifera
|
|
4
|
153 cat("\nCreating specimen_id_field_call_data_table...\n");
|
|
0
|
154 specimen_id_field_call_data_table <- data.table(user_specimen_ids, affy_metadata_data_frame$field_call);
|
|
|
155 # Rename the user_specimen_ids column.
|
|
|
156 setnames(specimen_id_field_call_data_table, c("user_specimen_ids"), c("user_specimen_id"));
|
|
|
157 # Rename the V2 column.
|
|
|
158 setnames(specimen_id_field_call_data_table, c("V2"), c("field_call"));
|
|
|
159
|
|
|
160 # Connect to database.
|
|
|
161 conn <- get_database_connection(opt$database_connection_string);
|
|
|
162 # Import the sample table.
|
|
|
163 sample_table <- tbl(conn, "sample");
|
|
|
164 # Import the genotype table.
|
|
|
165 genotype_table <- tbl(conn, "genotype");
|
|
2
|
166 # Import the probe_annotation table.
|
|
|
167 probe_annotation_table <- tbl(conn, "probe_annotation");
|
|
0
|
168 # Select columns from the sample table and the
|
|
|
169 # genotype table joined by genotype_id.
|
|
|
170 sample_table_columns <- sample_table %>% select(user_specimen_id, affy_id, bcoral_genet_id, genotype_id);
|
|
|
171 smlg <- sample_table_columns %>%
|
|
|
172 left_join(genotype_table %>%
|
|
|
173 select("id", "coral_mlg_clonal_id", "coral_mlg_rep_sample_id", "genetic_coral_species_call"),
|
|
|
174 by=c("genotype_id"="id"));
|
|
|
175 # Name the columns.
|
|
|
176 smlg_data_frame <- as.data.frame(smlg);
|
|
|
177 colnames(smlg_data_frame) <- c("user_specimen_id", "affy_id", "bcoral_genet_id", "genotype_id",
|
|
9
|
178 "coral_mlg_clonal_id", "coral_mlg_rep_sample_id", "genetic_coral_species_call");
|
|
4
|
179 log_data_frame("smlg_data_frame", smlg_data_frame);
|
|
0
|
180 # Missing GT in samples submitted.
|
|
|
181 start_time <- time_start("Discovering missing GT in samples");
|
|
|
182 gt <- extract.gt(vcf, element="GT", as.numeric=FALSE);
|
|
|
183 missing_gt <- apply(gt, MARGIN=2, function(x){ sum(is.na(x))});
|
|
|
184 missing_gt <- (missing_gt / nrow(vcf)) * 100;
|
|
|
185 missing_gt_data_frame <- data.frame(missing_gt);
|
|
4
|
186 log_data_frame("missing_gt_data_frame", missing_gt_data_frame);
|
|
0
|
187 # The specimen_id_field_call_data_table looks like this:
|
|
|
188 # rn missing_gt
|
|
|
189 # a100000-4368120-060520-256_I07.CEL 0.06092608
|
|
|
190 # a100000-4368120-060520-256_K07.CEL 0.05077173
|
|
4
|
191 cat("\nCreating missing_gt_data_table...\n");
|
|
|
192 missing_gt_data_table <- setDT(missing_gt_data_frame, keep.rownames=TRUE)[];
|
|
0
|
193 # Rename the rn column.
|
|
|
194 setnames(missing_gt_data_table, c("rn"), c("affy_id"));
|
|
|
195 # Rename the missing_gt column.
|
|
|
196 setnames(missing_gt_data_table, c("missing_gt"), c("percent_missing_data_coral"));
|
|
|
197 # Round data to two digits.
|
|
|
198 missing_gt_data_table$percent_missing_data_coral <- round(missing_gt_data_table$percent_missing_data_coral, digits=2);
|
|
|
199 time_elapsed(start_time);
|
|
|
200
|
|
2
|
201 # Subset genotypes for the fixed SNPs by probe id.
|
|
|
202 # Select columns from the probe_annotation table.
|
|
|
203 probe_annotation_table_columns <- probe_annotation_table %>% select(probe_set_id, custid, fixed_status, acerv_allele);
|
|
|
204 # Convert to data frame.
|
|
|
205 fixed_snp_data_frame <- as.data.frame(probe_annotation_table_columns);
|
|
|
206 # Name the columns.
|
|
|
207 colnames(fixed_snp_data_frame) <- c("probe_set_id", "custid", "fixed_status", "acerv_allele");
|
|
|
208 # Filter unwanted rows.
|
|
|
209 fixed_snp_data_frame <- subset(fixed_snp_data_frame, fixed_snp_data_frame$fixed_status=="keep");
|
|
4
|
210 log_data_frame("fixed_snp_data_frame", fixed_snp_data_frame);
|
|
2
|
211 gt_fixed <- gt[rownames(gt) %in% fixed_snp_data_frame$probe_set_id, ];
|
|
4
|
212 log_data_frame("gt_fixed", gt_fixed);
|
|
2
|
213
|
|
|
214 # Missing GT in fixed SNPs.
|
|
|
215 missing_gt_fixed <- apply(gt_fixed, MARGIN=2, function(x){ sum(is.na(x))});
|
|
|
216 missing_gt_fixed <- (missing_gt_fixed / nrow(gt_fixed)) * 100;
|
|
|
217 missing_gt_fixed_data_frame <- data.frame(missing_gt_fixed);
|
|
4
|
218 log_data_frame("missing_gt_fixed_data_frame", missing_gt_fixed_data_frame);
|
|
|
219 cat("\nCreating missing_gt_fixed_data_table...\n");
|
|
2
|
220 missing_gt_fixed_data_table <- setDT(missing_gt_fixed_data_frame, keep.rownames=TRUE)[];
|
|
|
221 # Rename the rn column.
|
|
|
222 setnames(missing_gt_fixed_data_table, c("rn"), c("affy_id"));
|
|
|
223 # Rename the missing_gt column.
|
|
|
224 setnames(missing_gt_fixed_data_table, c("missing_gt_fixed"), c("percent_missing_data_fixed"));
|
|
|
225 # Round data to two digits.
|
|
|
226 missing_gt_fixed_data_table$percent_missing_data_fixed <- round(missing_gt_fixed_data_table$percent_missing_data_fixed, digits=2);
|
|
|
227
|
|
|
228 # Heterozygous alleles for fixed SNPs.
|
|
0
|
229 start_time <- time_start("Discovering heterozygous alleles");
|
|
2
|
230 heterozygous_alleles <- apply(gt_fixed, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("0/1", x))))});
|
|
|
231 heterozygous_alleles <- (heterozygous_alleles / nrow(gt_fixed)) * 100;
|
|
0
|
232 heterozygous_alleles_data_frame <- data.frame(heterozygous_alleles);
|
|
4
|
233 log_data_frame("heterozygous_alleles_data_frame", heterozygous_alleles_data_frame);
|
|
0
|
234 # The heterozygous_alleles_data_table looks like this:
|
|
|
235 # rn heterozygous_alleles
|
|
|
236 # a100000-4368120-060520-256_I07.CEL 73.94903
|
|
|
237 # a100000-4368120-060520-256_K07.CEL 74.40089
|
|
4
|
238 cat("\nCreating heterozygous_alleles_data_table...\n");
|
|
0
|
239 heterozygous_alleles_data_table <- setDT(heterozygous_alleles_data_frame, keep.rownames=TRUE)[];
|
|
|
240 # Rename the rn column.
|
|
|
241 setnames(heterozygous_alleles_data_table, c("rn"), c("affy_id"));
|
|
|
242 # Rename the heterozygous_alleles column.
|
|
|
243 setnames(heterozygous_alleles_data_table, c("heterozygous_alleles"), c("percent_heterozygous_coral"));
|
|
|
244 # Round data to two digits.
|
|
|
245 heterozygous_alleles_data_table$percent_heterozygous_coral <- round(heterozygous_alleles_data_table$percent_heterozygous_coral, digits=2);
|
|
|
246 time_elapsed(start_time);
|
|
|
247
|
|
2
|
248 # Create list of Acerv reference and alternative probes.
|
|
|
249 rAC <- subset(fixed_snp_data_frame, fixed_snp_data_frame$acerv_allele=="reference");
|
|
|
250 aAC <- subset(fixed_snp_data_frame, fixed_snp_data_frame$acerv_allele=="alternative");
|
|
|
251
|
|
|
252 # Subset probes for the reference and alternative SNPs in Acerv.
|
|
|
253 ref_ac <- gt_fixed[rownames(gt_fixed) %in% rAC$probe,];
|
|
|
254 alt_ac <- gt_fixed[rownames(gt_fixed) %in% aAC$probe,];
|
|
|
255
|
|
|
256 # Reference alleles for each species.
|
|
|
257 reference_alleles_ac <- apply(ref_ac, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("0/0", x))))});
|
|
|
258 reference_alleles_ap <- apply(alt_ac, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("0/0", x))))});
|
|
|
259
|
|
|
260 # Alternative alleles for each species.
|
|
|
261 alternative_alleles_ac <- apply(alt_ac, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("1/1", x))))});
|
|
|
262 alternative_alleles_ap <- apply(ref_ac, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("1/1", x))))});
|
|
|
263
|
|
|
264 # Apalm alleles.
|
|
0
|
265 start_time <- time_start("Discovering reference alleles");
|
|
2
|
266 ap_sum <- rowSums(cbind(reference_alleles_ap,alternative_alleles_ap));
|
|
|
267 ap_alleles <- (ap_sum / nrow(gt_fixed)) * 100;
|
|
|
268 ap_alleles_data_frame <- data.frame(ap_alleles);
|
|
4
|
269 log_data_frame("ap_alleles_data_frame", ap_alleles_data_frame);
|
|
0
|
270 # The reference_alleles_data_table looks like this:
|
|
|
271 # rn reference_alleles
|
|
|
272 # a100000-4368120-060520-256_I07.CEL 11.60642
|
|
|
273 # a100000-4368120-060520-256_K07.CEL 11.45918
|
|
4
|
274 cat("\nCreating ap_alleles_data_table...\n");
|
|
2
|
275 ap_alleles_data_table <- setDT(ap_alleles_data_frame, keep.rownames=TRUE)[];
|
|
0
|
276 # Rename the rn column.
|
|
2
|
277 setnames(ap_alleles_data_table, c("rn"), c("affy_id"));
|
|
0
|
278 # Rename the reference_alleles column.
|
|
2
|
279 setnames(ap_alleles_data_table, c("ap_alleles"), c("percent_apalm_coral"));
|
|
0
|
280 # Round data to two digits.
|
|
2
|
281 ap_alleles_data_table$percent_apalm_coral <- round(ap_alleles_data_table$percent_apalm_coral, digits=2);
|
|
0
|
282 time_elapsed(start_time);
|
|
|
283
|
|
2
|
284 # Acerv alleles.
|
|
0
|
285 start_time <- time_start("Discovering alternative alleles");
|
|
2
|
286 ac_sum <- rowSums(cbind(reference_alleles_ac,alternative_alleles_ac));
|
|
|
287 ac_alleles <- (ac_sum / nrow(gt_fixed)) * 100;
|
|
|
288 ac_alleles_data_frame <- data.frame(ac_alleles);
|
|
4
|
289 log_data_frame("ac_alleles_data_frame", ac_alleles_data_frame);
|
|
0
|
290 # The alternative_alleles_data_table looks like this:
|
|
|
291 # rn alternative_alleles
|
|
|
292 # a100000-4368120-060520-256_I07.CEL 14.38363
|
|
|
293 # a100000-4368120-060520-256_K07.CEL 14.08916
|
|
4
|
294 cat("\nCreating ac_alleles_data_table...\n");
|
|
2
|
295 ac_alleles_data_table <- setDT(ac_alleles_data_frame, keep.rownames=TRUE)[];
|
|
0
|
296 # Rename the rn column.
|
|
2
|
297 setnames(ac_alleles_data_table, c("rn"), c("affy_id"));
|
|
0
|
298 # Rename the alternative_alleles column.
|
|
2
|
299 setnames(ac_alleles_data_table, c("ac_alleles"), c("percent_acerv_coral"));
|
|
0
|
300 # Round data to two digits.
|
|
2
|
301 ac_alleles_data_table$percent_acerv_coral <- round(ac_alleles_data_table$percent_acerv_coral, digits=2);
|
|
0
|
302 time_elapsed(start_time);
|
|
|
303
|
|
|
304 # The mlg_ids_data_table looks like this:
|
|
|
305 # mlg_ids
|
|
|
306 # a550962-4368120-060520-500_M23.CEL
|
|
|
307 # a550962-4368120-060520-256_A19.CEL
|
|
4
|
308 cat("\nCreating mlg_ids_data_table...\n");
|
|
0
|
309 mlg_ids_data_table <- data.table(mlg_ids, keep.rownames=TRUE);
|
|
|
310 # Rename the mlg_ids column.
|
|
|
311 setnames(mlg_ids_data_table, c("mlg_ids"), c("affy_id"));
|
|
|
312
|
|
|
313 # sample_mlg_tibble looks like this:
|
|
|
314 # A tibble: 262 x 3
|
|
|
315 # Groups: group [?]
|
|
|
316 # group affy_id coral_mlg_clonal_id coral_mlg_rep_sample_id
|
|
|
317 # <int> <chr> <chr> <chr>
|
|
|
318 # 1 a550962-4368.CEL NA 13905
|
|
|
319 sample_mlg_tibble <- mlg_ids_data_table %>%
|
|
|
320 group_by(row_number()) %>%
|
|
|
321 dplyr::rename(group="row_number()") %>%
|
|
|
322 unnest (affy_id) %>%
|
|
|
323 # Join with mlg table.
|
|
|
324 left_join(smlg_data_frame %>%
|
|
8
|
325 select("affy_id","coral_mlg_clonal_id", "coral_mlg_rep_sample_id",
|
|
|
326 "genetic_coral_species_call", "bcoral_genet_id"),
|
|
0
|
327 by="affy_id");
|
|
|
328
|
|
|
329 # If found in database, group members on previous mlg id.
|
|
|
330 uniques <- unique(sample_mlg_tibble[c("group", "coral_mlg_clonal_id")]);
|
|
|
331 uniques <- uniques[!is.na(uniques$coral_mlg_clonal_id),];
|
|
|
332 na.mlg <- which(is.na(sample_mlg_tibble$coral_mlg_clonal_id));
|
|
|
333 na.group <- sample_mlg_tibble$group[na.mlg];
|
|
|
334 sample_mlg_tibble$coral_mlg_clonal_id[na.mlg] <- uniques$coral_mlg_clonal_id[match(na.group, uniques$group)];
|
|
|
335
|
|
|
336 # Find out if the sample mlg matched a previous genotyped sample.
|
|
|
337 # sample_mlg_match_tibble looks like this:
|
|
|
338 # A tibble: 262 x 4
|
|
|
339 # Groups: group [230]
|
|
|
340 # group affy_id coral_mlg_clonal_id db_match
|
|
|
341 # <int> <chr> <chr> <chr>
|
|
|
342 # 1 a550962-436.CEL NA no_match
|
|
|
343 sample_mlg_match_tibble <- sample_mlg_tibble %>%
|
|
|
344 group_by(group) %>%
|
|
|
345 mutate(db_match = ifelse(is.na(coral_mlg_clonal_id), "no_match", "match"));
|
|
|
346
|
|
|
347 # Create new mlg id for samples with no matches in the database.
|
|
|
348 none <- unique(sample_mlg_match_tibble[c("group", "coral_mlg_clonal_id")]);
|
|
|
349 none <- none[is.na(none$coral_mlg_clonal_id),];
|
|
|
350 na.mlg2 <- which(is.na(sample_mlg_match_tibble$coral_mlg_clonal_id));
|
|
|
351 n.g <- sample_mlg_match_tibble$group[na.mlg2];
|
|
|
352 ct <- length(unique(n.g));
|
|
|
353
|
|
|
354 # List of new group ids, the sequence starts at the number of
|
|
|
355 # ids present in sample_mlg_match_tibble$coral_mlg_clonal_ids
|
|
|
356 # plus 1.
|
|
|
357 n.g_ids <- sprintf("HG%04d", seq((sum(!is.na(unique(sample_mlg_match_tibble["coral_mlg_clonal_id"]))) + 1), by=1, length=ct));
|
|
|
358
|
|
|
359 # Assign the new id iteratively for all that have NA.
|
|
|
360 for (i in 1:length(na.mlg2)) {
|
|
|
361 sample_mlg_match_tibble$coral_mlg_clonal_id[na.mlg2[i]] <- n.g_ids[match(sample_mlg_match_tibble$group[na.mlg2[i]], unique(n.g))];
|
|
|
362 }
|
|
|
363
|
|
|
364 # Subset population_info_data_table for all samples.
|
|
|
365 # affy_id_user_specimen_id_vector looks like this:
|
|
|
366 # affy_id user_specimen_id
|
|
|
367 # a100000-432.CEL 13704
|
|
|
368 affy_id_user_specimen_id_vector <- population_info_data_table[c(2, 3)];
|
|
|
369
|
|
|
370 # Merge data frames for final table.
|
|
|
371 start_time <- time_start("Merging data frames");
|
|
|
372 stag_db_report <- specimen_id_field_call_data_table %>%
|
|
|
373 left_join(affy_id_user_specimen_id_vector %>%
|
|
|
374 select("affy_id", "user_specimen_id"),
|
|
|
375 by="user_specimen_id") %>%
|
|
|
376 mutate(db_record = ifelse(affy_id %in% smlg_data_frame$affy_id, "genotyped", "new")) %>%
|
|
|
377 filter(db_record=="new") %>%
|
|
|
378 left_join(sample_mlg_match_tibble %>%
|
|
|
379 select("affy_id", "coral_mlg_clonal_id", "db_match"),
|
|
|
380 by="affy_id") %>%
|
|
|
381 left_join(missing_gt_data_table %>%
|
|
|
382 select("affy_id", "percent_missing_data_coral"),
|
|
|
383 by="affy_id") %>%
|
|
2
|
384 left_join(missing_gt_fixed_data_table %>%
|
|
|
385 select("affy_id", "percent_missing_data_fixed"),
|
|
|
386 by="affy_id") %>%
|
|
0
|
387 left_join(heterozygous_alleles_data_table %>%
|
|
|
388 select("affy_id", "percent_heterozygous_coral"),
|
|
|
389 by="affy_id") %>%
|
|
2
|
390 left_join(ac_alleles_data_table %>%
|
|
1
|
391 select("affy_id", "percent_acerv_coral"),
|
|
0
|
392 by="affy_id") %>%
|
|
2
|
393 left_join(ap_alleles_data_table %>%
|
|
1
|
394 select("affy_id", "percent_apalm_coral"),
|
|
0
|
395 by="affy_id") %>%
|
|
|
396 mutate(db_match = ifelse(is.na(db_match), "failed", db_match))%>%
|
|
|
397 mutate(coral_mlg_clonal_id = ifelse(is.na(coral_mlg_clonal_id), "failed", coral_mlg_clonal_id)) %>%
|
|
2
|
398 mutate(genetic_coral_species_call = ifelse(percent_apalm_coral >= 85 & percent_acerv_coral <= 10, "A.palmata", "other")) %>%
|
|
|
399 mutate(genetic_coral_species_call = ifelse(percent_acerv_coral >= 85 & percent_apalm_coral <= 10, "A.cervicornis", genetic_coral_species_call)) %>%
|
|
0
|
400 mutate(genetic_coral_species_call = ifelse(percent_heterozygous_coral > 40, "A.prolifera", genetic_coral_species_call)) %>%
|
|
|
401 ungroup() %>%
|
|
2
|
402 select(-group, -db_record);
|
|
0
|
403 time_elapsed(start_time);
|
|
|
404
|
|
|
405 start_time <- time_start("Writing csv output");
|
|
|
406 write.csv(stag_db_report, file=opt$output_stag_db_report, quote=FALSE);
|
|
|
407 time_elapsed(start_time);
|
|
|
408
|
|
|
409 # Representative clone for genotype table.
|
|
|
410 start_time <- time_start("Creating representative clone for genotype table");
|
|
2
|
411 no_dup_genotypes_genind <- clonecorrect(genind_clone, strata=~pop.genind_obj.);
|
|
0
|
412 id_rep <- mlg.id(no_dup_genotypes_genind);
|
|
4
|
413 cat("\nCreating id_data_table...\n");
|
|
0
|
414 id_data_table <- data.table(id_rep, keep.rownames=TRUE);
|
|
|
415 # Rename the id_rep column.
|
|
|
416 setnames(id_data_table, c("id_rep"), c("affy_id"));
|
|
|
417 time_elapsed(start_time);
|
|
7
|
418 # Remove clonecorrect genind from memory.
|
|
|
419 rm(no_dup_genotypes_genind);
|
|
0
|
420
|
|
|
421 # Table of alleles for the new samples subset to new plate data.
|
|
|
422 # Create vector indicating number of individuals desired from
|
|
|
423 # affy_id column of stag_db_report data table.
|
|
|
424 i <- ifelse(is.na(stag_db_report[3]), "", stag_db_report[[3]]);
|
|
2
|
425 i <- i[!apply(i== "", 1, all), ];
|
|
0
|
426
|
|
|
427 # Subset VCF to the user samples.
|
|
|
428 start_time <- time_start("Subsetting vcf to the user samples");
|
|
6
|
429 affy_list <- append(stag_db_report$affy_id,"FORMAT");
|
|
|
430 svcf <- vcf[,colnames(vcf@gt) %in% affy_list];
|
|
0
|
431 write.vcf(svcf, "subset.vcf.gz");
|
|
7
|
432
|
|
|
433 # Remove original and subset VCFs written to file from R memory.
|
|
|
434 rm(svcf);
|
|
|
435 rm(vcf);
|
|
|
436
|
|
|
437 # Load in subset VCF.
|
|
0
|
438 vcf.fn <- "subset.vcf.gz";
|
|
|
439 snpgdsVCF2GDS(vcf.fn, "test3.gds", method="biallelic.only");
|
|
|
440 genofile <- snpgdsOpen(filename="test3.gds", readonly=FALSE);
|
|
|
441 gds_array <- read.gdsn(index.gdsn(genofile, "sample.id"));
|
|
|
442 # gds_array looks like this:
|
|
|
443 # [1] "a550962-4368120-060520-500_A03.CEL" "a550962-4368120-060520-500_A05.CEL"
|
|
|
444 gds_data_frame <- data.frame(gds_array);
|
|
4
|
445 log_data_frame("gds_data_frame", gds_data_frame);
|
|
0
|
446 # gds_data_frame looks like this:
|
|
|
447 # gds_array
|
|
|
448 # a550962-4368120-060520-500_A03.CEL
|
|
|
449 # a550962-4368120-060520-500_A05.CEL
|
|
4
|
450 cat("\nCreating gds_data_table...\n");
|
|
0
|
451 gds_data_table <- setDT(gds_data_frame, keep.rownames=FALSE)[];
|
|
|
452 # Rename the gds_array column.
|
|
|
453 setnames(gds_data_table, c("gds_array"), c("affy_id"));
|
|
|
454 # affy_id_region_list looks like this:
|
|
|
455 # affy_id region
|
|
|
456 # a100000-4368120-060520-256_I07.CEL USVI
|
|
|
457 # a100000-4368120-060520-256_K07.CEL USVI
|
|
|
458 affy_id_region_list <- population_info_data_table[c(2,3,4)];
|
|
|
459 gds_data_table_join <- gds_data_table %>%
|
|
|
460 left_join(affy_id_region_list %>%
|
|
2
|
461 select("affy_id", "user_specimen_id", "region"),
|
|
0
|
462 by='affy_id')%>%
|
|
|
463 drop_na();
|
|
|
464 samp.annot <- data.frame(pop.group=c(gds_data_table_join$region));
|
|
|
465 add.gdsn(genofile, "sample.annot", samp.annot);
|
|
|
466 # population_code looks like this:
|
|
|
467 # [1] 18.361733 18.361733 18.361733 18.361733 18.361733 18.361733
|
|
|
468 # [7] 25.11844009 25.11844009 25.11844009 25.11844009 25.11844009 25.11844009
|
|
|
469 population_code <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"));
|
|
|
470 pop.group <- as.factor(read.gdsn(index.gdsn(genofile, "sample.annot/pop.group")));
|
|
|
471 # pop.group looks like this:
|
|
|
472 # [1] 18.361733 18.361733 18.361733 18.361733 18.361733 18.361733
|
|
|
473 # [7] 25.11844009 25.11844009 25.11844009 25.11844009 25.11844009 25.11844009
|
|
|
474 time_elapsed(start_time);
|
|
|
475
|
|
|
476 # Distance matrix calculation and sample labels change to user specimen ids.
|
|
|
477 start_time <- time_start("Calculating distance matrix");
|
|
|
478 ibs <- snpgdsIBS(genofile, num.thread=2, autosome.only=FALSE);
|
|
|
479 ibs$sample.id <-gds_data_table_join$user_specimen_id;
|
|
|
480 time_elapsed(start_time);
|
|
|
481
|
|
|
482 # Cluster analysis on the genome-wide IBS pairwise distance matrix.
|
|
|
483 start_time <- time_start("Clustering the genome-wide IBS pairwise distance matrix");
|
|
|
484 set.seed(100);
|
|
2
|
485 par(cex=0.6, cex.lab=1, cex.axis=1.5, cex.main=2);
|
|
0
|
486 ibs.hc <- snpgdsHCluster(ibs);
|
|
|
487 time_elapsed(start_time);
|
|
|
488
|
|
|
489 # cols looks like this:
|
|
|
490 # blue1 red green pink orange blue2
|
|
|
491 # "#0C5BB0FF" "#EE0011FF" "#15983DFF" "#EC579AFF" "#FA6B09FF" "#149BEDFF"
|
|
|
492 # green2 yellow turquoise poop
|
|
|
493 # "#A1C720FF" "#FEC10BFF" "#16A08CFF" "#9A703EFF"
|
|
|
494 cols <- piratepal("basel");
|
|
|
495 set.seed(999);
|
|
|
496
|
|
|
497 # Generate plots.
|
|
|
498 # Default clustering.
|
|
|
499 start_time <- time_start("Creating ibs_default.pdf");
|
|
|
500 # Start PDF device driver.
|
|
|
501 dev.new(width=40, height=20);
|
|
|
502 file_path = get_file_path(output_plots_dir, "ibs_default.pdf");
|
|
|
503 pdf(file=file_path, width=40, height=20);
|
|
|
504 rv <- snpgdsCutTree(ibs.hc, col.list=cols, pch.list=15);
|
|
|
505 snpgdsDrawTree(rv, main="Color by Cluster", leaflab="perpendicular", yaxis.kinship=FALSE);
|
|
2
|
506 abline(h=0.032, lty=2);
|
|
0
|
507 legend("topleft", legend=levels(rv$samp.group), xpd=T, col=cols[1:nlevels(rv$samp.group)], pch=15, ncol=4, cex=1.2);
|
|
|
508 dev.off()
|
|
|
509 time_elapsed(start_time);
|
|
|
510
|
|
|
511 # Color cluster by region.
|
|
|
512 start_time <- time_start("Creating ibs_region.pdf");
|
|
|
513 # Start PDF device driver.
|
|
|
514 dev.new(width=40, height=20);
|
|
|
515 file_path = get_file_path(output_plots_dir, "ibs_region.pdf");
|
|
|
516 pdf(file=file_path, width=40, height=20);
|
|
|
517 race <- as.factor(population_code);
|
|
2
|
518 rv2 <- snpgdsCutTree(ibs.hc, samp.group=race, col.list=cols, pch.list=15);
|
|
0
|
519 snpgdsDrawTree(rv2, main="Color by Region", leaflab="perpendicular", yaxis.kinship=FALSE);
|
|
|
520 legend("topleft", legend=levels(race), xpd=T, col=cols[1:nlevels(race)], pch=15, ncol=4, cex=1.2);
|
|
|
521 dev.off()
|
|
|
522 time_elapsed(start_time);
|
|
|
523
|
|
|
524 # Missing data barplot.
|
|
|
525 start_time <- time_start("Creating missing_data.pdf");
|
|
|
526 population_info_data_table$miss <- stag_db_report$percent_missing_data_coral[match(missing_gt_data_frame$affy_id, stag_db_report$affy_id)];
|
|
|
527 test2 <- which(!is.na(population_info_data_table$miss));
|
|
|
528 miss96 <- population_info_data_table$miss[test2];
|
|
|
529 name96 <- population_info_data_table$user_specimen_id[test2];
|
|
|
530 # Start PDF device driver.
|
|
|
531 dev.new(width=20, height=10);
|
|
|
532 file_path = get_file_path(output_plots_dir, "missing_data.pdf");
|
|
|
533 pdf(file=file_path, width=20, height=10);
|
|
|
534 par(mar = c(8, 4, 4, 2));
|
|
|
535 x <- barplot(miss96, las=2, col=cols, ylim=c(0, 3), cex.axis=0.8, space=0.8, ylab="Missingness (%)", xaxt="n");
|
|
|
536 text(cex=0.8, x=x-0.25, y=-.05, name96, xpd=TRUE, srt=60, adj=1);
|
|
|
537 dev.off()
|
|
|
538 time_elapsed(start_time);
|
|
|
539
|
|
|
540 # Sample MLG on a map.
|
|
|
541 start_time <- time_start("Creating mlg_map.pdf");
|
|
|
542 # Get the lattitude and longitude boundaries for rendering
|
|
|
543 # the map. Tese boundaries will restrict the map to focus
|
|
|
544 # (i.e., zoom) on the region of the world map from which
|
|
|
545 # the samples were taken.
|
|
|
546 max_latitude <- max(affy_metadata_data_frame$latitude, na.rm=TRUE);
|
|
|
547 min_latitude <- min(affy_metadata_data_frame$latitude, na.rm=TRUE);
|
|
|
548 latitude_range_vector <- c(min_latitude-3, max_latitude+3);
|
|
|
549 max_longitude <- max(affy_metadata_data_frame$longitude, na.rm=TRUE);
|
|
|
550 min_longitude <- min(affy_metadata_data_frame$longitude, na.rm=TRUE);
|
|
|
551 longitude_range_vector <- c(min_longitude-3, max_longitude+3);
|
|
|
552 # Get the palette colors for rendering plots.
|
|
|
553 colors <- length(unique(stag_db_report$coral_mlg_clonal_id));
|
|
|
554 # Get a color palette.
|
|
|
555 palette <- colorRampPalette(piratepal("basel"));
|
|
|
556 # Start PDF device driver.
|
|
|
557 dev.new(width=20, height=20);
|
|
|
558 file_path = get_file_path(output_plots_dir, "mlg_map.pdf");
|
|
|
559 pdf(file=file_path, width=20, height=20);
|
|
|
560 world_data = map_data("world");
|
|
|
561 # Add the coral_mlg_clonal_id column from the stag_db_report
|
|
|
562 # data fram to the affy_metadata_data_frame.
|
|
|
563 affy_metadata_data_frame$mlg <- stag_db_report$coral_mlg_clonal_id;
|
|
|
564 # Get the number of colors needed from the palette for plotting
|
|
|
565 # the sample locations on the world map.
|
|
|
566 num_colors = length(unique(affy_metadata_data_frame$mlg));
|
|
|
567 # Get a color palette.
|
|
|
568 palette = colorRampPalette(piratepal("basel"));
|
|
|
569 ggplot() +
|
|
|
570 geom_map(data=world_data, map=world_data, aes(x=long, y=lat, group=group, map_id=region), fill="white", colour="#7f7f7f") +
|
|
|
571 coord_quickmap(xlim=longitude_range_vector, ylim=latitude_range_vector) +
|
|
|
572 geom_point(data=affy_metadata_data_frame, aes(x=longitude, y=latitude, group=mlg, colour=mlg), alpha=.7, size=3) +
|
|
|
573 scale_color_manual(values=palette(num_colors)) +
|
|
|
574 theme(legend.position="bottom") +
|
|
|
575 guides(color=guide_legend(nrow=8, byrow=F));
|
|
|
576
|
|
|
577 # Sample MLG on a map for each region.
|
|
|
578 for (i in unique(affy_metadata_data_frame$region)) {
|
|
2
|
579 m <- i;
|
|
|
580 num_colors_2 = length(unique(affy_metadata_data_frame$mlg[which(affy_metadata_data_frame$region == m)]));
|
|
|
581 max_latitude_region <- max(affy_metadata_data_frame$latitude[which(affy_metadata_data_frame$region == m)], na.rm=TRUE);
|
|
|
582 min_latitude_region <- min(affy_metadata_data_frame$latitude[which(affy_metadata_data_frame$region == m)], na.rm=TRUE);
|
|
|
583 latitude_range_vector_region <- c(min_latitude_region-0.5, max_latitude_region+0.5);
|
|
|
584 max_longitude_region <- max(affy_metadata_data_frame$longitude[which(affy_metadata_data_frame$region == m)], na.rm=TRUE);
|
|
|
585 min_longitude_region <- min(affy_metadata_data_frame$longitude[which(affy_metadata_data_frame$region == m)], na.rm=TRUE);
|
|
|
586 longitude_range_vector_region <- c(min_longitude_region-0.5, max_longitude_region+0.5);
|
|
|
587 print(ggplot() +
|
|
|
588 geom_map(data=world_data, map=world_data, aes(x=long, y=lat, group=group, map_id=region),
|
|
|
589 fill="grey", colour="#7f7f7f") +
|
|
|
590 coord_quickmap(xlim=longitude_range_vector_region, ylim=latitude_range_vector_region, clip = "on") +
|
|
|
591 geom_point(data=affy_metadata_data_frame[which(affy_metadata_data_frame$region == m),],
|
|
|
592 aes(x=longitude, y=latitude, group=mlg, colour=mlg), alpha=.5, size=3) +
|
|
|
593 scale_color_manual(values=palette(num_colors_2)) +
|
|
|
594 theme(legend.position="bottom") + labs(title=paste("MLG assignments for", m)) +
|
|
|
595 guides(color=guide_legend(nrow=8, byrow=F)));
|
|
0
|
596 }
|
|
|
597 dev.off()
|
|
|
598 time_elapsed(start_time);
|
|
|
599
|
|
|
600 if (!is.null(opt$output_nj_phylogeny_tree)) {
|
|
|
601 # Create a phylogeny tree of samples based on distance matrices.
|
|
|
602 # Start PDF device driver.
|
|
|
603 start_time <- time_start("Creating nj_phylogeny_tree.pdf");
|
|
|
604 # Table of alleles for the new samples subset to new plate data.
|
|
|
605 # Create vector indicating number of individuals desired from
|
|
|
606 # affy_id column of stag_db_report data table.
|
|
|
607 i <- ifelse(is.na(stag_db_report[1]), "", stag_db_report[[1]]);
|
|
|
608 i <- i[!apply(i== "", 1, all),];
|
|
|
609 sample_alleles_vector <- genind_clone[i, mlg.reset=FALSE, drop=FALSE];
|
|
|
610 dev.new(width=40, height=80);
|
|
|
611 file_path = get_file_path(output_plots_dir, "nj_phylogeny_tree.pdf");
|
|
|
612 pdf(file=file_path, width=40, height=80);
|
|
|
613 # Organize branches by clade.
|
|
|
614 nj_phylogeny_tree <- sample_alleles_vector %>%
|
|
|
615 aboot(dist=provesti.dist, sample=100, tree="nj", cutoff=50, quiet=TRUE, showtree = FALSE) %>%
|
|
|
616 ladderize();
|
|
|
617 nj_phylogeny_tree$tip.label <- stag_db_report$user_specimen_id[match(nj_phylogeny_tree$tip.label, stag_db_report$affy_id)];
|
|
|
618 plot.phylo(nj_phylogeny_tree, tip.color=cols[sample_alleles_vector$pop], label.offset=0.0025, cex=0.6, font=2, lwd=4, align.tip.label=F, no.margin=T);
|
|
|
619 # Add a scale bar showing 5% difference.
|
|
|
620 add.scale.bar(0, 0.95, length=0.05, cex=0.65, lwd=2);
|
|
|
621 nodelabels(nj_phylogeny_tree$node.label, cex=.5, adj=c(1.5, -0.1), frame="n", font=3, xpd=TRUE);
|
|
|
622 legend("topright", legend=c(levels(sample_alleles_vector$pop)), text.col=cols, xpd=T, cex=0.8);
|
|
|
623 dev.off()
|
|
|
624 time_elapsed(start_time);
|
|
|
625 }
|
|
|
626
|
|
|
627 # Generate a pie chart for each sample with a genotype.
|
|
|
628 # Store the numerical and user_specimen_id values from
|
|
|
629 # stag_db_report for the charts (user_specimen_id names
|
|
|
630 # will be used to label each chart).
|
|
|
631 start_time <- time_start("Creating percent_breakdown.pdf");
|
|
3
|
632 stag_db_report_data_table <- stag_db_report[c(-2, -3, -4, -5)];
|
|
0
|
633 # Remove NA and NaN values.
|
|
|
634 stag_db_report_data_table <- na.omit(stag_db_report_data_table);
|
|
|
635 # Translate to N (i.e., number of samples with a genotype)
|
|
|
636 # columns and 5 rows.
|
|
|
637 translated_stag_db_report_data_table <- t(stag_db_report_data_table);
|
|
|
638 translated_stag_db_report_matrix <- as.matrix(translated_stag_db_report_data_table[-1,]);
|
|
|
639 # Set the storage mode of the matrix to numeric. In some
|
|
|
640 # cases this could result in the following:
|
|
|
641 # Warning message:
|
|
|
642 # In mde(x) : NAs introduced by coercion
|
|
|
643 mode(translated_stag_db_report_matrix) <- "numeric";
|
|
|
644 # Remove NA and NaN values that may have been introduced
|
|
|
645 # by coercion.
|
|
|
646 translated_stag_db_report_matrix <- na.omit(translated_stag_db_report_matrix);
|
|
2
|
647 translated_stag_db_report_matrix<-translated_stag_db_report_matrix[-c(1),];
|
|
|
648 #tsdbrm_row_means <- rowMeans(translated_stag_db_report_matrix, na.rm=TRUE);
|
|
0
|
649 dev.new(width=10, height=7);
|
|
|
650 file_path = get_file_path(output_plots_dir, "percent_breakdown.pdf");
|
|
|
651 pdf(file=file_path, width=10, height=7);
|
|
|
652 par(mfrow=c(3, 2));
|
|
2
|
653 col <- c("#A6A6A6","#FFA626","#EB0ACF", "#80FF00" );
|
|
0
|
654 # Generate a pie chart for each sample with genotypes.
|
|
|
655 for (i in 1:ncol(translated_stag_db_report_matrix)) {
|
|
2
|
656 tmp_labels <- paste(c("no call", "heterozygous", "A. cervicornis", "A. palmata"), " (", round(translated_stag_db_report_matrix[,i], 1), "%)", sep="");
|
|
0
|
657 main <- paste("Breakdown of SNP assignments for", translated_stag_db_report_data_table[1, i]);
|
|
|
658 pie(translated_stag_db_report_matrix[,i], labels=tmp_labels, radius=0.90, col=col, main=main, cex.main=.85, cex=0.75);
|
|
|
659 }
|
|
|
660 dev.off()
|
|
|
661 time_elapsed(start_time);
|
|
|
662
|
|
|
663 # close GDS file.
|
|
|
664 snpgdsClose(genofile);
|
|
|
665
|
|
|
666 # Prepare to output data frames for input to a downstream
|
|
|
667 # tool that will use them to update the stag database.
|
|
|
668 start_time <- time_start("Building data frames for insertion into database tables");
|
|
|
669 # sample_prep_data_frame looks like this (split across comment lines):
|
|
|
670 # user_specimen_id field_call bcoral_genet_id bsym_genet_id reef
|
|
|
671 # test_002 prolifera NA NA JohnsonsReef
|
|
|
672 # region latitude longitude geographic_origin colony_location
|
|
|
673 # Bahamas 18.36173 -64.77430 Reef NA
|
|
|
674 # depth disease_resist bleach_resist
|
|
|
675 # 5 NA N
|
|
|
676 # mortality tle spawning collector_last_name collector_first_name organization
|
|
|
677 # NA NA False Kitchen Sheila Penn State
|
|
|
678 # collection_date email seq_facility array_version public
|
|
|
679 # 2018-11-08 k89@psu.edu Affymetrix 1 True
|
|
|
680 # public_after_date sperm_motility healing_time dna_extraction_method
|
|
|
681 # NA -9 -9 NA
|
|
|
682 # dna_concentration registry_id result_folder_name plate_barcode mlg
|
|
|
683 # NA NA PRO100175_PSU175_SAX_b02 P9SR10074 HG0227
|
|
|
684 # affy_id percent_missing_data_coral percent_heterozygous_coral
|
|
|
685 # a550962-436.CEL 1.06 19.10
|
|
1
|
686 # percent_acerv_coral percent_apalm_coral
|
|
0
|
687 # 40.10459 39.73396
|
|
|
688 sample_prep_data_frame <- affy_metadata_data_frame %>%
|
|
|
689 left_join(stag_db_report %>%
|
|
|
690 select("user_specimen_id", "affy_id", "percent_missing_data_coral", "percent_heterozygous_coral",
|
|
1
|
691 "percent_acerv_coral", "percent_apalm_coral"),
|
|
0
|
692 by='user_specimen_id');
|
|
|
693 # Get the number of rows for all data frames.
|
|
|
694 num_rows <- nrow(sample_prep_data_frame);
|
|
|
695 # Set the column names so that we can extract only those columns
|
|
|
696 # needed for the sample table.
|
|
|
697 colnames(sample_prep_data_frame) <- c("user_specimen_id", "field_call", "bcoral_genet_id", "bsym_genet_id", "reef",
|
|
|
698 "region", "latitude", "longitude", "geographic_origin", "colony_location",
|
|
|
699 "depth", "disease_resist", "bleach_resist", "mortality", "tle",
|
|
|
700 "spawning", "collector_last_name", "collector_first_name", "organization",
|
|
|
701 "collection_date", "email", "seq_facility", "array_version", "public",
|
|
|
702 "public_after_date", "sperm_motility", "healing_time", "dna_extraction_method",
|
|
|
703 "dna_concentration", "registry_id", "result_folder_name", "plate_barcode",
|
|
|
704 "mlg", "affy_id", "percent_missing_data_coral", "percent_heterozygous_coral",
|
|
1
|
705 "percent_acerv_coral", "percent_apalm_coral");
|
|
0
|
706
|
|
|
707 # Output the data frame for updating the alleles table.
|
|
|
708 # Subset to only the new plate data.
|
|
|
709 i <- ifelse(is.na(stag_db_report[3]), "", stag_db_report[[3]]);
|
|
|
710 # Create a vector indicating the number of individuals desired
|
|
|
711 # from the affy_id collumn in the report_user data table.
|
|
|
712 i <- i[!apply(i=="", 1, all),];
|
|
|
713 # Subset the genclone object to the user data.
|
|
|
714 allele_vector <- genind_clone[i, mlg.reset=FALSE, drop=FALSE];
|
|
|
715 # Convert the subset genclone to a data frame.
|
|
|
716 allele_data_frame <- genind2df(allele_vector, sep="");
|
|
|
717 allele_data_frame <- allele_data_frame %>%
|
|
|
718 select(-pop);
|
|
|
719 # Allele string for Allele.table in database.
|
|
|
720 allele_table_data_frame <- unite(allele_data_frame, alleles, 1:19696, sep=" ", remove=TRUE);
|
|
|
721 allele_table_data_frame <- setDT(allele_table_data_frame, keep.rownames=TRUE)[];
|
|
|
722 setnames(allele_table_data_frame, c("rn"), c("affy_id"));
|
|
|
723 # write.csv(concat_sample_alleles,file=paste("Seed_genotype_alleles.csv",sep = ""),quote=FALSE,row.names=FALSE);
|
|
|
724 write_data_frame(output_data_dir, "allele.tabular", allele_table_data_frame);
|
|
|
725
|
|
|
726 # Output the data frame for updating the experiment table.
|
|
|
727 experiment_table_data_frame <- data.frame(matrix(ncol=4, nrow=num_rows));
|
|
|
728 colnames(experiment_table_data_frame) <- c("seq_facility", "array_version", "result_folder_name", "plate_barcode");
|
|
|
729 for (i in 1:num_rows) {
|
|
|
730 experiment_table_data_frame$seq_facility[i] <- sample_prep_data_frame$seq_facility[i];
|
|
|
731 experiment_table_data_frame$array_version[i] <- sample_prep_data_frame$array_version[i];
|
|
|
732 experiment_table_data_frame$result_folder_name[i] <- sample_prep_data_frame$result_folder_name[i];
|
|
|
733 experiment_table_data_frame$plate_barcode[i] <- sample_prep_data_frame$plate_barcode[i];
|
|
|
734 }
|
|
|
735 write_data_frame(output_data_dir, "experiment.tabular", experiment_table_data_frame);
|
|
|
736
|
|
|
737 # Output the data frame for updating the colony table.
|
|
|
738 # The geographic_origin value is used for deciding into which table
|
|
|
739 # to insert the latitude and longitude values. If the geographic_origin
|
|
|
740 # is "reef", the values will be inserted into the reef table, and if it is
|
|
|
741 # "colony", the values will be inserted into the colony table. We insert
|
|
|
742 # these values in both data frames so that the downstream tool that parses
|
|
|
743 # them can determine the appropriate table.
|
|
|
744 colony_table_data_frame <- data.frame(matrix(ncol=4, nrow=num_rows));
|
|
|
745 colnames(colony_table_data_frame) <- c("latitude", "longitude", "depth", "geographic_origin");
|
|
|
746 for (i in 1:num_rows) {
|
|
|
747 colony_table_data_frame$latitude[i] <- sample_prep_data_frame$latitude[i];
|
|
|
748 colony_table_data_frame$longitude[i] <- sample_prep_data_frame$longitude[i];
|
|
|
749 colony_table_data_frame$depth[i] <- sample_prep_data_frame$depth[i];
|
|
|
750 colony_table_data_frame$geographic_origin[i] <- sample_prep_data_frame$geographic_origin[i];
|
|
|
751 }
|
|
|
752 write_data_frame(output_data_dir, "colony.tabular", colony_table_data_frame);
|
|
|
753
|
|
|
754 # Output the data frame for populating the genotype table.
|
|
|
755 # Combine with previously genotyped samples.
|
|
|
756 # prep_genotype_tibble looks like this:
|
|
|
757 # A tibble: 220 x 7
|
|
|
758 # Groups: group [?]
|
|
|
759 # group affy_id coral_mlg_clona… user_specimen_id db_match
|
|
|
760 # <int> <chr> <chr> <chr> <chr>
|
|
|
761 # 1 a10000… 13905 HG0048 match
|
|
|
762 # genetic_coral_species_call coral_mlg_rep_sample_id
|
|
|
763 # <chr> <chr>
|
|
|
764 # A.palmata 1104
|
|
|
765 prep_genotype_tibble <- id_data_table %>%
|
|
|
766 group_by(row_number()) %>%
|
|
|
767 dplyr::rename(group='row_number()') %>%
|
|
|
768 unnest(affy_id) %>%
|
|
8
|
769 left_join(sample_mlg_match_tibble %>%
|
|
|
770 select("affy_id", "coral_mlg_rep_sample_id", "coral_mlg_clonal_id",
|
|
|
771 "genetic_coral_species_call", "bcoral_genet_id", "db_match"),
|
|
|
772 by='affy_id') %>%
|
|
|
773 right_join(sample_mlg_match_tibble %>%
|
|
|
774 select("coral_mlg_rep_sample_id", "coral_mlg_clonal_id"),
|
|
|
775 by='coral_mlg_clonal_id') %>%
|
|
|
776 mutate(coral_mlg_rep_sample_id=ifelse(is.na(coral_mlg_rep_sample_id.x),coral_mlg_rep_sample_id.y,coral_mlg_rep_sample_id.x)) %>%
|
|
|
777 ungroup() %>%
|
|
|
778 dplyr::select(-coral_mlg_rep_sample_id.x,-coral_mlg_rep_sample_id.y, -group.x,-group.y) %>%
|
|
9
|
779 group_by(coral_mlg_clonal_id) %>%
|
|
|
780 arrange(coral_mlg_rep_sample_id) %>%
|
|
|
781 slice(1);
|
|
8
|
782
|
|
0
|
783 # Confirm that the representative mlg is the same between runs.
|
|
|
784 uniques2 <- unique(prep_genotype_tibble[c("group", "coral_mlg_rep_sample_id")]);
|
|
|
785 uniques2 <- uniques2[!is.na(uniques2$coral_mlg_rep_sample_id),];
|
|
|
786 na.mlg3 <- which(is.na(prep_genotype_tibble$coral_mlg_rep_sample_id));
|
|
|
787 na.group2 <- prep_genotype_tibble$group[na.mlg3];
|
|
|
788 prep_genotype_tibble$coral_mlg_rep_sample_id[na.mlg3] <- uniques2$coral_mlg_rep_sample_id[match(na.group2, uniques2$group)];
|
|
|
789 # Transform the representative mlg column with new genotyped samples.
|
|
|
790 # representative_mlg_tibble looks like this:
|
|
|
791 # A tibble: 220 x 5
|
|
|
792 # affy_id coral_mlg_rep_sa… coral_mlg_clona… user_specimen_id
|
|
|
793 # <chr> <chr> <chr> <chr>
|
|
|
794 # a100000-… 13905 HG0048 13905
|
|
|
795 # genetic_coral_species_call bcoral_genet_id
|
|
|
796 # <chr> <chr>
|
|
|
797 # A.palmata C1651
|
|
|
798 representative_mlg_tibble <- prep_genotype_tibble %>%
|
|
8
|
799 mutate(coral_mlg_rep_sample_id=ifelse(is.na(coral_mlg_rep_sample_id) & (db_match =="no_match"), affy_id, coral_mlg_rep_sample_id)) %>%
|
|
0
|
800 ungroup() %>%
|
|
9
|
801 select(-group);
|
|
0
|
802 # prep_genotype_table_tibble looks like this:
|
|
|
803 # affy_id coral_mlg_clonal_id user_specimen_id db_match
|
|
|
804 # a550962...CEL HG0120 1090 match
|
|
|
805 # genetic_coral_species_call coral_mlg_rep_sample_id
|
|
|
806 # A.palmata 1104
|
|
|
807 prep_genotype_table_tibble <- stag_db_report %>%
|
|
|
808 select("affy_id", "coral_mlg_clonal_id", "user_specimen_id", "db_match", "genetic_coral_species_call") %>%
|
|
|
809 left_join(representative_mlg_tibble %>%
|
|
9
|
810 select("coral_mlg_rep_sample_id", "coral_mlg_clonal_id"),
|
|
|
811 by='coral_mlg_clonal_id');
|
|
0
|
812 # genotype_table_tibble looks like this:
|
|
|
813 # affy_id coral_mlg_clonal_id user_specimen_id db_match
|
|
|
814 # a550962-436.CEL HG0120 1090 match
|
|
|
815 # genetic_coral_species_call coral_mlg_rep_sample_id bcoral_genet_id
|
|
|
816 # A.palmata 1104 <NA>
|
|
|
817 genotype_table_tibble <- prep_genotype_table_tibble %>%
|
|
|
818 left_join(affy_metadata_data_frame %>%
|
|
|
819 select("user_specimen_id", "bcoral_genet_id"),
|
|
9
|
820 by='user_specimen_id');
|
|
0
|
821 write_data_frame(output_data_dir, "genotype.tabular", genotype_table_tibble);
|
|
|
822
|
|
|
823 # Output the file needed for populating the person table.
|
|
|
824 person_table_data_frame <- data.frame(matrix(ncol=4, nrow=num_rows));
|
|
|
825 colnames(person_table_data_frame) <- c("last_name", "first_name", "organization", "email");
|
|
|
826 # person_table_data_frame looks like this:
|
|
|
827 # last_name first_name organization email
|
|
|
828 # Kitchen Sheila Penn State s89@psu.edu
|
|
|
829 for (i in 1:num_rows) {
|
|
|
830 person_table_data_frame$last_name[i] <- sample_prep_data_frame$collector_last_name[i];
|
|
|
831 person_table_data_frame$first_name[i] <- sample_prep_data_frame$collector_first_name[i];
|
|
|
832 person_table_data_frame$organization[i] <- sample_prep_data_frame$organization[i];
|
|
|
833 person_table_data_frame$email[i] <- sample_prep_data_frame$email[i];
|
|
|
834 }
|
|
|
835 write_data_frame(output_data_dir, "person.tabular", person_table_data_frame);
|
|
|
836
|
|
|
837 # Output the file needed for populating the phenotype table.
|
|
|
838 phenotype_table_data_frame <- data.frame(matrix(ncol=7, nrow=num_rows));
|
|
|
839 colnames(phenotype_table_data_frame) <- c("disease_resist", "bleach_resist", "mortality", "tle",
|
|
|
840 "spawning", "sperm_motility", "healing_time");
|
|
|
841 # phenotype_table_data_frame looks like this:
|
|
|
842 # disease_resist bleach_resist mortality tle spawning sperm_motility healing_time
|
|
|
843 # NA NA NA NA False NA NA
|
|
|
844 for (i in 1:num_rows) {
|
|
|
845 phenotype_table_data_frame$disease_resist[i] <- sample_prep_data_frame$disease_resist[i];
|
|
|
846 phenotype_table_data_frame$bleach_resist[i] <- sample_prep_data_frame$bleach_resist[i];
|
|
|
847 phenotype_table_data_frame$mortality[i] <- sample_prep_data_frame$mortality[i];
|
|
|
848 phenotype_table_data_frame$tle[i] <- sample_prep_data_frame$tle[i];
|
|
|
849 phenotype_table_data_frame$spawning[i] <- sample_prep_data_frame$spawning[i];
|
|
|
850 phenotype_table_data_frame$sperm_motility[i] <- sample_prep_data_frame$sperm_motility[i];
|
|
|
851 phenotype_table_data_frame$healing_time[i] <- sample_prep_data_frame$healing_time[i];
|
|
|
852 }
|
|
|
853 write_data_frame(output_data_dir, "phenotype.tabular", phenotype_table_data_frame);
|
|
|
854
|
|
|
855 # Output the file needed for populating the reef table.
|
|
|
856 reef_table_data_frame <- data.frame(matrix(ncol=5, nrow=num_rows));
|
|
|
857 colnames(reef_table_data_frame) <- c("name", "region", "latitude", "longitude", "geographic_origin");
|
|
|
858 # The geographic_origin value is used for deciding into which table
|
|
|
859 # to insert the latitude and longitude values. If the geographic_origin
|
|
|
860 # is "reef", the values will be inserted into the reef table, and if it is
|
|
|
861 # "colony", the values will be inserted into the colony table. We insert
|
|
|
862 # these values in both data frames so that the downstream tool that parses
|
|
|
863 # them can determine the appropriate table.
|
|
|
864 # reef_table_data_frame looks like this:
|
|
|
865 # name region latitude longitude geographic_origin
|
|
|
866 # JohnsonsReef Bahamas 18.361733 -64.7743 Reef
|
|
|
867 for (i in 1:num_rows) {
|
|
|
868 reef_table_data_frame$name[i] <- sample_prep_data_frame$reef[i];
|
|
|
869 reef_table_data_frame$region[i] <- sample_prep_data_frame$region[i];
|
|
|
870 reef_table_data_frame$latitude[i] <- sample_prep_data_frame$latitude[i];
|
|
|
871 reef_table_data_frame$longitude[i] <- sample_prep_data_frame$longitude[i];
|
|
|
872 reef_table_data_frame$geographic_origin[i] <- sample_prep_data_frame$geographic_origin[i];
|
|
|
873 }
|
|
|
874 write_data_frame(output_data_dir, "reef.tabular", reef_table_data_frame);
|
|
|
875
|
|
|
876 # Output the file needed for populating the sample table.
|
|
|
877 sample_table_data_frame <- data.frame(matrix(ncol=20, nrow=num_rows));
|
|
|
878 colnames(sample_table_data_frame) <- c("affy_id", "colony_location", "collection_date", "user_specimen_id",
|
|
|
879 "registry_id", "depth", "dna_extraction_method", "dna_concentration",
|
|
|
880 "public", "public_after_date", "percent_missing_data_coral",
|
|
1
|
881 "percent_missing_data_sym", "percent_acerv_coral",
|
|
|
882 "percent_reference_sym", "percent_apalm_coral",
|
|
0
|
883 "percent_alternative_sym", "percent_heterozygous_coral",
|
|
|
884 "percent_heterozygous_sym", "field_call", "bcoral_genet_id");
|
|
|
885 for (i in 1:num_rows) {
|
|
|
886 sample_table_data_frame$affy_id[i] <- sample_prep_data_frame$affy_id[i];
|
|
|
887 sample_table_data_frame$colony_location[i] <- sample_prep_data_frame$colony_location[i];
|
|
|
888 sample_table_data_frame$collection_date[i] <- sample_prep_data_frame$collection_date[i];
|
|
|
889 sample_table_data_frame$user_specimen_id[i] <- sample_prep_data_frame$user_specimen_id[i];
|
|
|
890 sample_table_data_frame$registry_id[i] <- sample_prep_data_frame$registry_id[i];
|
|
|
891 sample_table_data_frame$depth[i] <- sample_prep_data_frame$depth[i];
|
|
|
892 sample_table_data_frame$dna_extraction_method[i] <- sample_prep_data_frame$dna_extraction_method[i];
|
|
|
893 sample_table_data_frame$dna_concentration[i] <- sample_prep_data_frame$dna_concentration[i];
|
|
|
894 sample_table_data_frame$public[i] <- sample_prep_data_frame$public[i];
|
|
|
895 sample_table_data_frame$public_after_date[i] <- sample_prep_data_frame$public_after_date[i];
|
|
|
896 sample_table_data_frame$percent_missing_data_coral[i] <- sample_prep_data_frame$percent_missing_data_coral[i];
|
|
|
897 sample_table_data_frame$percent_missing_data_sym[i] <- DEFAULT_MISSING_NUMERIC_VALUE;
|
|
1
|
898 sample_table_data_frame$percent_acerv_coral[i] <- sample_prep_data_frame$percent_acerv_coral[i];
|
|
0
|
899 sample_table_data_frame$percent_reference_sym[i] <- DEFAULT_MISSING_NUMERIC_VALUE;
|
|
1
|
900 sample_table_data_frame$percent_apalm_coral[i] <- sample_prep_data_frame$percent_apalm_coral[i];
|
|
0
|
901 sample_table_data_frame$percent_alternative_sym[i] <- DEFAULT_MISSING_NUMERIC_VALUE;
|
|
|
902 sample_table_data_frame$percent_heterozygous_coral[i] <- sample_prep_data_frame$percent_heterozygous_coral[i];
|
|
|
903 sample_table_data_frame$percent_heterozygous_sym[i] <- DEFAULT_MISSING_NUMERIC_VALUE;
|
|
|
904 sample_table_data_frame$field_call[i] <- sample_prep_data_frame$field_call[i];
|
|
9
|
905 sample_table_data_frame$bcoral_genet_id[i] <- sample_prep_data_frame$bcoral_genet_id[i];
|
|
0
|
906 }
|
|
|
907 write_data_frame(output_data_dir, "sample.tabular", sample_table_data_frame);
|
|
|
908
|
|
|
909 # Output the file needed for populating the taxonomy table.
|
|
|
910 # taxonomy_table_data_frame looks like this:
|
|
|
911 # genetic_coral_species_call affy_id genus_name species_name
|
|
|
912 # A.palmata a550962-4368120-060520-500_A05.CEL Acropora palmata
|
|
|
913 taxonomy_table_data_frame <- stag_db_report %>%
|
|
|
914 select(genetic_coral_species_call, affy_id) %>%
|
|
2
|
915 mutate(genus_name = ifelse(grepl("^A.*", genetic_coral_species_call), "Acropora",ifelse(!is.na(genetic_coral_species_call), "other", NA))) %>%
|
|
0
|
916 mutate(species_name = ifelse(genetic_coral_species_call == "A.palmata", "palmata", "other")) %>%
|
|
|
917 mutate(species_name = ifelse(genetic_coral_species_call == "A.cervicornis", "cervicornis", species_name)) %>%
|
|
|
918 mutate(species_name = ifelse(genetic_coral_species_call == "A.prolifera", "prolifera", species_name));
|
|
|
919 colnames(taxonomy_table_data_frame) <- c("genetic_coral_species_call", "affy_id", "genus_name", "species_name");
|
|
|
920 write_data_frame(output_data_dir, "taxonomy.tabular", taxonomy_table_data_frame);
|
|
|
921 time_elapsed(start_time);
|
