annotate coral_multilocus_genotype.R @ 0:adaf89535d2e draft

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