############################################################# # EasyECD v2018.0 # # First Written by CN Stephan & R Rollo, 2016 # # Automate DNA data analysis for eye colour SNPS using qPCR # ############################################################# ## Instructions: ## qPCR output file must be of the type ".csv" ## The qPCR file must contain data in the first four columns with the following headers on the first row: C1 = Well Position; C2 = Sample Name; C3 = SNP Assay Name; C4 = Call. ## "Well Position" should be an alpha-numeric code, e.g., A1. ## "Sample Name" data can be a number or text. ## "SNP Assay Name" should be one of: HERC2, OCA2, SLC24A4, MATP, TYR, IRF4 ## "Call" should be one of: Negative Control (NC), Undetermined, Homozygous G/G, Heterozygous A/G, Homozygous A/A, Homozygous C/C, Heterozygous C/T, Homozygous T/T, Heterozygous G/T, Homozygous G/G, Heterozygous C/G. ## For exemplar formatting, see CRANIOFACIALidentification.com ## Load libraries required.packages <- c("tcltk") new.packages <- required.packages[!(required.packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) library(tcltk) ## Display Logo Window R.base.dir <- system.file() setwd(paste(R.base.dir,"/../../doc/html",sep="")) image <- as.numeric(file.access(paste(R.base.dir,"/../../doc/html/EasyECD_logo.gif",sep=""), mode=0)) if (image==-1) { tkmessageBox(title="EasyECD", message="You are missing the EasyECD logo. \n\nTo use EasyECD you must save the EasyECD logo to the following path: \n\nC:/Program Files/R/R-[version]/doc/html. \n\n See CRANIOFACIALidentification.com for more information.",icon="warning") q() } if (image==0) { tt <- tktoplevel() tkwm.title(tt, "EasyECD") tkwm.geometry(tt,"600x280+550+350") image1<-tclVar() tcl("image","create","photo",image1,file="EasyECD_logo.gif") imgAsLabel<-tklabel(tt,image=image1,bg="white") tkpack(imgAsLabel) tkconfigure(tt,cursor="watch") Sys.sleep (10) tkconfigure(tt,cursor="arrow") tkdestroy(tt) ## Terms of use tkmessageBox(title="EasyECD", message="By using this program you agree to CRANIOFACIALidentification.com's 'Terms of Use'. \n\nNO WARRANTIES OR GUARANTEES OF ANY KIND are provided with EasyECD. You use EasyECD completely, and entirely, at your own risk. \n\nSee www.CRANIOFACIALidentification.com for further details. \n ------------------------------------------------------------------------\n\nInstructions: \n- The qPCR file must be of the type '.cvs' with four columns of data:\n C1 - Well position; C2 - Sample; C3 - SNP Assay; and C4- Call.\n",icon="info",type="ok") ## Set working directory directory <- tclvalue(tkchooseDirectory(title="Select the WORKING DIRECTORY, i.e., where your qPCR .cvs file is located and where you want the new data files written to.")) setwd(directory) ## Load the raw FSTD data, pull out landmarks, and make numeric matrix tddirectory <- tclvalue(tkgetOpenFile(title="Open the .csv data file containing the qPCR results you wish to analyse")) data <- read.csv(tddirectory, header = TRUE, na.strings = "") #Remove controls nafind <- as.numeric(as.character(data[,2])) data[,2] <- nafind data[,5:9] <- 10 data <- na.omit(data) #Snip-ID codes #rs12913832 = HERC #rs1800407 = OCA #rs12896399 = SLC24A4 #rs16891982 = MATP #rs1393350 = TYR #rs12203592 = IRF4 # SNP to minor allele #HERC; minor allele = A for (i in 1:length(data[,3])) { if (data[i,3] == "HERC2" & data[i,4] == "Homozygous G/G") {data[i,5] <- "0"} else {if (data[i,3] == "HERC2" & data[i,4] == "Heterozygous A/G") {data[i,5] <- "1"} else {if (data[i,3] == "HERC2" & data[i,4] == "Homozygous A/A") {data[i,5] <- "2"} }} #OCA; minor allele = T if (data[i,3] == "OCA2" & data[i,4] == "Homozygous C/C") {data[i,5] <- "0"} else {if (data[i,3] == "OCA2" & data[i,4] == "Heterozygous C/T") {data[i,5] <- "1"} else {if (data[i,3] == "OCA2" & data[i,4] == "Homozygous T/T") {data[i,5] <- "2"} }} #SLC24A4; minor allele = G if (data[i,3] == "SLC24A4" & data[i,4] == "Homozygous T/T") {data[i,5] <- "0"} else {if (data[i,3] == "SLC24A4" & data[i,4] == "Heterozygous G/T") {data[i,5] <- "1"} else {if (data[i,3] == "SLC24A4" & data[i,4] == "Homozygous G/G") {data[i,5] <- "2"} }} #MATP; minor allele = C if (data[i,3] == "MATP" & data[i,4] == "Homozygous G/G") {data[i,5] <- "0"} else {if (data[i,3] == "MATP" & data[i,4] == "Heterozygous C/G") {data[i,5] <- "1"} else {if (data[i,3] == "MATP" & data[i,4] == "Homozygous C/C") {data[i,5] <- "2"} }} #TYR; minor allele = A if (data[i,3] == "TYR" & data[i,4] == "Homozygous G/G") {data[i,5] <- "0"} else {if (data[i,3] == "TYR" & data[i,4] == "Heterozygous A/G") {data[i,5] <- "1"} else {if (data[i,3] == "TYR" & data[i,4] == "Homozygous A/A") {data[i,5] <- "2"} }} #IRF4; minor allele = T if (data[i,3] == "IRF4" & data[i,4] == "Homozygous C/C") {data[i,5] <- "0"} else {if (data[i,3] == "IRF4" & data[i,4] == "Heterozygous C/T") {data[i,5] <- "1"} else {if (data[i,3] == "IRF4" & data[i,4] == "Homozygous T/T") {data[i,5] <- "2"} }} } ## order by samples data <- data[order(data[,2], data[,3]),] del <- which(data[,5]==10) samp <- unique(data[del,2]) for (i in 1:length(samp)) { del2 <- which(data[,2]==samp[i]) data <- data[-c(del2),] } colnames(data) <- c("Well", "Sample", "SNP", "Call", "NumMA", "Beta1_Const", "Beta2_Cons","Beta1_Cal", "Beta2_Cal") ind <- (length(data[,1])/6)-1 data[,5] <- as.numeric(data[,5]) for (i in 0:ind) { data[1+6*i,6] <- -4.8727 data[2+6*i,6] <- 0.5951 data[3+6*i,6] <- -1.5257 data[4+6*i,6] <- 1.1536 data[5+6*i,6] <- -0.5253 data[6+6*i,6] <- 0.4431 data[1+6*i,7] <- -1.99 data[2+6*i,7] <- 0.6872 data[3+6*i,7] <- -0.7363 data[4+6*i,7] <- 1.052 data[5+6*i,7] <- -0.0077 data[6+6*i,7] <- 0.259 } data[,8] <- data[,5]*data[,6] data[,9] <- data[,5]*data[,7] ## Setup and Calculate probabilities prob <- matrix(NA,ind+1,8) colnames(prob) <- c("Sample", "Alpha1_Cal", "Alpha2_Cal", "Prob_Blue", "Prob_Interm", "Prob_Brown", "Predicted_Colour_0.5", "Predicted_Colour_0.7") for (i in 0:ind) { prob[i+1, 1] <- data[1+6*i,2] } #Calculate probabilities alpha1 <- 3.8402 alpha2 <- 0.372 for (i in 0:ind) { prob[1+i,2] <- exp(alpha1+sum(data[(1+6*i):(6+6*i),8])) prob[1+i,3] <- exp(alpha2+sum(data[(1+6*i):(6+6*i),9])) } prob[,4] <- prob[,2]/(1+prob[,2]+prob[,3]) ## Prob Blue prob[,5] <- prob[,3]/(1+prob[,2]+prob[,3]) ## Prob Interm prob[,6] <- 1-prob[,4]-prob[,5] ## Prob Brown prob <- as.data.frame(prob) for (i in 1:(ind+1)) { if (prob[i,4]>=0.5) {prob[i,7]<-"Blue"} if (prob[i,5]>=0.5) {prob[i,7]<-"Intermediate"} if (prob[i,6]>=0.5) {prob[i,7]<-"Brown"} if (max(prob[i,4:6])<0.5) {prob[i,7]<-"Inconclusive"} if (prob[i,4]>=0.7) {prob[i,8]<-"Blue"} if (prob[i,5]>=0.7) {prob[i,8]<-"Intermediate"} if (prob[i,6]>=0.7) {prob[i,8]<-"Brown"} if (max(prob[i,4:6])<0.7) {prob[i,8]<-"Inconclusive"} } prob[,2:6] <- round(prob[,2:6],3) ## Export write.csv(prob, file = "EasyECD_results.csv") ##Accuracy calculations # PPV<-tp/(tp+fp) # NPV<-tn/(fn+tn) # sensitivity<-tp/(tp+fn) # specificity<-tn/(fp+tn) ## Update dlg <- tktoplevel() tkwm.title(dlg,"EasyECD") labelText <- tclVar("EasyECD has generated a comprehensive overview of the data that you loaded. \n\n Your output files have been successfully saved to your chosen directory. \n\nEnjoy your results and keep an eye out for new updates to EasyECD @ CRANIOFACIALidentification.com! \n") label1 <- tklabel(dlg,text=tclvalue(labelText)) tkconfigure(label1,textvariable=labelText) tkgrid(label1) Cont.but <- tkbutton(dlg,text="Close EasyECD",command=function()tkdestroy(dlg)) tkgrid(Cont.but) tkfocus(dlg) tkbind(dlg, "", function() {tkgrab.release(dlg)}) tkwait.window(dlg) rm(list = ls()) print("EasyECD analysis complete") } ## End of EasyECD