|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("adegenet"))
|
|
|
4 suppressPackageStartupMessages(library("ape"))
|
|
7
|
5 suppressPackageStartupMessages(library("data.table"))
|
|
12
|
6 suppressPackageStartupMessages(library("dbplyr"))
|
|
7
|
7 suppressPackageStartupMessages(library("dplyr"))
|
|
0
|
8 suppressPackageStartupMessages(library("ggplot2"))
|
|
|
9 suppressPackageStartupMessages(library("knitr"))
|
|
4
|
10 suppressPackageStartupMessages(library("optparse"))
|
|
|
11 suppressPackageStartupMessages(library("poppr"))
|
|
|
12 suppressPackageStartupMessages(library("RColorBrewer"))
|
|
17
|
13 suppressPackageStartupMessages(library("rnaturalearth"))
|
|
|
14 suppressPackageStartupMessages(library("rnaturalearthdata"))
|
|
7
|
15 suppressPackageStartupMessages(library("RPostgres"))
|
|
17
|
16 suppressPackageStartupMessages(library("sf"))
|
|
|
17 suppressPackageStartupMessages(library(SNPRelate))
|
|
12
|
18 suppressPackageStartupMessages(library("tidyr"))
|
|
4
|
19 suppressPackageStartupMessages(library("vcfR"))
|
|
|
20 suppressPackageStartupMessages(library("vegan"))
|
|
17
|
21 suppressPackageStartupMessages(library("yarrr"))
|
|
|
22 theme_set(theme_bw())
|
|
0
|
23
|
|
|
24 option_list <- list(
|
|
7
|
25 make_option(c("--database_connection_string"), action="store", dest="database_connection_string", help="Corals (stag) database connection string"),
|
|
|
26 make_option(c("--input_affy_metadata"), action="store", dest="input_affy_metadata", help="Affymetrix 96 well plate input file"),
|
|
4
|
27 make_option(c("--input_pop_info"), action="store", dest="input_pop_info", help="Population information input file"),
|
|
7
|
28 make_option(c("--input_vcf"), action="store", dest="input_vcf", help="VCF input file"),
|
|
17
|
29 make_option(c("--output_stag_db_report"), action="store", dest="output_stag_db_report", help="stag db report output file"),
|
|
|
30 make_option(c("--nj_tree"), action="store", dest="nj_tree", help="neighbor-joining tree output file")
|
|
0
|
31 )
|
|
|
32
|
|
|
33 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
|
|
|
34 args <- parse_args(parser, positional_arguments=TRUE);
|
|
|
35 opt <- args$options;
|
|
|
36
|
|
|
37 get_file_path = function(file_name) {
|
|
|
38 file_path = paste("output_plots_dir", file_name, sep="/");
|
|
|
39 return(file_path);
|
|
|
40 }
|
|
|
41
|
|
8
|
42 get_database_connection <- function(db_conn_string) {
|
|
|
43 # Instantiate database connection.
|
|
|
44 # The connection string has this format:
|
|
|
45 # postgresql://user:password@host/dbname
|
|
|
46 conn_items <- strsplit(db_conn_string, "://")[[1]];
|
|
|
47 string_needed <- conn_items[2];
|
|
|
48 items_needed <- strsplit(string_needed, "@")[[1]];
|
|
|
49 user_pass_string <- items_needed[1];
|
|
|
50 host_dbname_string <- items_needed[2];
|
|
|
51 user_pass_items <- strsplit(user_pass_string, ":")[[1]];
|
|
|
52 host_dbname_items <- strsplit(host_dbname_string, "/")[[1]];
|
|
|
53 user <- user_pass_items[1];
|
|
|
54 pass <- user_pass_items[2];
|
|
|
55 host <- host_dbname_items[1];
|
|
|
56 dbname <- host_dbname_items[2];
|
|
|
57 # FIXME: is there a way to not hard-code the port?
|
|
16
|
58 conn <- DBI::dbConnect(RPostgres::Postgres(), host=host, port="5432", dbname=dbname, user=user, password=pass);
|
|
8
|
59 return (conn);
|
|
|
60 }
|
|
|
61
|
|
3
|
62 # Read in VCF input file.
|
|
2
|
63 vcf <- read.vcfR(opt$input_vcf);
|
|
0
|
64
|
|
7
|
65 # Convert VCF file into a genind for the Poppr package.
|
|
|
66 # TODO: probably should not hard-code 2 cores.
|
|
17
|
67 # changed to genind format for extracting alleles later
|
|
|
68 # trade-off is it is a bit slower to import data
|
|
|
69 # gl <- vcfR2genlight(vcf, n.cores=2)
|
|
|
70 # gind <- new("genind", (as.matrix(gl)))
|
|
|
71 gind <- vcfR2genind(vcf);
|
|
7
|
72
|
|
0
|
73 # Add population information to the genind object.
|
|
10
|
74 poptab <- read.table(opt$input_pop_info, check.names=FALSE, header=F, na.strings=c("", "NA"), stringsAsFactors=FALSE, sep="\t");
|
|
9
|
75 colnames(poptab) <- c("row_id", "affy_id", "user_specimen_id", "region");
|
|
8
|
76 gind@pop <- as.factor(poptab$region);
|
|
17
|
77 strata(gind)<-data.frame(pop(gind));
|
|
7
|
78
|
|
|
79 # Convert genind object to a genclone object.
|
|
8
|
80 obj2 <- as.genclone(gind);
|
|
7
|
81
|
|
|
82 # Calculate the bitwise distance between individuals.
|
|
8
|
83 xdis <- bitwise.dist(obj2);
|
|
0
|
84
|
|
17
|
85 # Multilocus genotypes (threshold of 3.2%).
|
|
|
86 # threshold doubled because of how the data is formatted in genind compared to genlight
|
|
|
87 mlg.filter(obj2, distance=xdis) <- 0.032;
|
|
8
|
88 m <- mlg.table(obj2, background=TRUE, color=TRUE);
|
|
0
|
89
|
|
|
90 # Create table of MLGs.
|
|
8
|
91 id <- mlg.id(obj2);
|
|
17
|
92 #dt <- data.table(id, keep.rownames=TRUE);
|
|
|
93 #setnames(dt, c("id"), c("affy_id"));
|
|
7
|
94
|
|
16
|
95 # Read user's Affymetrix 96 well plate tabular file.
|
|
17
|
96 pinfo <- read.table(opt$input_affy_metadata, header=FALSE, stringsAsFactors=FALSE, sep="\t", na.strings = c("", "NA"));
|
|
16
|
97 colnames(pinfo) <- c("user_specimen_id", "field_call", "bcoral_genet_id", "bsym_genet_id", "reef",
|
|
|
98 "region", "latitude", "longitude", "geographic_origin", "sample_location",
|
|
17
|
99 "latitude_outplant", "longitude_outplant", "depth", "disease_resist",
|
|
16
|
100 "bleach_resist", "mortality","tle", "spawning", "collector_last_name",
|
|
17
|
101 "collector_first_name", "organization", "collection_date", "email", "seq_facility",
|
|
16
|
102 "array_version", "public", "public_after_date", "sperm_motility", "healing_time",
|
|
17
|
103 "dna_extraction_method", "dna_concentration", "registry_id");
|
|
9
|
104 pinfo$user_specimen_id <- as.character(pinfo$user_specimen_id);
|
|
|
105 pinfo2 <- as.character(pinfo$user_specimen_id);
|
|
17
|
106 pi <- data.table(pinfo2, pinfo$field_call);
|
|
9
|
107 setnames(pi, c("pinfo2"), c("user_specimen_id"));
|
|
17
|
108 setnames(pi, c("V2"), c("field_call"));
|
|
7
|
109
|
|
8
|
110 # Connect to database.
|
|
|
111 conn <- get_database_connection(opt$database_connection_string);
|
|
7
|
112
|
|
|
113 # Import the sample table.
|
|
16
|
114 sample_table <- tbl(conn, "sample");
|
|
|
115
|
|
|
116 # Import the genotype table.
|
|
|
117 genotype_table <- tbl(conn, "genotype");
|
|
7
|
118
|
|
16
|
119 # Select columns from the sample table and the
|
|
|
120 # genotype table joined by genotype_id.
|
|
|
121 sample_table_columns <- sample_table %>% select(user_specimen_id, affy_id, genotype_id);
|
|
|
122 smlg <- sample_table_columns %>%
|
|
|
123 left_join(genotype_table %>%
|
|
|
124 select("id", "coral_mlg_clonal_id", "symbio_mlg_clonal_id"),
|
|
|
125 by=c("genotype_id" = "id"));
|
|
7
|
126
|
|
|
127 # Convert to dataframe.
|
|
|
128 sm <- data.frame(smlg);
|
|
16
|
129 # Name the columns.
|
|
|
130 colnames(sm) <- c("user_specimen_id", "affy_id", "genotype_id", "coral_mlg_clonal_id", "symbio_mlg_clonal_id");
|
|
7
|
131
|
|
12
|
132 # Missing GT in samples submitted.
|
|
|
133 gt <- extract.gt(vcf, element="GT", as.numeric=FALSE);
|
|
|
134 myMiss <- apply(gt, MARGIN=2, function(x){ sum(is.na(x))});
|
|
|
135 myMiss <- (myMiss / nrow(vcf)) * 100;
|
|
|
136 miss <- data.frame(myMiss);
|
|
|
137
|
|
7
|
138 # Convert missing data into data table.
|
|
8
|
139 mi <-setDT(miss, keep.rownames=TRUE)[];
|
|
9
|
140 setnames(mi, c("rn"), c("affy_id"));
|
|
7
|
141 setnames(mi, c("myMiss"), c("percent_missing_data_coral"));
|
|
|
142 # Round missing data to two digits.
|
|
9
|
143 mi$percent_missing_data_coral <- round(mi$percent_missing_data_coral, digits=2);
|
|
7
|
144
|
|
17
|
145 #heterozygous alleles
|
|
12
|
146 hets <- apply(gt, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("0/1", x))))} );
|
|
|
147 hets <- (hets / nrow(vcf)) * 100;
|
|
|
148 ht <- data.frame(hets);
|
|
|
149
|
|
7
|
150 # Convert heterozygosity data into data table.
|
|
|
151 ht <-setDT(ht, keep.rownames=TRUE)[];
|
|
9
|
152 setnames(ht, c("rn"), c("affy_id"));
|
|
7
|
153 setnames(ht, c("hets"), c("percent_mixed_coral"));
|
|
|
154 # Round missing data to two digits.
|
|
|
155 ht$percent_mixed<-round(ht$percent_mixed, digits=2);
|
|
|
156
|
|
17
|
157 #reference alleles
|
|
12
|
158 refA <- apply(gt, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("0/0", x))))} );
|
|
|
159 refA <- (refA / nrow(vcf)) * 100;
|
|
|
160 rA <- data.frame(refA);
|
|
|
161
|
|
7
|
162 # Convert refA data into data.table.
|
|
|
163 rA <-setDT(rA, keep.rownames=TRUE)[];
|
|
9
|
164 setnames(rA, c("rn"), c("affy_id"));
|
|
7
|
165 setnames(rA, c("refA"), c("percent_reference_coral"));
|
|
|
166 # round missing data to two digits.
|
|
|
167 rA$percent_reference<-round(rA$percent_reference, digits=2);
|
|
|
168
|
|
17
|
169 #alternative alleles
|
|
12
|
170 altB <- apply(gt, MARGIN=2, function(x) {sum(lengths(regmatches(x, gregexpr("1/1", x))))} );
|
|
|
171 altB <- (altB / nrow(vcf)) * 100;
|
|
|
172 aB <- data.frame(altB);
|
|
|
173
|
|
7
|
174 # Convert altB data into data table.
|
|
|
175 aB <-setDT(aB, keep.rownames=TRUE)[];
|
|
9
|
176 setnames(aB, c("rn"), c("affy_id"));
|
|
7
|
177 setnames(aB, c("altB"), c("percent_alternative_coral"));
|
|
|
178 # Round missing data to two digits.
|
|
|
179 aB$percent_alternative<-round(aB$percent_alternative, digits=2);
|
|
|
180
|
|
|
181 #convert mlg id to data.table format
|
|
|
182 dt <- data.table(id, keep.rownames=TRUE);
|
|
9
|
183 setnames(dt, c("id"), c("affy_id"));
|
|
7
|
184
|
|
|
185 # Transform.
|
|
|
186 df3 <- dt %>%
|
|
|
187 group_by(row_number()) %>%
|
|
16
|
188 dplyr::rename(group="row_number()") %>%
|
|
12
|
189 unnest (affy_id) %>%
|
|
7
|
190 # Join with mlg table.
|
|
|
191 left_join(sm %>%
|
|
9
|
192 select("affy_id","coral_mlg_clonal_id"),
|
|
16
|
193 by="affy_id");
|
|
7
|
194
|
|
|
195 # If found in database, group members on previous mlg id.
|
|
|
196 uniques <- unique(df3[c("group", "coral_mlg_clonal_id")]);
|
|
|
197 uniques <- uniques[!is.na(uniques$coral_mlg_clonal_id),];
|
|
|
198 na.mlg <- which(is.na(df3$coral_mlg_clonal_id));
|
|
|
199 na.group <- df3$group[na.mlg];
|
|
|
200 df3$coral_mlg_clonal_id[na.mlg] <- uniques$coral_mlg_clonal_id[match(na.group, uniques$group)];
|
|
|
201
|
|
|
202 # Determine if the sample mlg matched previous genotyped sample.
|
|
|
203 df4<- df3 %>%
|
|
|
204 group_by(group) %>%
|
|
12
|
205 mutate(DB_match = ifelse(is.na(coral_mlg_clonal_id),"no_match", "match"));
|
|
7
|
206
|
|
|
207 # Create new mlg id for samples that did not match those in the database.
|
|
|
208 none <- unique(df4[c("group", "coral_mlg_clonal_id")]);
|
|
|
209 none <- none[is.na(none$coral_mlg_clonal_id),];
|
|
|
210 na.mlg2 <- which(is.na(df4$coral_mlg_clonal_id));
|
|
|
211 n.g <- df4$group[na.mlg2];
|
|
|
212 ct <- length(unique(n.g));
|
|
|
213
|
|
|
214 # List of new group ids, the sequence starts at the number of
|
|
|
215 # ids present in df4$coral_mlg_clonal_ids plus 1. Not sure if
|
|
|
216 # the df4 file contains all ids. If it doesn't then look below
|
|
17
|
217 # to change the seq() function.
|
|
7
|
218 n.g_ids <- sprintf("HG%04d", seq((sum(!is.na(unique(df4["coral_mlg_clonal_id"]))) + 1), by=1, length=ct));
|
|
16
|
219 # Pair group with new ids.
|
|
7
|
220 rat <- cbind(unique(n.g), n.g_ids);
|
|
16
|
221 # Assign the new id iteratively for all that have NA.
|
|
7
|
222 for (i in 1:length(na.mlg2)) {
|
|
|
223 df4$coral_mlg_clonal_id[na.mlg2[i]] <- n.g_ids[match(df4$group[na.mlg2[i]], unique(n.g))];
|
|
|
224 }
|
|
|
225
|
|
16
|
226 # Subset poptab for all samples.
|
|
9
|
227 subpop <- poptab[c(2, 3)];
|
|
|
228
|
|
7
|
229 # Merge data frames for final table.
|
|
|
230 report_user <- pi %>%
|
|
9
|
231 left_join(subpop %>%
|
|
|
232 select("affy_id", "user_specimen_id"),
|
|
16
|
233 by="user_specimen_id") %>%
|
|
9
|
234 left_join(df4 %>%
|
|
|
235 select("affy_id", "coral_mlg_clonal_id", "DB_match"),
|
|
16
|
236 by="affy_id") %>%
|
|
7
|
237 left_join(mi %>%
|
|
9
|
238 select("affy_id", "percent_missing_data_coral"),
|
|
16
|
239 by="affy_id") %>%
|
|
7
|
240 left_join(ht %>%
|
|
9
|
241 select("affy_id", "percent_mixed_coral"),
|
|
16
|
242 by="affy_id") %>%
|
|
7
|
243 left_join(rA %>%
|
|
9
|
244 select("affy_id", "percent_reference_coral"),
|
|
16
|
245 by="affy_id") %>%
|
|
7
|
246 left_join(aB %>%
|
|
9
|
247 select("affy_id", "percent_alternative_coral"),
|
|
16
|
248 by="affy_id") %>%
|
|
7
|
249 mutate(DB_match = ifelse(is.na(DB_match), "failed", DB_match))%>%
|
|
12
|
250 mutate(coral_mlg_clonal_id = ifelse(is.na(coral_mlg_clonal_id), "failed", coral_mlg_clonal_id)) %>%
|
|
17
|
251 mutate(genetic_coral_species_call=ifelse(percent_alternative_coral >= 40 & percent_alternative_coral<= 44.5,"A.palmata","other")) %>%
|
|
|
252 mutate(genetic_coral_species_call=ifelse(percent_alternative_coral >= 45.5 & percent_alternative_coral<= 50,"A.cervicornis",genetic_coral_species_call)) %>%
|
|
|
253 mutate(genetic_coral_species_call=ifelse(percent_heterozygous_coral > 40,"A.prolifera",genetic_coral_species_call)) %>%
|
|
7
|
254 ungroup() %>%
|
|
|
255 select(-group);
|
|
|
256
|
|
12
|
257 write.csv(report_user, file=opt$output_stag_db_report, quote=FALSE);
|
|
0
|
258
|
|
17
|
259 # Database tables
|
|
|
260 ## Sample.table
|
|
|
261 sample_db <- pinfo %>%
|
|
|
262 left_join(
|
|
|
263 report_user %>%
|
|
|
264 select("user_specimen_id","affy_id",
|
|
|
265 "percent_missing_data_coral","percent_heterozygous_coral","percent_reference_coral",
|
|
|
266 "percent_alternative_coral"),
|
|
|
267 by='user_specimen_id');
|
|
|
268
|
|
|
269 ###representative clone for genotype.table
|
|
|
270 cc<-clonecorrect(obj2, strata= ~pop.gind.);
|
|
|
271 id_rep<-mlg.id(cc);
|
|
|
272 dt_cc<-data.table(id_rep,keep.rownames = TRUE);
|
|
|
273 setnames(dt_cc, c("id_rep"), c("affy_id"));
|
|
|
274
|
|
|
275 ###transform mlg data.table
|
|
|
276 df_cc <- dt_cc %>%
|
|
|
277 group_by(row_number()) %>%
|
|
|
278 rename(group='row_number()') %>%
|
|
|
279 unnest(affy_id) %>%
|
|
|
280 left_join(report_user %>%
|
|
|
281 select("coral_mlg_clonal_id","user_specimen_id","affy_id"),
|
|
|
282 by='affy_id') %>%
|
|
|
283 mutate(coral_mlg_rep_sample_id=ifelse(is.na(coral_mlg_clonal_id),"",affy_id)) %>%
|
|
|
284 ungroup() %>%
|
|
|
285 select(-group);
|
|
7
|
286
|
|
17
|
287 ##geno.table
|
|
|
288 geno_db <- df4 %>%
|
|
|
289 left_join(df_cc %>%
|
|
|
290 select("affy_id","coral_mlg_rep_sample_id","user_specimen_id"),
|
|
|
291 by='affy_id') %>%
|
|
|
292 ungroup() %>%
|
|
|
293 select(-group);
|
|
|
294
|
|
|
295 ##taxonomy.table
|
|
|
296
|
|
|
297 tax_db <- report_user %>%
|
|
|
298 select(genetic_coral_species_call, affy_id) %>%
|
|
|
299 mutate(genus_name =ifelse(genetic_coral_species_call==
|
|
|
300 genetic_coral_species_call[grep("^A.*",genetic_coral_species_call)],"Acropora","other")) %>%
|
|
|
301 mutate(species_name=ifelse(genetic_coral_species_call=="A.palmata","palmata","other"))%>%
|
|
|
302 mutate(species_name=ifelse(genetic_coral_species_call =="A.cervicornis","cervicornis",species_name))%>%
|
|
|
303 mutate(species_name=ifelse(genetic_coral_species_call=="A.prolifera","prolifera", species_name));
|
|
|
304
|
|
|
305
|
|
|
306 # Table of alleles for the new samples
|
|
|
307 ## subset to new plate data
|
|
|
308 ### create vector indicating number of individuals desired
|
|
|
309 ### made from affy_id collumn from report_user data table
|
|
|
310 i<-ifelse(is.na(report_user[1]),"",report_user[[1]]);
|
|
|
311 i<-i[!apply(i == "", 1, all),];
|
|
|
312 sub96<-obj2[i, mlg.reset = FALSE, drop = FALSE];
|
|
|
313
|
|
|
314 # convert to data frame
|
|
|
315 at_96<-genind2df(sub96, sep="");
|
|
|
316 at_96<- at_96 %>%
|
|
|
317 select(-pop);
|
|
|
318
|
|
|
319 # allele string for Allele.table in database
|
|
|
320 uat_96<-unite(at_96, alleles, 1:19696, sep = " ", remove = TRUE);
|
|
|
321 uat_96<-setDT(uat_96, keep.rownames = TRUE)[];
|
|
|
322 setnames(uat_96, c("rn"), c("user_specimen_id"));
|
|
|
323 # write.csv(uat_96,file=paste("Seed_genotype_alleles.csv",sep = ""),quote=FALSE,row.names=FALSE);
|
|
0
|
324
|
|
4
|
325 # Create a phylogeny of samples based on distance matrices.
|
|
17
|
326 cols <- piratepal("basel");
|
|
4
|
327 set.seed(999);
|
|
|
328 # Start PDF device driver.
|
|
|
329 dev.new(width=10, height=7);
|
|
|
330 file_path = get_file_path("nj_phylogeny.pdf");
|
|
|
331 pdf(file=file_path, width=10, height=7);
|
|
|
332 # Organize branches by clade.
|
|
9
|
333 theTree <- sub96 %>%
|
|
17
|
334 aboot(dist=provesti.dist, sample=100, tree="nj", cutoff=50, quiet=TRUE) %>%
|
|
7
|
335 ladderize();
|
|
9
|
336 theTree$tip.label <- report_user$user_specimen_id[match(theTree$tip.label, report_user$affy_id)];
|
|
|
337 plot.phylo(theTree, tip.color=cols[sub96$pop], label.offset=0.0125, cex=0.3, font=2, lwd=4, align.tip.label=F, no.margin=T);
|
|
4
|
338 # Add a scale bar showing 5% difference..
|
|
9
|
339 add.scale.bar(0, 0.95, length=0.05, cex=0.65, lwd=3);
|
|
|
340 nodelabels(theTree$node.label, cex=.5, adj=c(1.5, -0.1), frame="n", font=3, xpd=TRUE);
|
|
17
|
341 legend("topright", legend=c(levels(sub96$pop)), text.col=cols, xpd=T, cex=0.8);
|
|
|
342 dev.off()
|
|
|
343
|
|
|
344 write.tree(theTree, file =opt$nj_tree, quote=FALSE);
|
|
|
345
|
|
|
346 # identity-by-state analysis
|
|
|
347 #if (!requireNamespace("BiocManager", quietly = TRUE))
|
|
|
348 # install.packages("BiocManager")
|
|
|
349 #BiocManager::install("SNPRelate", version = "3.8")
|
|
|
350
|
|
|
351 #subset VCF to the user samples
|
|
|
352 l<-length(i);
|
|
|
353 n<-ncol(vcf@gt);
|
|
|
354 s<-n-l;
|
|
|
355 svcf<-vcf[,s:n];
|
|
|
356 write.vcf(svcf, "subset.vcf.gz");
|
|
|
357
|
|
|
358 vcf.fn <- "subset.vcf.gz";
|
|
|
359 snpgdsVCF2GDS(vcf.fn, "test3.gds", method="biallelic.only");
|
|
|
360
|
|
|
361 genofile <- snpgdsOpen(filename="test3.gds", readonly=FALSE);
|
|
|
362 hd<-read.gdsn(index.gdsn(genofile, "sample.id"));
|
|
|
363 hd<-data.frame(hd);
|
|
|
364 hd<-setDT(hd, keep.rownames = FALSE)[];
|
|
|
365 setnames(hd, c("hd"), c("user_specimen_id"));
|
|
|
366
|
|
|
367 subpop2<- poptab[c(2,4)];
|
|
|
368 poptab_sub <- hd %>%
|
|
|
369 left_join(
|
|
|
370 subpop2 %>%
|
|
|
371 select("affy_id","region"),
|
|
|
372 by='affy_id')%>%
|
|
|
373 drop_na();
|
|
|
374
|
|
|
375 samp.annot <- data.frame(pop.group = c(poptab_sub$region));
|
|
|
376 add.gdsn(genofile, "sample.annot", samp.annot);
|
|
|
377
|
|
|
378 pop_code <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"));
|
|
|
379 pop.group <- as.factor(read.gdsn(index.gdsn(genofile, "sample.annot/pop.group")));
|
|
|
380
|
|
|
381 # Identity-By-State Analysis - distance matrix calculation
|
|
|
382 ibs <- snpgdsIBS(genofile, num.thread=2, autosome.only=FALSE);
|
|
|
383
|
|
|
384 # cluster analysis on the genome-wide IBS pairwise distance matrix
|
|
|
385 set.seed(100);
|
|
|
386 par(cex=0.6, cex.lab=1, cex.axis=1.5,cex.main=2);
|
|
|
387 ibs.hc <- snpgdsHCluster(snpgdsIBS(genofile, autosome.only=FALSE));
|
|
|
388
|
|
|
389 # default clustering.
|
|
|
390 dev.new(width=10, height=7);
|
|
|
391 file_path = get_file_path("IBS_default.pdf");
|
|
|
392 pdf (file=file_path, width=10, height=7);
|
|
|
393 rv <- snpgdsCutTree(ibs.hc, col.list=cols, pch.list=15);
|
|
|
394 snpgdsDrawTree(rv, main="Color by Cluster", leaflab="perpendicular",y.label=0.2);
|
|
|
395 legend("topleft", legend=levels(rv$samp.group), xpd=T, col=cols[1:nlevels(rv$samp.group)], pch=15, ncol=4, cex=1.2);
|
|
|
396 dev.off()
|
|
|
397
|
|
|
398 # color cluster by region.
|
|
|
399 dev.new(width=10, height=7);
|
|
|
400 file_path = get_file_path("IBS_Region.pdf");
|
|
|
401 pdf (file=file_path, width=10, height=7);
|
|
|
402 race <- as.factor(pop_code);
|
|
|
403 rv2 <- snpgdsCutTree(ibs.hc,samp.group=race,col.list=cols,pch.list=15);
|
|
|
404 snpgdsDrawTree(rv2, main="Color by Region", leaflab="perpendicular",y.label=0.2);
|
|
|
405 legend("topleft", legend=levels(race), xpd=T, col=cols[1:nlevels(race)], pch=15, ncol=4, cex=1.2);
|
|
|
406 dev.off()
|
|
|
407
|
|
|
408 #close GDS file
|
|
|
409 snpgdsClose(genofile);
|
|
|
410
|
|
|
411 # Sample MLG on a map.
|
|
|
412 world <- ne_countries(scale = "medium", returnclass = "sf");
|
|
|
413 class(world);
|
|
|
414
|
|
|
415 pinfo$mlg<-report_user$coral_mlg_clonal_id;
|
|
|
416 n <- nrow(pinfo);
|
|
|
417
|
|
|
418 mxlat<-max(pinfo$latitude,na.rm = TRUE);
|
|
|
419 mnlat<-min(pinfo$latitude,na.rm = TRUE);
|
|
|
420 mxlong<-max(pinfo$longitude,na.rm = TRUE);
|
|
|
421 mnlong<-min(pinfo$longitude,na.rm = TRUE);
|
|
|
422
|
|
|
423 p5<-ggplot(data = world) +
|
|
|
424 geom_sf() +
|
|
|
425 coord_sf(xlim = c(mnlong-3, mxlong+3), ylim = c(mnlat-3,mxlat+3), expand = FALSE);
|
|
|
426
|
|
|
427 colourCount = length(unique(pinfo$mlg));
|
|
|
428 getPalette = colorRampPalette(piratepal("basel"));
|
|
|
429 dev.new(width=10, height=7);
|
|
|
430 file_path = get_file_path("mlg_map.pdf");
|
|
|
431 pdf (file=file_path, width=10, height=7);
|
|
|
432 p6<-p5+ geom_point(data = pinfo,aes(x =longitude, y=latitude, group=mlg, color = mlg), alpha=.7, size=3)+
|
|
|
433 scale_color_manual(values=getPalette(colourCount))+
|
|
|
434 theme(legend.position="bottom")+
|
|
|
435 guides(color=guide_legend(nrow=8,byrow=F));
|
|
|
436 p6;
|
|
|
437 dev.off()
|
|
0
|
438
|
|
9
|
439 # Missing data barplot.
|
|
|
440 poptab$miss <- report_user$percent_missing_data_coral[match(miss$affy_id, report_user$affy_id)];
|
|
|
441 test2 <- which(!is.na(poptab$miss));
|
|
|
442 miss96 <- poptab$miss[test2];
|
|
|
443 name96 <- poptab$user_specimen_id[test2];
|
|
|
444 dev.new(width=10, height=7);
|
|
|
445 file_path = get_file_path("missing_data.pdf");
|
|
|
446 pdf (file=file_path, width=10, height=7);
|
|
|
447 par(mar = c(8, 4, 4, 2));
|
|
|
448 x <- barplot(miss96, las=2, col=cols, ylim=c(0, 3), cex.axis=0.8, space=0.8, ylab="Missingness (%)", xaxt="n");
|
|
|
449 text(cex=0.6, x=x-0.25, y=-.05, name96, xpd=TRUE, srt=60, adj=1);
|
|
|
450 dev.off()
|
|
|
451
|
|
15
|
452 # Generate a pie chart for each sample with a genotype.
|
|
|
453 # Store the numerical and user_specimen_id values from
|
|
|
454 # report_user for the charts (user_specimen_id names
|
|
|
455 # will be used to label each chart).
|
|
12
|
456 dt1 <- data.table(report_user);
|
|
|
457 dt1 <- report_user[c(-2, -3, -4)];
|
|
|
458 dt1 <- na.omit(dt1);
|
|
15
|
459 # Translate to N (i.e., number of samples with a
|
|
|
460 # genotype) columns and 5 rows.
|
|
12
|
461 tdt1 <- t(dt1);
|
|
|
462 # Make another data table and transpose it the same as dt1 to
|
|
15
|
463 # get numerics. These will feed into the creation of N vectors.
|
|
12
|
464 dt2 <- data.table(report_user);
|
|
|
465 dt2 <- report_user[c(-1, -2, -3, -4)];
|
|
15
|
466 # Translate to N columns and 5 rows.
|
|
12
|
467 tdt2 <- t(dt2);
|
|
|
468 tdt1_matrix <- as.matrix(tdt1[-1,]);
|
|
15
|
469 # The number of columns is the number of samples with genotypes.
|
|
|
470 nc <- ncol(tdt1_matrix);
|
|
12
|
471 mode(tdt1_matrix) <- "numeric";
|
|
|
472 spy <- rowMeans(tdt1_matrix);
|
|
|
473 dev.new(width=10, height=7);
|
|
|
474 file_path = get_file_path("percent_breakdown.pdf");
|
|
|
475 pdf(file=file_path, width=10, height=7);
|
|
|
476 # Average pie of all samples.
|
|
|
477 labels <- paste(c("missing data", "mixed", "reference", "alternative"), " (", round(spy, 1), "%)", sep="");
|
|
|
478 col <- c("GREY", "#006DDB", "#24FF24", "#920000");
|
|
|
479 main <- "Average breakdown of SNP assignments across all samples";
|
|
|
480 pie(spy, labels=labels, radius=0.60, col=col, main=main, cex.main=.75);
|
|
|
481 par(mfrow=c(3, 2));
|
|
14
|
482 col <- c("GREY", "#006DDB", "#24FF24", "#920000");
|
|
15
|
483 for (i in 1:nc) {
|
|
14
|
484 tmp_labels <- paste(labels, " (", round(tdt1_matrix[,i], 1), "%)", sep="");
|
|
12
|
485 main <- paste("Breakdown of SNP assignments for", tdt1[1, i]);
|
|
14
|
486 pie(tdt1_matrix[,i], labels=tmp_labels, radius=0.90, col=col, main=main, cex.main=.85, cex=0.75);
|
|
12
|
487 }
|
|
|
488 dev.off()
|
