Analysis of GWAS results of S. pimpinellifolium screen performed in The Plant Accelerator
The data for this analysis was collected by Dr. Mitchell Morton at The Plant Accelerator in Adelaide, Australia. The data was analyzed as described in Dr. Morton’s PhD Thesis (KAUST). The data was subsequently used for the Genome Wide Association, using the ASReml script (similar to the method performed in (Awlia et al., 2021)[https://onlinelibrary.wiley.com/doi/10.1111/tpj.15310]). The GWAS analysis was performed by M.M.Julkowska at Boyce Thompson Institute. The raw data output of the GWAS script, as well as Manhattan and QQ Plots can be found at 
.
This analysis is meant to select the loci that would be of interest for subsequent studies - including genetic validation and such.
Libraries necessary for this analyses are:
#install.packages("qqman")
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("qvalue")
library(qqman)## ## For example usage please run: vignette('qqman')## ## Citation appreciated but not required:## Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731, https://doi.org/10.21105/joss.00731.## library(qvalue)
library(cowplot)
library(reshape2)GWAS plots - Manhattan and QQ plots
We need to generate the Manhattan plots from the Rdata files. We will do it by following the steps:
First - let’s load the data that we need:
setwd("/Users/magda/Dropbox/DataAndAnalysis/Tomato/ASReml/TPA_all_ASReml/")
my_files <- list.files(pattern = ".rda")
SNP <- read.csv("SNP_freq_TPA_Mitch.csv")
head(SNP)load(my_files[1])
head(results)Now - let’s prepare the table for graphing and determine the Bonferroni threshold (p value / # SNPs used)
# Prepare table for graphing:
args <- commandArgs(trailingOnly=TRUE)
results$Pval <-as.numeric(as.character(results$Pval))
dim(results)## [1] 708545 5head(results)a <- nrow(results)
b <- 0.05/a
b## [1] 7.056715e-08a## [1] 708545Awesome - now let’s plot the graphs:
title <- as.character(my_files[1])
title <- gsub("_gwas2029.rda", "", title)
# tiff(paste(title, "_Manhattan_plot.tiff"))
manhattan(results, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title,
ylim=c(0,12), col = c("coral3", "cornsilk4"), suggestiveline=FALSE, genomewideline=-log10(b))
# dev.off()
# tiff(paste(title, "_QQplot_MAF0.01.tiff"))
qq(results$Pval, main = title, cex = 1.5, las = 1)
# dev.off()In case you are interested in Manhattan / QQ plots with increased MAF (Minor Allele Frequency) or MAC (Minor Allele Count), we would need to fuse the results table with the SNP frequency data, subset the data for desired MAF/MAC, and adjust the Bonferroni threshold. Since we are considering only the SNPs with MAC > 2, let’s fo the same plot for MAC > 2:
head(results)results2 <- merge(SNP, results)
results2dim(results)## [1] 708545 5dim(SNP)## [1] 708545 3dim(results2)## [1] 708545 7After we checked that the fusion worked well - let’s calculate Bonferroni threshold:
results_new <- subset(results2, results2$MAC > 2)
a <- nrow(results_new)
b <- 0.05/a
b## [1] 1.303553e-07a## [1] 383567And now - let’s recalculate the graphs:
manhattan(results_new, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title,
ylim=c(0,12), col = c("coral3", "cornsilk4"), suggestiveline=FALSE, genomewideline=-log10(b))
# dev.off()
# tiff(paste(title, "_QQplot_MAF0.01.tiff"))
qq(results_new$Pval, main = title, cex = 1.5, las = 1)
# dev.off()Let’s loop this last bit - to generate the graphs for all traits with MAC > 2
# for(i in 1:length(my_files)){
# load(my_files[i])
# args <- commandArgs(trailingOnly=TRUE)
# results$Pval <-as.numeric(as.character(results$Pval))
# results2 <- merge(SNP, results)
# results_new <- subset(results2, results2$MAC > 2)
# title <- as.character(my_files[i])
# title <- gsub("_gwas2029.rda", "", title)
# tiff(paste(title, "_Manhattan_MAC3.tiff"))
# manhattan(results_new, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title, ylim=c(0,12),
# col = c("coral3", "cornsilk4"), suggestiveline=FALSE, genomewideline=-log10(b))
# dev.off()
#
# tiff(paste(title, "_QQplot_MAC3.tiff"))
# qq(results_new$Pval, main = title, cex = 1.5, las = 1)
# dev.off()
# }Done.
Loci identification
Let’s now move to identifying the loci - since from the Manhattan plots it is incredibly hard to see how many SNPs are in individual “trains”. For this - we will first generate the file with semi-significant loci - where -log10(pvalue) > 5. Just to be generous - we can later switch up to Bonferroni threshold.
The only traits that I would like to exclude - since the QQ plots look rubbish - are:
to_exclude <- c("WaterNa_me_SALT_gwas2029.rda", "STI2.AGR_gwas2029.rda", "STI1.RGR_gwas2029.rda",
"SFM_me_SALT_gwas2029.rda", "SDM_me_SALT_gwas2029.rda", "NaKratio_me_SALT_gwas2029.rda",
"Area.smooth.TUE.p3top1.BLUE_SALT_gwas2029.rda", "Area.smooth.TUE.p3top1.BLUE_CONTROL_gwas2029.rda",
"Area.smooth.AGR.p3top1.BLUE_SALT_gwas2029.rda", "Area.smooth.AGR.6to9.BLUE_SALT_gwas2029.rda",
"Area.smooth.AGR.6to9.BLUE_CONTROL_gwas2029.rda", "Area.smooth.AGR.0to9.BLUE_SALT_gwas2029.rda",
"Area.smooth.AGR.0to5.BLUE_SALT_gwas2029.rda", "STI1.RGR_gwas2029.rda")Let’s exclude these traits here:
my_files2 <- subset(my_files, !(my_files %in% to_exclude))
my_files2## [1] "Area.smooth.AGR.0to5.BLUE_CONTROL_gwas2029.rda"
## [2] "Area.smooth.AGR.0to9.BLUE_CONTROL_gwas2029.rda"
## [3] "Area.smooth.AGR.10to14.BLUE_CONTROL_gwas2029.rda"
## [4] "Area.smooth.AGR.10to14.BLUE_SALT_gwas2029.rda"
## [5] "Area.smooth.AGR.p3top1.BLUE_CONTROL_gwas2029.rda"
## [6] "Area.smooth.RGR.0to5.BLUE_CONTROL_gwas2029.rda"
## [7] "Area.smooth.RGR.0to5.BLUE_SALT_gwas2029.rda"
## [8] "Area.smooth.RGR.0to9.BLUE_CONTROL_gwas2029.rda"
## [9] "Area.smooth.RGR.0to9.BLUE_SALT_gwas2029.rda"
## [10] "Area.smooth.RGR.10to14.BLUE_CONTROL_gwas2029.rda"
## [11] "Area.smooth.RGR.10to14.BLUE_SALT_gwas2029.rda"
## [12] "Area.smooth.RGR.6to9.BLUE_CONTROL_gwas2029.rda"
## [13] "Area.smooth.RGR.6to9.BLUE_SALT_gwas2029.rda"
## [14] "Area.smooth.RGR.p3top1.BLUE_CONTROL_gwas2029.rda"
## [15] "Area.smooth.RGR.p3top1.BLUE_SALT_gwas2029.rda"
## [16] "Area.smooth.TUE.0to5.BLUE_CONTROL_gwas2029.rda"
## [17] "Area.smooth.TUE.0to5.BLUE_SALT_gwas2029.rda"
## [18] "Area.smooth.TUE.0to9.BLUE_CONTROL_gwas2029.rda"
## [19] "Area.smooth.TUE.0to9.BLUE_SALT_gwas2029.rda"
## [20] "Area.smooth.TUE.10to14.BLUE_CONTROL_gwas2029.rda"
## [21] "Area.smooth.TUE.10to14.BLUE_SALT_gwas2029.rda"
## [22] "Area.smooth.TUE.6to9.BLUE_CONTROL_gwas2029.rda"
## [23] "Area.smooth.TUE.6to9.BLUE_SALT_gwas2029.rda"
## [24] "NaKratio_me_CONTROL_gwas2029.rda"
## [25] "SDM_me_CONTROL_gwas2029.rda"
## [26] "SFM_me_CONTROL_gwas2029.rda"
## [27] "STI1.AGR_gwas2029.rda"
## [28] "STI2.RGR_gwas2029.rda"
## [29] "STI3.AGR_gwas2029.rda"
## [30] "STI3.RGR_gwas2029.rda"
## [31] "TissueK_me_CONTROL_gwas2029.rda"
## [32] "TissueK_me_SALT_gwas2029.rda"
## [33] "TissueNa_me_CONTROL_gwas2029.rda"
## [34] "TissueNa_me_SALT_gwas2029.rda"
## [35] "Transpiration.smooth.rate.0to5.BLUE_CONTROL_gwas2029.rda"
## [36] "Transpiration.smooth.rate.0to5.BLUE_SALT_gwas2029.rda"
## [37] "Transpiration.smooth.rate.0to9.BLUE_CONTROL_gwas2029.rda"
## [38] "Transpiration.smooth.rate.0to9.BLUE_SALT_gwas2029.rda"
## [39] "Transpiration.smooth.rate.10to14.BLUE_CONTROL_gwas2029.rda"
## [40] "Transpiration.smooth.rate.10to14.BLUE_SALT_gwas2029.rda"
## [41] "Transpiration.smooth.rate.6to9.BLUE_CONTROL_gwas2029.rda"
## [42] "Transpiration.smooth.rate.6to9.BLUE_SALT_gwas2029.rda"
## [43] "Transpiration.smooth.rate.p3top1.BLUE_CONTROL_gwas2029.rda"
## [44] "Transpiration.smooth.rate.p3top1.BLUE_SALT_gwas2029.rda"
## [45] "WaterK_me_CONTROL_gwas2029.rda"
## [46] "WaterK_me_SALT_gwas2029.rda"
## [47] "WaterNa_me_CONTROL_gwas2029.rda"length(my_files)## [1] 60length(my_files2)## [1] 47OK - looks good. So now - I would like to load the files one by one, and save all significant SNPs into one file:
i=1
load(my_files2[i])
args <- commandArgs(trailingOnly=TRUE)
results$Pval <-as.numeric(as.character(results$Pval))
title <- as.character(my_files[i])
title <- gsub("_gwas2029.rda", "", title)
results$trait <- title
results$LOD <- -log10(results$Pval)
head(results)LOD5_SNPs <- subset(results, results$LOD > 5)
LOD5_SNPsGreat - it seems to work well. Now - let’s loop it for all the other files:
for(i in 2:length(my_files2)){
load(my_files2[i])
args <- commandArgs(trailingOnly=TRUE)
results$Pval <-as.numeric(as.character(results$Pval))
title <- as.character(my_files[i])
title <- gsub("_gwas2029.rda", "", title)
results$trait <- title
results$LOD <- -log10(results$Pval)
head(results)
temporary <- subset(results, results$LOD > 5)
LOD5_SNPs <- rbind(LOD5_SNPs, temporary)
}
dim(LOD5_SNPs)## [1] 5412 7LOD5_SNPsunique(LOD5_SNPs$trait)## [1] "Area.smooth.AGR.0to5.BLUE_CONTROL" "Area.smooth.AGR.0to5.BLUE_SALT"
## [3] "Area.smooth.AGR.0to9.BLUE_CONTROL" "Area.smooth.AGR.0to9.BLUE_SALT"
## [5] "Area.smooth.AGR.10to14.BLUE_SALT" "Area.smooth.AGR.6to9.BLUE_CONTROL"
## [7] "Area.smooth.AGR.6to9.BLUE_SALT" "Area.smooth.AGR.p3top1.BLUE_CONTROL"
## [9] "Area.smooth.RGR.0to5.BLUE_CONTROL" "Area.smooth.RGR.0to5.BLUE_SALT"
## [11] "Area.smooth.RGR.0to9.BLUE_CONTROL" "Area.smooth.RGR.0to9.BLUE_SALT"
## [13] "Area.smooth.RGR.10to14.BLUE_CONTROL" "Area.smooth.RGR.6to9.BLUE_CONTROL"
## [15] "Area.smooth.RGR.6to9.BLUE_SALT" "Area.smooth.RGR.p3top1.BLUE_SALT"
## [17] "Area.smooth.TUE.0to5.BLUE_CONTROL" "Area.smooth.TUE.0to5.BLUE_SALT"
## [19] "Area.smooth.TUE.0to9.BLUE_CONTROL" "Area.smooth.TUE.10to14.BLUE_CONTROL"
## [21] "Area.smooth.TUE.10to14.BLUE_SALT" "Area.smooth.TUE.6to9.BLUE_SALT"
## [23] "Area.smooth.TUE.p3top1.BLUE_CONTROL" "Area.smooth.TUE.p3top1.BLUE_SALT"
## [25] "NaKratio_me_SALT" "SDM_me_SALT"
## [27] "SFM_me_CONTROL" "SFM_me_SALT"
## [29] "STI1.AGR" "STI1.RGR"
## [31] "STI2.AGR" "STI2.RGR"
## [33] "STI3.AGR" "STI3.RGR"
## [35] "TissueK_me_CONTROL" "TissueK_me_SALT"
## [37] "TissueNa_me_SALT"because we would also like to consider in our SNP selections other features - such as SNP Minor Allele Frequency - we should add it to this file now:
head(SNP)head(LOD5_SNPs)LOD5_SNPs$SNP %in% SNP$SNP## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [113] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [127] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [141] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [155] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [183] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [197] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [225] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [239] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [253] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [267] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [281] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [295] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [309] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [323] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [337] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [351] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [365] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [379] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [393] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [407] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [435] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [449] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [463] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [477] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [491] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [505] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [519] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [533] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [547] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [561] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [575] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [589] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [603] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [617] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [631] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [645] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [659] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [673] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [687] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [701] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [715] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [729] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [743] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [757] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [771] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [785] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [799] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [813] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [827] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [841] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [855] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [869] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [883] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [897] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [911] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [925] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [939] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [953] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [967] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [981] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [995] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1009] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1023] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1037] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1051] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1065] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1079] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1093] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1107] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1135] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1149] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1163] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1177] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1191] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1205] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1219] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1233] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1247] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1261] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1275] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1289] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1303] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1317] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1331] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1345] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1359] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1373] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1387] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1401] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1415] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1429] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1443] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1457] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1471] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1485] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1499] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1513] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1527] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1541] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1555] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1569] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1583] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1597] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1611] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1625] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1639] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1653] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1667] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1681] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1695] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1709] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1723] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1737] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1751] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1765] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1779] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1793] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1807] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1821] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1835] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1849] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1863] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1877] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1891] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1905] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1919] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1933] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1947] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1961] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1975] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1989] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2003] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2017] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2031] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2045] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2059] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2073] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2087] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2101] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2115] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2129] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2143] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2157] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2171] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2185] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2199] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2213] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2227] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2255] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2269] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2283] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2297] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2311] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2325] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2339] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2353] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2367] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2381] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2395] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2409] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2423] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2437] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2451] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2465] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2479] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2493] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2507] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2521] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2535] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2549] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2563] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2577] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2591] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2605] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2619] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2633] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2647] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2661] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2675] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2689] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2703] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2717] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2731] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2745] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2759] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2773] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2787] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2801] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2815] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2829] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2843] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2857] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2871] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2885] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2899] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2913] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2927] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2941] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2955] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2969] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2983] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2997] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3011] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3025] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3039] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3053] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3067] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3081] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3095] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3109] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3123] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3137] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3165] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3179] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3193] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3207] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3221] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3235] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3249] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3263] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3277] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3291] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3305] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3319] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3333] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3347] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3361] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3375] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3389] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3403] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3417] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3431] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3445] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3459] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3473] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3487] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3501] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3515] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3529] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3543] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3557] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3571] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3585] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3599] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3613] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3627] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3641] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3655] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3669] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3683] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3697] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3711] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3725] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3739] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3753] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3767] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3781] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3795] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3809] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3823] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3837] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3851] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3865] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3879] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3893] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3907] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3921] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3935] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3949] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3963] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3977] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3991] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4005] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4019] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4033] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4047] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4061] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4075] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4089] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4103] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4117] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4131] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4145] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4159] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4173] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4187] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4201] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4215] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4229] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4243] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4257] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4285] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4299] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4313] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4327] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4341] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4355] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4369] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4383] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4397] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4411] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4425] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4439] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4453] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4467] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4481] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4495] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4509] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4523] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4537] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4551] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4565] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4579] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4593] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4607] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4621] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4635] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4649] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4663] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4677] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4691] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4705] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4719] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4733] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4747] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4761] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4775] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4789] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4803] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4817] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4831] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4845] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4859] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4873] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4887] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4901] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4915] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4929] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4943] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4957] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4971] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4985] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4999] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5013] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5027] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5041] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5055] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5069] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5083] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5097] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5111] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5125] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5139] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5153] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5167] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5195] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5209] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5223] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5237] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5251] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5265] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5279] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5293] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5307] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5321] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5335] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5349] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5363] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5377] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5391] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5405] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEBonf2<- merge(LOD5_SNPs, SNP, by="SNP")
dim(LOD5_SNPs)## [1] 5412 7dim(Bonf2)## [1] 5412 9head(Bonf2)The files have different dimentions - since we loaded ALL the SNPs now - even the ones that were NEVER found to be significant. For this purpose - we will just subset the file again for LOD score above 5:
Bonf3 <- subset(Bonf2, Bonf2$LOD > 5)
dim(Bonf3)## [1] 5412 9head(Bonf3)Looks much better. And now - I want to organize this file - so we can actually see where the interesting loci are. For this purpose - we need to order SNPs per Chromosome and Position:
Bonf4 <- Bonf3[order(Bonf3$Chr, Bonf3$Pos),]
head(Bonf4)Looks gravy!
Cool - since we said we are not going to consider associations represented only by 2 accessions - we will now subset all of our associations for MAC > 2:
Bonf5 <- subset(Bonf4, Bonf4$MAC > 2)
dim(Bonf5)## [1] 3225 9head(Bonf5)OK - now in order to go through the SNPs - we need to calculate the distance between the neighboring SNPs:
Bonf5$dist <- 0
for(i in 2:nrow(Bonf5)){
Bonf5$dist[i] <- Bonf5$Pos[i] - Bonf5$Pos[i-1]
}
head(Bonf5)mean(Bonf5$dist)## [1] 21663.31min(Bonf5$dist)## [1] -87043188write.csv(Bonf5, "TPA_LOD5_loci_Mitch.csv", row.names = FALSE)OK - looks like it worked well.
Now - in order to get more insight into what exciting associations are there, we will generate a table where we make an overview of the loci. We need following info: - 1st SNP - last SNP - how many SNPs within the locus - how many traits are within this locus - mean LOD - max LOD - mean MAC - max MAC
First - let’s generate empty table:
names <- c(text="Chr", "First_SNP", "Last_SNP", "NoSNPs", "NoTraits", "meanLOD", "maxLOD", "meanMAC", "maxMAC", "Locus_number")
Loci_overview <- data.frame()
for (k in names) Loci_overview[[k]] <- as.character()
Loci_overviewThen - we need to create a loop where we go through the SNPs and combine following rows, as long as the distance is reasonable (let’s say 10 000):
Loci_overview[1,1] <- Bonf5$Chr[1]
Loci_overview[1,2] <- Bonf5$Pos[1]
i <- 1
x <- Bonf5[1,]
repeat {
i = i+1
y <- Bonf5[i,]
if(y$dist[1] < 10000){
x <- rbind(x, y)
}
if(y$dist[1] > 10000){
break
}
}
Loci_overview[1,1]<- min(x$Chr)
Loci_overview[1,2]<- min(x$Pos)
Loci_overview[1,3]<- max(x$Pos)
Loci_overview[1,4]<- length(unique(x$SNP))
Loci_overview[1,5]<- length(unique(x$trait))
Loci_overview[1,6]<- mean(x$LOD)
Loci_overview[1,7]<- max(x$LOD)
Loci_overview[1,8]<- mean(x$MAC)
Loci_overview[1,9]<- max(x$MAC)
lc <- 1
Loci_overview[1,10] <- paste("Locus", lc)
Loci_overviewGreat - before we put it into a loop - let’s adjust just few of the parameters:
lc <- lc + 1
lc## [1] 2i## [1] 2x <- y
xOK - now - let’s loop it up!
# nrow(Bonf5)
# while (i < 10836) {
# repeat {
# i = i+1
# y <- Bonf5[i,]
# if(y[1,10] < 10000 & y[1,10] > -1){
# x <- rbind(x, y)
# }
# else{
# Loci_overview[lc,1]<- min(x$Chr)
# Loci_overview[lc,2]<- min(x$Pos)
# Loci_overview[lc,3]<- max(x$Pos)
# Loci_overview[lc,4]<- length(unique(x$SNP))
# Loci_overview[lc,5]<- length(unique(x$trait))
# Loci_overview[lc,6]<- mean(x$LOD)
# Loci_overview[lc,7]<- max(x$LOD)
# Loci_overview[lc,8]<- mean(x$MAC)
# Loci_overview[lc,9]<- max(x$MAC)
# Loci_overview[lc,10] <- paste("Locus", lc)
# lc <- lc + 1
# x <- y
# }
# }
# }
# Loci_overview
# write.csv(Loci_overview, "Mitch_Loci_overview.csv", row.names = FALSE)For whatever reason - the above is always ending with an error - so let’s hash it out and upload the resulting file through CSV
Loci_overview <- read.csv("Mitch_Loci_overview.csv")
Loci_overviewOK- looks good! So we have 2286 associations with LOD > 5.
Let’s have a look how many of these loci were actually Bonferroni significant:
Bonf_associations <- subset(Loci_overview, Loci_overview$maxLOD > -log10(b))
Bonf_associationsCool - so we have 345 associations that were Bonferroni significant, now let’s have a look at how many LOCI (associations with more than 1 SNP) we have among these Bonferroni guys:
tru_loci <- subset(Bonf_associations, Bonf_associations$NoSNPs > 1)
tru_locidim(tru_loci)## [1] 68 10write.csv(tru_loci, "TPA_Mitch_Bonf_loci.csv", row.names = FALSE)We have total of 68 loci!
Awesome - let’s have a look at the BEST loci:
tru_loci$NoSNPs <- as.numeric(as.character(tru_loci$NoSNPs))
tru_loci$NoTraits <- as.numeric(as.character(tru_loci$NoTraits))
tru_loci$First_SNP <- as.numeric(as.character(tru_loci$First_SNP))
tru_loci$Last_SNP <- as.numeric(as.character(tru_loci$Last_SNP))
best_loci <- subset(tru_loci, tru_loci$NoSNPs > 4)
best_loci <- subset(best_loci, best_loci$NoTraits > 3)
best_loci$dist <- best_loci$Last_SNP - best_loci$First_SNP
dim(best_loci)## [1] 4 11best_lociCool - we have 4 BEST loci. Let’s have a look at them individually:
library(ggpubr)## Loading required package: ggplot2##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legendfor(i in 1:nrow(best_loci)){
best_loci$Chr[i]
max(best_loci$Last_SNP)
Chr_i_want <- best_loci$Chr[i]
Pos_min <- min(best_loci$First_SNP[i])
Pos_max <- max(best_loci$Last_SNP[i])
temp <- subset(Bonf5, Bonf5$Chr == Chr_i_want)
temp <- subset(temp, temp$Pos > Pos_min)
temp <- subset(temp, temp$Pos < Pos_max)
temp
plot(ggscatter(temp, x = "Pos", y = "LOD", size = "MAC", color = "trait") + theme(legend.position = "right"))
}
Locus 1
Chr_temp <- as.numeric(as.character(best_loci[1,1]))
from_SNP <- best_loci[1,2]
to_SNP <- best_loci[1,3]
Chr_temp## [1] 1from_SNP## [1] 77167543to_SNP## [1] 77181398BestLocus_1 <- subset(Bonf5, Bonf5$Chr == Chr_temp)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos >= from_SNP)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos <= to_SNP)
BestLocus_1write.csv(BestLocus_1, "Best_Locus1.csv", row.names = FALSE)Locus 2
Chr_temp <- as.numeric(as.character(best_loci[2,1]))
from_SNP <- best_loci[2,2]
to_SNP <- best_loci[2,3]
Chr_temp## [1] 1from_SNP## [1] 77311777to_SNP## [1] 77330707BestLocus_1 <- subset(Bonf5, Bonf5$Chr == Chr_temp)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos >= from_SNP)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos <= to_SNP)
BestLocus_1write.csv(BestLocus_1, "Best_Locus2.csv", row.names = FALSE)Locus 3
Chr_temp <- as.numeric(as.character(best_loci[3,1]))
from_SNP <- best_loci[3,2]
to_SNP <- best_loci[3,3]
Chr_temp## [1] 1from_SNP## [1] 79863262to_SNP## [1] 79871863BestLocus_1 <- subset(Bonf5, Bonf5$Chr == Chr_temp)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos >= from_SNP)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos <= to_SNP)
BestLocus_1write.csv(BestLocus_1, "Best_Locus3.csv", row.names = FALSE)Locus 4
Chr_temp <- as.numeric(as.character(best_loci[4,1]))
from_SNP <- best_loci[4,2]
to_SNP <- best_loci[4,3]
Chr_temp## [1] 3from_SNP## [1] 4621370to_SNP## [1] 4653226BestLocus_1 <- subset(Bonf5, Bonf5$Chr == Chr_temp)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos >= from_SNP)
BestLocus_1 <- subset(BestLocus_1, BestLocus_1$Pos <= to_SNP)
BestLocus_1write.csv(BestLocus_1, "Best_Locus4.csv", row.names = FALSE)Further inspection of most interesting loci
OK - so now we have established which loci were associated by multiple traits - let’s have a look at their overall context and LD within the identified loci. We are NOT going to look at the associations with the traits that we previously established to have odd QQ plots (STI.RGR and SFM_me_SALT).
There are few traits of interest that have interesting loci:
my_files2## [1] "Area.smooth.AGR.0to5.BLUE_CONTROL_gwas2029.rda"
## [2] "Area.smooth.AGR.0to9.BLUE_CONTROL_gwas2029.rda"
## [3] "Area.smooth.AGR.10to14.BLUE_CONTROL_gwas2029.rda"
## [4] "Area.smooth.AGR.10to14.BLUE_SALT_gwas2029.rda"
## [5] "Area.smooth.AGR.p3top1.BLUE_CONTROL_gwas2029.rda"
## [6] "Area.smooth.RGR.0to5.BLUE_CONTROL_gwas2029.rda"
## [7] "Area.smooth.RGR.0to5.BLUE_SALT_gwas2029.rda"
## [8] "Area.smooth.RGR.0to9.BLUE_CONTROL_gwas2029.rda"
## [9] "Area.smooth.RGR.0to9.BLUE_SALT_gwas2029.rda"
## [10] "Area.smooth.RGR.10to14.BLUE_CONTROL_gwas2029.rda"
## [11] "Area.smooth.RGR.10to14.BLUE_SALT_gwas2029.rda"
## [12] "Area.smooth.RGR.6to9.BLUE_CONTROL_gwas2029.rda"
## [13] "Area.smooth.RGR.6to9.BLUE_SALT_gwas2029.rda"
## [14] "Area.smooth.RGR.p3top1.BLUE_CONTROL_gwas2029.rda"
## [15] "Area.smooth.RGR.p3top1.BLUE_SALT_gwas2029.rda"
## [16] "Area.smooth.TUE.0to5.BLUE_CONTROL_gwas2029.rda"
## [17] "Area.smooth.TUE.0to5.BLUE_SALT_gwas2029.rda"
## [18] "Area.smooth.TUE.0to9.BLUE_CONTROL_gwas2029.rda"
## [19] "Area.smooth.TUE.0to9.BLUE_SALT_gwas2029.rda"
## [20] "Area.smooth.TUE.10to14.BLUE_CONTROL_gwas2029.rda"
## [21] "Area.smooth.TUE.10to14.BLUE_SALT_gwas2029.rda"
## [22] "Area.smooth.TUE.6to9.BLUE_CONTROL_gwas2029.rda"
## [23] "Area.smooth.TUE.6to9.BLUE_SALT_gwas2029.rda"
## [24] "NaKratio_me_CONTROL_gwas2029.rda"
## [25] "SDM_me_CONTROL_gwas2029.rda"
## [26] "SFM_me_CONTROL_gwas2029.rda"
## [27] "STI1.AGR_gwas2029.rda"
## [28] "STI2.RGR_gwas2029.rda"
## [29] "STI3.AGR_gwas2029.rda"
## [30] "STI3.RGR_gwas2029.rda"
## [31] "TissueK_me_CONTROL_gwas2029.rda"
## [32] "TissueK_me_SALT_gwas2029.rda"
## [33] "TissueNa_me_CONTROL_gwas2029.rda"
## [34] "TissueNa_me_SALT_gwas2029.rda"
## [35] "Transpiration.smooth.rate.0to5.BLUE_CONTROL_gwas2029.rda"
## [36] "Transpiration.smooth.rate.0to5.BLUE_SALT_gwas2029.rda"
## [37] "Transpiration.smooth.rate.0to9.BLUE_CONTROL_gwas2029.rda"
## [38] "Transpiration.smooth.rate.0to9.BLUE_SALT_gwas2029.rda"
## [39] "Transpiration.smooth.rate.10to14.BLUE_CONTROL_gwas2029.rda"
## [40] "Transpiration.smooth.rate.10to14.BLUE_SALT_gwas2029.rda"
## [41] "Transpiration.smooth.rate.6to9.BLUE_CONTROL_gwas2029.rda"
## [42] "Transpiration.smooth.rate.6to9.BLUE_SALT_gwas2029.rda"
## [43] "Transpiration.smooth.rate.p3top1.BLUE_CONTROL_gwas2029.rda"
## [44] "Transpiration.smooth.rate.p3top1.BLUE_SALT_gwas2029.rda"
## [45] "WaterK_me_CONTROL_gwas2029.rda"
## [46] "WaterK_me_SALT_gwas2029.rda"
## [47] "WaterNa_me_CONTROL_gwas2029.rda"finalists <- c("Area.smooth.AGR.0to5.BLUE_CONTROL_gwas2029.rda", "Area.smooth.AGR.0to5.BLUE_SALT_gwas2029.rda", "Area.smooth.AGR.0to9.BLUE_SALT_gwas2029.rda",
"Area.smooth.AGR.10to14.BLUE_CONTROL_gwas2029.rda", "Area.smooth.TUE.10to14.BLUE_SALT_gwas2029.rda")Let’s start with the first trait
Let’s isolate significant SNPs within that locus from the rda results file - to observe the association pattern within that locus:
load(finalists[1])
finalists[1]## [1] "Area.smooth.AGR.0to5.BLUE_CONTROL_gwas2029.rda"head(results)results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)best_loci[1,]STI <- subset(results2, results2$Chr == best_loci[1,1])
STI <- subset(STI, STI$Pos > (min(best_loci[1,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[1,3]) + 50000))
STILocus1_tr1 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cornflowerblue", title = "Area.smooth.AGR.0to5.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus1_tr1
best_loci[2,]STI <- subset(results2, results2$Chr == best_loci[2,1])
STI <- subset(STI, STI$Pos > (min(best_loci[2,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[2,3]) + 50000))
STILocus2_tr1<- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cornflowerblue", title = "Area.smooth.AGR.0to5.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus2_tr1
best_loci[3,]STI <- subset(results2, results2$Chr == best_loci[3,1])
STI <- subset(STI, STI$Pos > (min(best_loci[3,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[3,3]) + 50000))
STILocus3_tr1 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cornflowerblue", title = "Area.smooth.AGR.0to5.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus3_tr1
best_loci[4,]STI <- subset(results2, results2$Chr == best_loci[4,1])
STI <- subset(STI, STI$Pos > (min(best_loci[4,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[4,3]) + 50000))
STILocus4_tr1 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cornflowerblue", title = "Area.smooth.AGR.0to5.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus4_tr1
load(finalists[2])
finalists[2]## [1] "Area.smooth.AGR.0to5.BLUE_SALT_gwas2029.rda"head(results)results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)best_loci[1,]STI <- subset(results2, results2$Chr == best_loci[1,1])
STI <- subset(STI, STI$Pos > (min(best_loci[1,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[1,3]) + 50000))
STILocus1_tr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "coral4", title = "Area.smooth.AGR.0to5.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus1_tr2
best_loci[2,]STI <- subset(results2, results2$Chr == best_loci[2,1])
STI <- subset(STI, STI$Pos > (min(best_loci[2,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[2,3]) + 50000))
STILocus2_tr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "coral4", title = "Area.smooth.AGR.0to5.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus2_tr2
best_loci[3,]STI <- subset(results2, results2$Chr == best_loci[3,1])
STI <- subset(STI, STI$Pos > (min(best_loci[3,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[3,3]) + 50000))
STILocus3_tr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "coral4", title = "Area.smooth.AGR.0to5.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus3_tr2
best_loci[4,]STI <- subset(results2, results2$Chr == best_loci[4,1])
STI <- subset(STI, STI$Pos > (min(best_loci[4,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[4,3]) + 50000))
STILocus4_tr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "coral4", title = "Area.smooth.AGR.0to5.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus4_tr2
load(finalists[3])
finalists[3]## [1] "Area.smooth.AGR.0to9.BLUE_SALT_gwas2029.rda"head(results)results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)best_loci[1,]STI <- subset(results2, results2$Chr == best_loci[1,1])
STI <- subset(STI, STI$Pos > (min(best_loci[1,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[1,3]) + 50000))
STILocus1_tr3 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "Area.smooth.AGR.0to9.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus1_tr3
best_loci[2,]STI <- subset(results2, results2$Chr == best_loci[2,1])
STI <- subset(STI, STI$Pos > (min(best_loci[2,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[2,3]) + 50000))
STILocus2_tr3 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "Area.smooth.AGR.0to9.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus2_tr3
best_loci[3,]STI <- subset(results2, results2$Chr == best_loci[3,1])
STI <- subset(STI, STI$Pos > (min(best_loci[3,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[3,3]) + 50000))
STILocus3_tr3 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "Area.smooth.AGR.0to9.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus3_tr3
best_loci[4,]STI <- subset(results2, results2$Chr == best_loci[4,1])
STI <- subset(STI, STI$Pos > (min(best_loci[4,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[4,3]) + 50000))
STILocus4_tr3 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "Area.smooth.AGR.0to9.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus4_tr3
load(finalists[4])
finalists[4]## [1] "Area.smooth.AGR.10to14.BLUE_CONTROL_gwas2029.rda"head(results)results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)best_loci[1,]STI <- subset(results2, results2$Chr == best_loci[1,1])
STI <- subset(STI, STI$Pos > (min(best_loci[1,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[1,3]) + 50000))
STILocus1_tr4 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cadetblue", title = "Area.smooth.AGR.10to14.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus1_tr4
best_loci[2,]STI <- subset(results2, results2$Chr == best_loci[2,1])
STI <- subset(STI, STI$Pos > (min(best_loci[2,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[2,3]) + 50000))
STILocus2_tr4 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cadetblue", title = "Area.smooth.AGR.10to14.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus2_tr4
best_loci[3,]STI <- subset(results2, results2$Chr == best_loci[3,1])
STI <- subset(STI, STI$Pos > (min(best_loci[3,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[3,3]) + 50000))
STILocus3_tr4 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cadetblue", title = "Area.smooth.AGR.10to14.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus3_tr4
best_loci[4,]STI <- subset(results2, results2$Chr == best_loci[4,1])
STI <- subset(STI, STI$Pos > (min(best_loci[4,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[4,3]) + 50000))
STILocus4_tr4 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "cadetblue", title = "Area.smooth.AGR.10to14.BLUE_CONTROL", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus4_tr4
load(finalists[5])
finalists[5]## [1] "Area.smooth.TUE.10to14.BLUE_SALT_gwas2029.rda"head(results)results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)best_loci[1,]STI <- subset(results2, results2$Chr == best_loci[1,1])
STI <- subset(STI, STI$Pos > (min(best_loci[1,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[1,3]) + 50000))
STILocus1_tr5 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "deeppink4", title = "Area.smooth.TUE.10to14.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus1_tr5
best_loci[2,]STI <- subset(results2, results2$Chr == best_loci[2,1])
STI <- subset(STI, STI$Pos > (min(best_loci[2,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[2,3]) + 50000))
STILocus2_tr5 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "deeppink4", title = "Area.smooth.TUE.10to14.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus2_tr5
best_loci[3,]STI <- subset(results2, results2$Chr == best_loci[3,1])
STI <- subset(STI, STI$Pos > (min(best_loci[3,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[3,3]) + 50000))
STILocus3_tr5 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "deeppink4", title = "Area.smooth.TUE.10to14.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus3_tr5
best_loci[4,]STI <- subset(results2, results2$Chr == best_loci[4,1])
STI <- subset(STI, STI$Pos > (min(best_loci[4,2]) - 50000))
STI <- subset(STI, STI$Pos < (max(best_loci[4,3]) + 50000))
STILocus4_tr5 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "deeppink4", title = "Area.smooth.TUE.10to14.BLUE_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Locus4_tr5
Calculating LD for most promising loci:
OK - so now we have narrowed it down to 4 loci - let’s calculate the actual LD for each of them. For this purpose - let’s load the Genotype information:
PimpGeno <- read.csv("Pimp_all_chromosomes.csv")
head(PimpGeno)dim(PimpGeno)## [1] 708546 205colnames(PimpGeno) <- PimpGeno[1,]
PimpGeno2 <- PimpGeno[2:708546,2:205]
head(PimpGeno2)unique(PimpGeno2$chrom)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"let’s subset to our locus of interest - starting with locus 1
locus <- subset(PimpGeno2, PimpGeno2$chrom == best_loci[1,1])
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (min(best_loci[1,2]) - 50000))
locus <- subset(locus, locus$pos < (max(best_loci[1,3]) + 50000))
head(locus)max(locus$pos)## [1] 77229599min(locus$pos)## [1] 77125406locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
head(locus2)locusT <- t(locus2)
dim(locusT)## [1] 193 71plot(locus$pos)
now we need to select only SNPs that show variance within our population - so more than one allele!
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1
}
}
colnames(good_locus)[1] <- colnames(locusT)[1]
head(good_locus)## S1_77125406 S1_77125450 S1_77127017 S1_77127149 S1_77127382 S1_77128151
## M186 "C" "A" "G" "C" "C" "C"
## M033 "C" "A" "G" "C" "C" "C"
## M062 "C" "A" "G" "C" "C" "C"
## M085 "C" "A" "G" "C" "C" "C"
## M102 "C" "A" "G" "C" "C" "C"
## M099 "C" "A" "G" "C" "C" "C"
## S1_77129226 S1_77131907 S1_77132962 S1_77133879 S1_77135338 S1_77135464
## M186 "C" "G" "G" "T" "A" "C"
## M033 "C" "G" "G" "T" "A" "C"
## M062 "C" "G" "G" "T" "A" "C"
## M085 "C" "G" "G" "T" "A" "T"
## M102 "C" "G" "G" "T" "A" "C"
## M099 "C" "G" "G" "T" "A" "C"
## S1_77137812 S1_77141507 S1_77144489 S1_77148751 S1_77148823 S1_77150662
## M186 "T" "T" "T" "G" "G" "C"
## M033 "T" "T" "T" "G" "G" "C"
## M062 "T" "T" "T" "G" "G" "C"
## M085 "T" "T" "T" "G" "G" "C"
## M102 "T" "T" "T" "G" "G" "C"
## M099 "T" "T" "T" "G" "G" "C"
## S1_77152211 S1_77153380 S1_77153640 S1_77154687 S1_77155851 S1_77156017
## M186 "T" "A" "C" "T" "C" "G"
## M033 "T" "A" "C" "T" "C" "G"
## M062 "T" "A" "C" "T" "C" "G"
## M085 "T" "A" "C" "T" "C" "G"
## M102 "T" "A" "C" "T" "C" "G"
## M099 "T" "A" "C" "T" "C" "G"
## S1_77162876 S1_77163984 S1_77164053 S1_77164203 S1_77164204 S1_77166348
## M186 "A" "A" "C" "C" "G" "A"
## M033 "A" "A" "C" "C" "G" "A"
## M062 "A" "A" "C" "C" "G" "A"
## M085 "A" "A" "T" "C" "G" "G"
## M102 "A" "A" "C" "C" "G" "A"
## M099 "A" "A" "C" "C" "G" "A"
## S1_77166757 S1_77167543 S1_77169248 S1_77172099 S1_77172136 S1_77172477
## M186 "G" "C" "C" "T" "C" "C"
## M033 "G" "C" "C" "T" "C" "C"
## M062 "G" "C" "C" "T" "C" "C"
## M085 "G" "C" "C" "T" "C" "C"
## M102 "G" "C" "C" "T" "C" "C"
## M099 "G" "C" "C" "T" "C" "C"
## S1_77172656 S1_77172983 S1_77173292 S1_77174599 S1_77174884 S1_77175097
## M186 "G" "A" "G" "T" "C" "A"
## M033 "G" "A" "G" "T" "C" "A"
## M062 "G" "A" "G" "T" "C" "A"
## M085 "G" "A" "G" "T" "C" "A"
## M102 "G" "A" "G" "T" "C" "A"
## M099 "G" "A" "G" "T" "C" "A"
## S1_77180581 S1_77181009 S1_77181328 S1_77181398 S1_77182430 S1_77183630
## M186 "G" "C" "C" "A" "A" "A"
## M033 "G" "C" "C" "A" "A" "A"
## M062 "G" "C" "A" "A" "A" "A"
## M085 "G" "C" "C" "A" "A" "A"
## M102 "G" "C" "C" "A" "A" "A"
## M099 "G" "C" "C" "A" "A" "A"
## S1_77185690 S1_77187172 S1_77187259 S1_77187538 S1_77191949 S1_77191981
## M186 "G" "A" "C" "C" "C" "T"
## M033 "G" "A" "C" "C" "C" "T"
## M062 "G" "A" "C" "C" "C" "T"
## M085 "G" "A" "C" "C" "C" "C"
## M102 "G" "A" "C" "C" "C" "T"
## M099 "G" "A" "C" "C" "C" "T"
## S1_77192410 S1_77192411 S1_77194647 S1_77195837 S1_77197111 S1_77198450
## M186 "G" "T" "G" "T" "G" "C"
## M033 "G" "T" "G" "T" "G" "C"
## M062 "G" "T" "G" "T" "G" "C"
## M085 "G" "T" "G" "T" "G" "C"
## M102 "G" "T" "G" "T" "G" "C"
## M099 "G" "T" "G" "T" "G" "C"
## S1_77210821 S1_77213576 S1_77219132 S1_77226356 S1_77227406 S1_77229599
## M186 "T" "T" "A" "C" "C" "T"
## M033 "T" "T" "A" "C" "C" "T"
## M062 "T" "T" "A" "C" "C" "T"
## M085 "T" "T" "A" "C" "C" "T"
## M102 "T" "T" "A" "C" "C" "T"
## M099 "T" "T" "A" "C" "C" "T"dim(locusT)## [1] 193 71dim(good_locus)## [1] 193 66ok - so the above is for calculating the LD between TWO loci. In order to make the correlation matrix for ALL loci of interest - we would need to make an empty matrix and put the values in there. Preferably with a loop!!!
library(genetics)## Loading required package: combinat##
## Attaching package: 'combinat'## The following object is masked from 'package:utils':
##
## combn## Loading required package: gdata## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.## ## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.##
## Attaching package: 'gdata'## The following object is masked from 'package:stats':
##
## nobs## The following object is masked from 'package:utils':
##
## object.size## The following object is masked from 'package:base':
##
## startsWith## Loading required package: gtools## Loading required package: MASS## Loading required package: mvtnorm## ## NOTE: THIS PACKAGE IS NOW OBSOLETE.## ## The R-Genetics project has developed an set of enhanced genetics## packages to replace 'genetics'. Please visit the project homepage## at http://rgenetics.org for informtion.## ##
## Attaching package: 'genetics'## The following objects are masked from 'package:base':
##
## %in%, as.factor, orderlibrary(LDheatmap)
# how many SNPs we have?
dim(good_locus)## [1] 193 66# 32 - so let;s make matrix that is 66 x 66
mymatrix_D <- matrix(, nrow = 66, ncol = 66)
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
#mymatrix_D
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`
}
}SNPdist <- gsub("S1_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
SNPdist## [1] 77125406 77125450 77127017 77127149 77127382 77128151 77129226 77131907
## [9] 77132962 77133879 77135338 77135464 77137812 77141507 77144489 77148751
## [17] 77148823 77150662 77152211 77153380 77153640 77154687 77155851 77156017
## [25] 77162876 77163984 77164053 77164203 77164204 77166348 77166757 77167543
## [33] 77169248 77172099 77172136 77172477 77172656 77172983 77173292 77174599
## [41] 77174884 77175097 77180581 77181009 77181328 77181398 77182430 77183630
## [49] 77185690 77187172 77187259 77187538 77191949 77191981 77192410 77192411
## [57] 77194647 77195837 77197111 77198450 77210821 77213576 77219132 77226356
## [65] 77227406 77229599length(SNPdist)## [1] 66dim(mymatrix_D)## [1] 66 66LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_77167543", "S1_77172099", "S1_77172136", "S1_77172477", "S1_77172656", "S1_77174884", "S1_77175097", "S1_77181398"), color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus1_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_77167543", "S1_77172099", "S1_77172136", "S1_77172477", "S1_77172656", "S1_77174884", "S1_77175097", "S1_77181398"), color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2Let’s save this before we forget into one figure with all the other zoomed-in plots from other interesting phenotypes:
library(cowplot)
mini_manhattans <- plot_grid(Locus1_tr1, Locus1_tr2, Locus1_tr3, Locus1_tr4,Locus1_tr5, ncol=5)
mini_manhattans
pdf("Mitch_Locus1_miniManhattan.pdf", height = 5, width = 30)
plot(mini_manhattans)
dev.off()## quartz_off_screen
## 2Locus 2:
locus <- subset(PimpGeno2, PimpGeno2$chrom == best_loci[2,1])
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (min(best_loci[2,2]) - 50000))
locus <- subset(locus, locus$pos < (max(best_loci[2,3]) + 50000))
head(locus)max(locus$pos)## [1] 77378720min(locus$pos)## [1] 77261914locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
head(locus2)locusT <- t(locus2)
dim(locusT)## [1] 193 64plot(locus$pos)
now we need to select only SNPs that show variance within our population - so more than one allele!
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1
}
}
colnames(good_locus)[1] <- colnames(locusT)[1]
head(good_locus)## S1_77261914 S1_77263740 S1_77266607 S1_77269804 S1_77273038 S1_77273990
## M186 "G" "A" "T" "T" "A" "T"
## M033 "G" "A" "T" "T" "A" "T"
## M062 "G" "A" "T" "T" "A" "T"
## M085 "A" "A" "T" "T" "A" "T"
## M102 "G" "A" "T" "T" "A" "T"
## M099 "G" "A" "T" "T" "A" "T"
## S1_77274707 S1_77278682 S1_77278922 S1_77279763 S1_77281698 S1_77283123
## M186 "G" "T" "A" "T" "G" "T"
## M033 "G" "T" "A" "T" "G" "T"
## M062 "G" "T" "A" "T" "G" "T"
## M085 "G" "T" "A" "T" "G" "T"
## M102 "G" "T" "A" "C" "G" "T"
## M099 "G" "T" "A" "T" "G" "T"
## S1_77286649 S1_77289095 S1_77290756 S1_77290878 S1_77291405 S1_77292185
## M186 "G" "G" "G" "A" "G" "T"
## M033 "G" "G" "G" "A" "G" "T"
## M062 "G" "G" "G" "A" "G" "T"
## M085 "G" "G" "G" "A" "G" "T"
## M102 "G" "G" "G" "A" "G" "T"
## M099 "G" "G" "G" "A" "G" "T"
## S1_77292373 S1_77293340 S1_77294114 S1_77298048 S1_77303244 S1_77308548
## M186 "C" "T" "C" "C" "T" "G"
## M033 "C" "T" "C" "C" "T" "G"
## M062 "C" "T" "C" "C" "T" "G"
## M085 "C" "T" "C" "C" "T" "G"
## M102 "C" "T" "C" "C" "T" "G"
## M099 "C" "T" "C" "C" "T" "G"
## S1_77309075 S1_77311777 S1_77313444 S1_77315990 S1_77317835 S1_77318711
## M186 "T" "A" "G" "G" "C" "C"
## M033 "T" "A" "G" "G" "C" "C"
## M062 "T" "A" "G" "G" "C" "C"
## M085 "T" "A" "G" "G" "C" "C"
## M102 "T" "A" "G" "G" "C" "C"
## M099 "T" "A" "G" "G" "C" "C"
## S1_77319515 S1_77320290 S1_77320338 S1_77320682 S1_77320749 S1_77321301
## M186 "G" "G" "G" "C" "C" "A"
## M033 "G" "G" "G" "C" "C" "A"
## M062 "G" "G" "G" "C" "C" "A"
## M085 "G" "G" "A" "C" "C" "A"
## M102 "G" "G" "G" "C" "C" "A"
## M099 "G" "G" "G" "C" "C" "A"
## S1_77321735 S1_77329652 S1_77330707 S1_77331119 S1_77332022 S1_77336243
## M186 "A" "G" "G" "T" "C" "A"
## M033 "A" "G" "G" "T" "C" "A"
## M062 "A" "G" "G" "T" "C" "A"
## M085 "A" "G" "G" "G" "C" "A"
## M102 "A" "G" "G" "T" "C" "A"
## M099 "A" "G" "G" "T" "C" "A"
## S1_77338855 S1_77342226 S1_77342547 S1_77342703 S1_77343111 S1_77343112
## M186 "G" "A" "C" "A" "A" "T"
## M033 "G" "A" "C" "A" "A" "T"
## M062 "G" "A" "C" "A" "A" "T"
## M085 "G" "A" "C" "A" "A" "T"
## M102 "G" "A" "C" "A" "A" "T"
## M099 "G" "A" "C" "A" "A" "T"
## S1_77348709 S1_77348899 S1_77356774 S1_77362711 S1_77363328 S1_77367456
## M186 "C" "C" "C" "G" "G" "T"
## M033 "C" "C" "C" "G" "G" "T"
## M062 "C" "C" "C" "G" "G" "T"
## M085 "C" "C" "C" "G" "G" "T"
## M102 "C" "C" "C" "G" "G" "T"
## M099 "C" "C" "C" "G" "G" "T"
## S1_77370767 S1_77370807 S1_77372952 S1_77374618 S1_77374662 S1_77374981
## M186 "G" "A" "G" "C" "C" "A"
## M033 "G" "A" "G" "C" "C" "A"
## M062 "G" "A" "G" "C" "C" "A"
## M085 "A" "A" "G" "C" "C" "A"
## M102 "G" "A" "G" "C" "C" "A"
## M099 "G" "A" "G" "C" "C" "A"
## S1_77375622 S1_77376630 S1_77376652 S1_77378720
## M186 "G" "T" "T" "T"
## M033 "G" "T" "T" "T"
## M062 "G" "T" "T" "T"
## M085 "G" "T" "T" "T"
## M102 "G" "T" "T" "T"
## M099 "G" "T" "T" "T"dim(locusT)## [1] 193 64dim(good_locus)## [1] 193 64ok - so the above is for calculating the LD between TWO loci. In order to make the correlation matrix for ALL loci of interest - we would need to make an empty matrix and put the values in there. Preferably with a loop!!!
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
#mymatrix_D
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`
}
}SNPdist <- gsub("S1_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
SNPdist## [1] 77261914 77263740 77266607 77269804 77273038 77273990 77274707 77278682
## [9] 77278922 77279763 77281698 77283123 77286649 77289095 77290756 77290878
## [17] 77291405 77292185 77292373 77293340 77294114 77298048 77303244 77308548
## [25] 77309075 77311777 77313444 77315990 77317835 77318711 77319515 77320290
## [33] 77320338 77320682 77320749 77321301 77321735 77329652 77330707 77331119
## [41] 77332022 77336243 77338855 77342226 77342547 77342703 77343111 77343112
## [49] 77348709 77348899 77356774 77362711 77363328 77367456 77370767 77370807
## [57] 77372952 77374618 77374662 77374981 77375622 77376630 77376652 77378720length(SNPdist)## [1] 64dim(mymatrix_D)## [1] 64 64LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_77311777", "S1_77313444", "S1_77315990", "S1_77319515", "S1_77320682", "S1_77329652", "S1_77330707"), color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus2_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_77311777", "S1_77313444", "S1_77315990", "S1_77319515", "S1_77320682", "S1_77329652", "S1_77330707"), color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2Let’s save this before we forget into one figure with all the other zoomed-in plots from other interesting phenotypes:
library(cowplot)
mini_manhattans <- plot_grid(Locus1_tr1, Locus1_tr2, Locus1_tr3, Locus1_tr4,Locus1_tr5, ncol=5)
mini_manhattans
pdf("Mitch_Locus2_miniManhattan.pdf", height = 5, width = 30)
plot(mini_manhattans)
dev.off()## quartz_off_screen
## 2Locus 3:
locus <- subset(PimpGeno2, PimpGeno2$chrom == best_loci[3,1])
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (min(best_loci[3,2]) - 50000))
locus <- subset(locus, locus$pos < (max(best_loci[3,3]) + 50000))
head(locus)max(locus$pos)## [1] 79921568min(locus$pos)## [1] 79814667locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
head(locus2)locusT <- t(locus2)
dim(locusT)## [1] 193 114plot(locus$pos)
now we need to select only SNPs that show variance within our population - so more than one allele!
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1
}
}
colnames(good_locus)[1] <- colnames(locusT)[1]
head(good_locus)## S1_79814667 S1_79816724 S1_79818053 S1_79819644 S1_79820061 S1_79820910
## M186 "C" "T" "A" "G" "A" "C"
## M033 "C" "T" "A" "G" "A" "C"
## M062 "C" "T" "A" "G" "A" "C"
## M085 "C" "T" "A" "G" "A" "C"
## M102 "C" "T" "A" "G" "A" "C"
## M099 "C" "T" "A" "G" "A" "C"
## S1_79820915 S1_79820936 S1_79821064 S1_79821475 S1_79822346 S1_79822368
## M186 "A" "T" "G" "A" "C" "C"
## M033 "A" "T" "G" "A" "C" "C"
## M062 "A" "T" "G" "A" "C" "C"
## M085 "A" "T" "G" "A" "C" "C"
## M102 "A" "T" "G" "A" "C" "C"
## M099 "A" "T" "G" "A" "C" "C"
## S1_79822701 S1_79825341 S1_79829963 S1_79830264 S1_79830399 S1_79830824
## M186 "G" "T" "C" "A" "G" "C"
## M033 "G" "T" "C" "A" "G" "C"
## M062 "G" "T" "C" "A" "G" "C"
## M085 "G" "T" "C" "A" "G" "C"
## M102 "G" "T" "C" "A" "G" "C"
## M099 "G" "T" "C" "A" "G" "C"
## S1_79831441 S1_79832755 S1_79834842 S1_79835019 S1_79835699 S1_79835884
## M186 "G" "G" "G" "C" "T" "A"
## M033 "G" "G" "G" "C" "T" "A"
## M062 "G" "G" "G" "C" "T" "A"
## M085 "G" "G" "G" "C" "T" "A"
## M102 "G" "G" "G" "C" "T" "A"
## M099 "G" "G" "G" "C" "T" "A"
## S1_79836004 S1_79836558 S1_79837198 S1_79838399 S1_79840108 S1_79840453
## M186 "T" "G" "T" "A" "A" "C"
## M033 "T" "G" "T" "A" "A" "C"
## M062 "T" "G" "T" "A" "A" "C"
## M085 "T" "G" "T" "A" "A" "C"
## M102 "T" "G" "T" "A" "A" "C"
## M099 "T" "G" "T" "A" "A" "C"
## S1_79840486 S1_79840582 S1_79840845 S1_79840943 S1_79849580 S1_79849618
## M186 "C" "C" "A" "T" "A" "C"
## M033 "C" "C" "A" "T" "A" "C"
## M062 "C" "C" "A" "T" "A" "C"
## M085 "C" "C" "A" "T" "A" "C"
## M102 "C" "C" "A" "T" "A" "C"
## M099 "C" "C" "A" "T" "A" "C"
## S1_79849770 S1_79850181 S1_79851040 S1_79855234 S1_79856019 S1_79856395
## M186 "C" "G" "A" "A" "A" "A"
## M033 "C" "G" "A" "A" "A" "A"
## M062 "C" "G" "A" "A" "A" "A"
## M085 "C" "G" "A" "A" "A" "A"
## M102 "C" "G" "A" "A" "A" "A"
## M099 "C" "G" "A" "A" "A" "A"
## S1_79856548 S1_79856686 S1_79856768 S1_79857255 S1_79859771 S1_79860052
## M186 "G" "A" "G" "T" "C" "G"
## M033 "G" "A" "G" "T" "C" "G"
## M062 "A" "A" "G" "T" "C" "G"
## M085 "G" "A" "G" "T" "C" "G"
## M102 "G" "A" "G" "T" "C" "G"
## M099 "G" "A" "G" "T" "C" "G"
## S1_79860310 S1_79861237 S1_79861268 S1_79862797 S1_79862862 S1_79863014
## M186 "A" "G" "C" "C" "C" "G"
## M033 "A" "G" "C" "C" "C" "G"
## M062 "A" "G" "C" "C" "C" "G"
## M085 "A" "G" "C" "C" "C" "G"
## M102 "A" "G" "C" "C" "C" "G"
## M099 "A" "G" "C" "C" "C" "G"
## S1_79863262 S1_79863620 S1_79864493 S1_79865185 S1_79865461 S1_79867264
## M186 "T" "T" "G" "C" "C" "A"
## M033 "T" "T" "G" "C" "C" "A"
## M062 "T" "T" "G" "C" "C" "A"
## M085 "T" "T" "G" "C" "C" "A"
## M102 "T" "T" "G" "C" "C" "A"
## M099 "T" "T" "G" "C" "C" "A"
## S1_79869569 S1_79870546 S1_79870562 S1_79871863 S1_79872042 S1_79872443
## M186 "C" "G" "A" "T" "T" "A"
## M033 "C" "G" "G" "T" "T" "A"
## M062 "C" "G" "A" "T" "T" "A"
## M085 "C" "G" "A" "T" "T" "A"
## M102 "C" "G" "A" "T" "T" "A"
## M099 "C" "G" "A" "T" "T" "A"
## S1_79874637 S1_79874651 S1_79876468 S1_79879230 S1_79879232 S1_79879578
## M186 "G" "C" "G" "C" "A" "G"
## M033 "G" "C" "G" "C" "A" "G"
## M062 "G" "C" "G" "C" "A" "G"
## M085 "G" "C" "G" "C" "A" "G"
## M102 "G" "C" "G" "C" "A" "G"
## M099 "G" "C" "G" "C" "A" "G"
## S1_79879712 S1_79881115 S1_79881360 S1_79881487 S1_79882057 S1_79882192
## M186 "G" "T" "G" "A" "A" "A"
## M033 "G" "T" "G" "A" "A" "A"
## M062 "G" "T" "G" "A" "A" "A"
## M085 "G" "T" "G" "A" "A" "A"
## M102 "G" "T" "G" "A" "A" "A"
## M099 "G" "T" "G" "A" "A" "A"
## S1_79882329 S1_79883394 S1_79888514 S1_79889159 S1_79891488 S1_79892519
## M186 "T" "A" "C" "T" "C" "A"
## M033 "T" "T" "C" "T" "C" "G"
## M062 "T" "T" "C" "T" "C" "G"
## M085 "T" "T" "C" "T" "C" "G"
## M102 "T" "T" "C" "T" "C" "G"
## M099 "T" "T" "C" "T" "C" "G"
## S1_79893043 S1_79893358 S1_79893449 S1_79894132 S1_79894668 S1_79895624
## M186 "T" "G" "T" "C" "T" "T"
## M033 "T" "G" "T" "C" "T" "T"
## M062 "T" "G" "T" "T" "T" "T"
## M085 "T" "G" "T" "C" "T" "T"
## M102 "T" "G" "T" "C" "T" "T"
## M099 "T" "G" "T" "C" "T" "T"
## S1_79896259 S1_79896990 S1_79897929 S1_79898654 S1_79901277 S1_79902546
## M186 "G" "G" "A" "G" "T" "G"
## M033 "G" "G" "A" "G" "T" "G"
## M062 "G" "G" "A" "G" "T" "G"
## M085 "G" "G" "A" "G" "T" "G"
## M102 "G" "G" "A" "G" "T" "G"
## M099 "G" "G" "A" "G" "T" "G"
## S1_79902815 S1_79903685 S1_79903717 S1_79903940 S1_79905989 S1_79907836
## M186 "C" "A" "C" "G" "C" "T"
## M033 "C" "A" "C" "G" "C" "T"
## M062 "C" "A" "C" "G" "C" "T"
## M085 "C" "A" "C" "G" "C" "T"
## M102 "C" "A" "C" "G" "C" "T"
## M099 "C" "A" "C" "G" "C" "T"
## S1_79911152 S1_79913140 S1_79917224 S1_79917734 S1_79920821 S1_79921488
## M186 "C" "A" "C" "C" "C" "G"
## M033 "C" "A" "C" "C" "C" "G"
## M062 "C" "A" "C" "C" "C" "G"
## M085 "C" "A" "C" "C" "C" "G"
## M102 "C" "A" "C" "C" "C" "G"
## M099 "C" "A" "C" "C" "C" "G"
## S1_79921568
## M186 "G"
## M033 "G"
## M062 "G"
## M085 "G"
## M102 "G"
## M099 "G"dim(locusT)## [1] 193 114dim(good_locus)## [1] 193 109ok - so the above is for calculating the LD between TWO loci. In order to make the correlation matrix for ALL loci of interest - we would need to make an empty matrix and put the values in there. Preferably with a loop!!!
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
#mymatrix_D
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`
}
}SNPdist <- gsub("S1_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
SNPdist## [1] 79814667 79816724 79818053 79819644 79820061 79820910 79820915 79820936
## [9] 79821064 79821475 79822346 79822368 79822701 79825341 79829963 79830264
## [17] 79830399 79830824 79831441 79832755 79834842 79835019 79835699 79835884
## [25] 79836004 79836558 79837198 79838399 79840108 79840453 79840486 79840582
## [33] 79840845 79840943 79849580 79849618 79849770 79850181 79851040 79855234
## [41] 79856019 79856395 79856548 79856686 79856768 79857255 79859771 79860052
## [49] 79860310 79861237 79861268 79862797 79862862 79863014 79863262 79863620
## [57] 79864493 79865185 79865461 79867264 79869569 79870546 79870562 79871863
## [65] 79872042 79872443 79874637 79874651 79876468 79879230 79879232 79879578
## [73] 79879712 79881115 79881360 79881487 79882057 79882192 79882329 79883394
## [81] 79888514 79889159 79891488 79892519 79893043 79893358 79893449 79894132
## [89] 79894668 79895624 79896259 79896990 79897929 79898654 79901277 79902546
## [97] 79902815 79903685 79903717 79903940 79905989 79907836 79911152 79913140
## [105] 79917224 79917734 79920821 79921488 79921568length(SNPdist)## [1] 109dim(mymatrix_D)## [1] 109 109LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_79867264", "S1_79869569", "S1_79870562", "S1_79871863"), color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus3_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S1_79867264", "S1_79869569", "S1_79870562", "S1_79871863"), color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2Let’s save this before we forget into one figure with all the other zoomed-in plots from other interesting phenotypes:
library(cowplot)
mini_manhattans <- plot_grid(Locus1_tr1, Locus1_tr2, Locus1_tr3, Locus1_tr4,Locus1_tr5, ncol=5)
mini_manhattans
pdf("Mitch_Locus3_miniManhattan.pdf", height = 5, width = 30)
plot(mini_manhattans)
dev.off()## quartz_off_screen
## 2Locus 4:
locus <- subset(PimpGeno2, PimpGeno2$chrom == best_loci[4,1])
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (min(best_loci[4,2]) - 50000))
locus <- subset(locus, locus$pos < (max(best_loci[4,3]) + 50000))
head(locus)max(locus$pos)## [1] 4690853min(locus$pos)## [1] 4571381locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
head(locus2)locusT <- t(locus2)
dim(locusT)## [1] 193 143plot(locus$pos)
now we need to select only SNPs that show variance within our population - so more than one allele!
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1
}
}
colnames(good_locus)[1] <- colnames(locusT)[1]
head(good_locus)## S3_4571381 S3_4571758 S3_4574165 S3_4575341 S3_4582900 S3_4582915
## M186 "T" "G" "A" "T" "G" "A"
## M033 "T" "G" "A" "T" "G" "A"
## M062 "T" "G" "A" "T" "G" "A"
## M085 "T" "G" "A" "T" "G" "A"
## M102 "T" "G" "A" "T" "G" "A"
## M099 "T" "G" "A" "T" "G" "A"
## S3_4583881 S3_4585421 S3_4585728 S3_4586470 S3_4586694 S3_4586731
## M186 "A" "C" "A" "T" "C" "A"
## M033 "A" "C" "A" "T" "C" "A"
## M062 "A" "C" "A" "T" "C" "A"
## M085 "A" "C" "A" "T" "C" "A"
## M102 "A" "C" "A" "T" "C" "A"
## M099 "A" "C" "A" "T" "C" "A"
## S3_4588214 S3_4588391 S3_4590830 S3_4590854 S3_4590926 S3_4590930
## M186 "T" "G" "T" "A" "G" "G"
## M033 "T" "G" "C" "A" "G" "G"
## M062 "T" "G" "C" "A" "G" "G"
## M085 "T" "G" "C" "A" "G" "G"
## M102 "T" "G" "C" "A" "G" "G"
## M099 "T" "G" "C" "A" "G" "G"
## S3_4592328 S3_4592891 S3_4595458 S3_4595932 S3_4596024 S3_4598149
## M186 "T" "G" "T" "G" "T" "C"
## M033 "T" "G" "T" "G" "T" "C"
## M062 "T" "G" "T" "G" "T" "C"
## M085 "T" "G" "T" "G" "T" "C"
## M102 "T" "G" "T" "G" "T" "C"
## M099 "T" "G" "T" "G" "T" "C"
## S3_4598569 S3_4599668 S3_4599702 S3_4599974 S3_4600723 S3_4602668
## M186 "A" "A" "A" "C" "C" "A"
## M033 "A" "A" "G" "C" "C" "A"
## M062 "A" "A" "G" "C" "C" "A"
## M085 "A" "A" "G" "C" "C" "A"
## M102 "A" "A" "G" "C" "C" "A"
## M099 "A" "A" "G" "C" "C" "A"
## S3_4602980 S3_4603950 S3_4604326 S3_4604859 S3_4605023 S3_4605188
## M186 "A" "T" "A" "A" "G" "C"
## M033 "A" "T" "A" "A" "G" "C"
## M062 "A" "T" "A" "A" "G" "C"
## M085 "A" "T" "A" "A" "G" "C"
## M102 "A" "T" "A" "A" "G" "C"
## M099 "A" "T" "A" "A" "G" "C"
## S3_4605315 S3_4605467 S3_4605524 S3_4605702 S3_4605865 S3_4606726
## M186 "A" "A" "C" "A" "T" "A"
## M033 "A" "A" "C" "A" "T" "A"
## M062 "A" "A" "C" "A" "T" "A"
## M085 "A" "A" "C" "A" "T" "A"
## M102 "A" "A" "C" "A" "T" "A"
## M099 "A" "A" "C" "A" "T" "A"
## S3_4607110 S3_4610035 S3_4610365 S3_4610928 S3_4611231 S3_4611653
## M186 "T" "T" "G" "T" "T" "G"
## M033 "T" "T" "G" "T" "T" "G"
## M062 "T" "T" "G" "T" "T" "G"
## M085 "T" "T" "G" "T" "T" "G"
## M102 "T" "T" "G" "T" "T" "G"
## M099 "T" "T" "G" "T" "T" "G"
## S3_4613096 S3_4617940 S3_4619977 S3_4621370 S3_4621908 S3_4622076
## M186 "A" "C" "A" "G" "A" "A"
## M033 "A" "C" "A" "G" "A" "A"
## M062 "A" "C" "A" "G" "A" "A"
## M085 "A" "C" "A" "G" "A" "A"
## M102 "A" "C" "A" "G" "A" "A"
## M099 "A" "C" "A" "G" "A" "A"
## S3_4623677 S3_4624465 S3_4625544 S3_4625558 S3_4625816 S3_4629034
## M186 "T" "T" "C" "C" "T" "A"
## M033 "T" "T" "C" "C" "T" "A"
## M062 "T" "T" "C" "C" "T" "A"
## M085 "T" "T" "C" "C" "T" "A"
## M102 "T" "T" "C" "C" "T" "A"
## M099 "T" "T" "C" "C" "T" "A"
## S3_4630077 S3_4630332 S3_4631533 S3_4632496 S3_4632635 S3_4634801
## M186 "T" "A" "T" "C" "C" "C"
## M033 "T" "A" "T" "C" "C" "C"
## M062 "T" "A" "T" "C" "C" "C"
## M085 "T" "A" "T" "C" "C" "C"
## M102 "T" "A" "T" "C" "C" "C"
## M099 "T" "A" "T" "C" "C" "C"
## S3_4635222 S3_4635228 S3_4636410 S3_4636546 S3_4641790 S3_4641953
## M186 "A" "T" "T" "C" "G" "T"
## M033 "A" "T" "T" "C" "G" "T"
## M062 "A" "T" "T" "T" "G" "T"
## M085 "A" "T" "T" "C" "G" "T"
## M102 "A" "T" "T" "C" "G" "T"
## M099 "A" "T" "T" "C" "G" "T"
## S3_4642022 S3_4642331 S3_4642791 S3_4643549 S3_4644699 S3_4645223
## M186 "T" "G" "C" "C" "A" "A"
## M033 "T" "G" "C" "C" "A" "A"
## M062 "T" "G" "C" "C" "A" "A"
## M085 "T" "G" "C" "C" "A" "A"
## M102 "T" "G" "C" "C" "A" "A"
## M099 "T" "G" "C" "C" "A" "A"
## S3_4645304 S3_4645905 S3_4646492 S3_4647642 S3_4647969 S3_4648341
## M186 "T" "G" "T" "G" "A" "A"
## M033 "T" "G" "T" "G" "A" "A"
## M062 "T" "G" "T" "G" "A" "A"
## M085 "T" "G" "T" "G" "A" "A"
## M102 "T" "G" "T" "G" "A" "A"
## M099 "T" "G" "T" "G" "A" "A"
## S3_4650075 S3_4652463 S3_4652773 S3_4653026 S3_4653088 S3_4653226
## M186 "G" "C" "A" "C" "C" "C"
## M033 "G" "C" "A" "C" "C" "C"
## M062 "G" "C" "A" "C" "C" "C"
## M085 "G" "C" "A" "C" "C" "C"
## M102 "G" "C" "A" "C" "C" "C"
## M099 "G" "C" "A" "C" "C" "C"
## S3_4653616 S3_4655844 S3_4655917 S3_4660084 S3_4664799 S3_4665277
## M186 "T" "G" "C" "G" "A" "A"
## M033 "T" "G" "C" "G" "A" "A"
## M062 "T" "G" "C" "G" "A" "A"
## M085 "T" "G" "C" "G" "A" "A"
## M102 "T" "G" "C" "G" "A" "A"
## M099 "T" "G" "C" "G" "A" "A"
## S3_4665405 S3_4665434 S3_4666942 S3_4667128 S3_4667921 S3_4669335
## M186 "T" "G" "T" "T" "C" "A"
## M033 "T" "G" "T" "T" "C" "A"
## M062 "T" "G" "T" "T" "C" "A"
## M085 "T" "G" "T" "T" "C" "A"
## M102 "T" "G" "T" "T" "C" "A"
## M099 "T" "G" "T" "T" "C" "A"
## S3_4670658 S3_4671361 S3_4671825 S3_4673462 S3_4674463 S3_4674531
## M186 "T" "G" "C" "A" "C" "A"
## M033 "T" "G" "C" "A" "C" "A"
## M062 "T" "G" "C" "A" "C" "A"
## M085 "T" "G" "C" "A" "C" "A"
## M102 "T" "G" "C" "A" "C" "A"
## M099 "T" "G" "C" "A" "C" "A"
## S3_4674974 S3_4675551 S3_4675644 S3_4675881 S3_4676898 S3_4678272
## M186 "A" "G" "G" "C" "C" "G"
## M033 "A" "G" "G" "C" "C" "G"
## M062 "A" "G" "G" "C" "C" "G"
## M085 "A" "G" "G" "C" "C" "G"
## M102 "A" "G" "G" "C" "C" "G"
## M099 "A" "G" "G" "C" "C" "G"
## S3_4678459 S3_4686304 S3_4687983 S3_4690853
## M186 "G" "G" "T" "G"
## M033 "G" "G" "T" "G"
## M062 "G" "G" "T" "G"
## M085 "G" "G" "T" "G"
## M102 "G" "G" "T" "G"
## M099 "G" "G" "T" "G"dim(locusT)## [1] 193 143dim(good_locus)## [1] 193 118ok - so the above is for calculating the LD between TWO loci. In order to make the correlation matrix for ALL loci of interest - we would need to make an empty matrix and put the values in there. Preferably with a loop!!!
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
#mymatrix_D
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`
}
}SNPdist <- gsub("S3_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
SNPdist## [1] 4571381 4571758 4574165 4575341 4582900 4582915 4583881 4585421 4585728
## [10] 4586470 4586694 4586731 4588214 4588391 4590830 4590854 4590926 4590930
## [19] 4592328 4592891 4595458 4595932 4596024 4598149 4598569 4599668 4599702
## [28] 4599974 4600723 4602668 4602980 4603950 4604326 4604859 4605023 4605188
## [37] 4605315 4605467 4605524 4605702 4605865 4606726 4607110 4610035 4610365
## [46] 4610928 4611231 4611653 4613096 4617940 4619977 4621370 4621908 4622076
## [55] 4623677 4624465 4625544 4625558 4625816 4629034 4630077 4630332 4631533
## [64] 4632496 4632635 4634801 4635222 4635228 4636410 4636546 4641790 4641953
## [73] 4642022 4642331 4642791 4643549 4644699 4645223 4645304 4645905 4646492
## [82] 4647642 4647969 4648341 4650075 4652463 4652773 4653026 4653088 4653226
## [91] 4653616 4655844 4655917 4660084 4664799 4665277 4665405 4665434 4666942
## [100] 4667128 4667921 4669335 4670658 4671361 4671825 4673462 4674463 4674531
## [109] 4674974 4675551 4675644 4675881 4676898 4678272 4678459 4686304 4687983
## [118] 4690853length(SNPdist)## [1] 118dim(mymatrix_D)## [1] 118 118LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S3_4621370", "S3_4623677", "S3_4631533", "S3_4635222", "S3_4635228", "S3_4636410", "S3_4641790", "S3_4642022", "S3_4642331", "S3_4647642", "S3_4652463", "S3_4653226"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus4_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S3_4621370", "S3_4623677", "S3_4631533", "S3_4635222", "S3_4635228", "S3_4636410", "S3_4641790", "S3_4642022", "S3_4642331", "S3_4647642", "S3_4652463", "S3_4653226"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2Let’s save this before we forget into one figure with all the other zoomed-in plots from other interesting phenotypes:
library(cowplot)
mini_manhattans <- plot_grid(Locus1_tr1, Locus1_tr2, Locus1_tr3, Locus1_tr4,Locus1_tr5, ncol=5)
mini_manhattans
pdf("Mitch_Locus4_miniManhattan.pdf", height = 5, width = 30)
plot(mini_manhattans)
dev.off()## quartz_off_screen
## 2single trait association loci
Few other loci to check - just in case come from single trait associations - that might be interesting nevertheless
load(my_files2[34])
my_files2[34]## [1] "TissueNa_me_SALT_gwas2029.rda"results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)
results2_L7 <- subset(results2, results2$LOD > 7)
results2_L7Looks like we have four loci of interest Chr. 2 Pos 9912134 Chr. 4 Pos 17995110 Chr. 4 Pos 18007088 Chr. 4 Pos 18486112 - 18514133
STI <- subset(results2, results2$Chr == 2)
STI <- subset(STI, STI$Pos > (9912134- 50000))
STI <- subset(STI, STI$Pos < (9912134 + 50000))
STINa_Chr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "TissueNa_me_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr2
pdf("Mitch_Locus_Na_Chr2.1_miniManhattan.pdf")
plot(Na_Chr2)
dev.off()## quartz_off_screen
## 2locus <- subset(PimpGeno2, PimpGeno2$chrom == "2")
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (9912134- 50000))
locus <- subset(locus, locus$pos < (9912134 + 50000))
locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
tail(locus2)locusT <- t(locus2)
head(locusT)## S2_9864039 S2_9864510 S2_9864600 S2_9864680 S2_9865743 S2_9866011
## M186 "G" "G" "T" "T" "G" "T"
## M033 "G" "G" "T" "T" "G" "T"
## M062 "G" "G" "T" "T" "G" "T"
## M085 "G" "G" "T" "T" "G" "T"
## M102 "G" "G" "T" "T" "G" "T"
## M099 "G" "G" "T" "T" "G" "T"
## S2_9866538 S2_9866617 S2_9867765 S2_9867877 S2_9868140 S2_9868905
## M186 "A" "G" "G" "G" "G" "C"
## M033 "A" "G" "G" "G" "G" "C"
## M062 "A" "G" "G" "G" "G" "C"
## M085 "A" "G" "G" "G" "G" "C"
## M102 "A" "G" "G" "G" "G" "C"
## M099 "A" "G" "G" "G" "G" "C"
## S2_9869345 S2_9869515 S2_9869516 S2_9870354 S2_9873014 S2_9873802
## M186 "G" "A" "T" "C" "C" "G"
## M033 "G" "A" "T" "C" "C" "G"
## M062 "G" "A" "T" "C" "C" "G"
## M085 "G" "A" "T" "C" "C" "G"
## M102 "G" "A" "T" "C" "C" "G"
## M099 "G" "A" "T" "C" "C" "G"
## S2_9874502 S2_9875049 S2_9877397 S2_9878194 S2_9881483 S2_9881625
## M186 "G" "C" "C" "C" "C" "T"
## M033 "G" "C" "C" "C" "C" "T"
## M062 "G" "C" "C" "C" "C" "T"
## M085 "G" "C" "C" "C" "C" "T"
## M102 "G" "C" "C" "C" "C" "T"
## M099 "G" "C" "C" "C" "C" "T"
## S2_9881763 S2_9883761 S2_9884714 S2_9884989 S2_9885321 S2_9885370
## M186 "G" "C" "A" "C" "G" "C"
## M033 "G" "C" "A" "C" "G" "C"
## M062 "G" "C" "A" "C" "G" "C"
## M085 "G" "T" "A" "C" "G" "C"
## M102 "G" "C" "A" "C" "G" "C"
## M099 "G" "C" "A" "C" "G" "C"
## S2_9885428 S2_9885670 S2_9886011 S2_9886267 S2_9886302 S2_9910564
## M186 "C" "C" "C" "G" "A" "A"
## M033 "C" "C" "C" "G" "A" "A"
## M062 "C" "C" "C" "G" "A" "A"
## M085 "C" "C" "C" "G" "A" "A"
## M102 "C" "C" "C" "G" "A" "A"
## M099 "C" "C" "C" "G" "A" "A"
## S2_9911113 S2_9911273 S2_9911502 S2_9912134 S2_9912150 S2_9912327
## M186 "T" "C" "G" "G" "G" "T"
## M033 "T" "C" "G" "G" "G" "T"
## M062 "T" "C" "G" "G" "G" "T"
## M085 "T" "C" "G" "G" "G" "T"
## M102 "T" "C" "G" "G" "G" "T"
## M099 "T" "C" "G" "G" "G" "T"
## S2_9912583 S2_9912725 S2_9913207 S2_9913486 S2_9913679 S2_9914307
## M186 "A" "C" "C" "C" "A" "G"
## M033 "A" "C" "C" "C" "A" "G"
## M062 "A" "C" "C" "C" "A" "G"
## M085 "A" "C" "C" "C" "A" "G"
## M102 "A" "C" "C" "C" "A" "G"
## M099 "A" "C" "C" "C" "A" "G"
## S2_9915860 S2_9916714 S2_9917030 S2_9917034 S2_9917217 S2_9917416
## M186 "C" "A" "G" "G" "G" "C"
## M033 "C" "A" "G" "G" "G" "C"
## M062 "C" "A" "G" "G" "G" "C"
## M085 "C" "A" "G" "G" "G" "C"
## M102 "C" "A" "G" "G" "G" "C"
## M099 "C" "A" "G" "G" "G" "C"
## S2_9926278 S2_9927586 S2_9927597 S2_9928443 S2_9929195 S2_9929471
## M186 "G" "C" "T" "T" "G" "C"
## M033 "G" "C" "T" "T" "G" "C"
## M062 "G" "C" "T" "T" "G" "C"
## M085 "G" "C" "T" "T" "G" "C"
## M102 "G" "C" "T" "T" "G" "C"
## M099 "G" "C" "T" "T" "G" "C"
## S2_9931061 S2_9931481 S2_9931701 S2_9932342 S2_9932547 S2_9932910
## M186 "C" "C" "C" "C" "C" "T"
## M033 "C" "C" "C" "C" "C" "T"
## M062 "C" "C" "C" "C" "C" "T"
## M085 "C" "C" "C" "C" "C" "T"
## M102 "C" "C" "C" "C" "C" "T"
## M099 "C" "C" "C" "C" "C" "T"
## S2_9933239 S2_9933906 S2_9934742 S2_9934827 S2_9935103 S2_9936032
## M186 "G" "C" "G" "A" "A" "G"
## M033 "G" "C" "G" "A" "A" "G"
## M062 "G" "C" "G" "A" "A" "G"
## M085 "G" "C" "G" "A" "A" "G"
## M102 "G" "C" "G" "A" "A" "G"
## M099 "G" "C" "G" "A" "A" "G"
## S2_9936138 S2_9937415 S2_9937750 S2_9937827 S2_9937998 S2_9938217
## M186 "C" "C" "T" "G" "G" "A"
## M033 "C" "C" "T" "G" "G" "A"
## M062 "C" "C" "T" "G" "G" "A"
## M085 "C" "C" "G" "G" "G" "A"
## M102 "C" "C" "T" "G" "G" "A"
## M099 "C" "C" "T" "G" "G" "A"
## S2_9939747 S2_9941225 S2_9941965 S2_9942959 S2_9943748 S2_9944395
## M186 "G" "C" "C" "C" "C" "G"
## M033 "G" "C" "C" "C" "C" "G"
## M062 "G" "C" "C" "C" "C" "G"
## M085 "G" "C" "C" "C" "C" "A"
## M102 "G" "C" "C" "C" "C" "G"
## M099 "G" "C" "C" "C" "C" "G"
## S2_9945460 S2_9945758 S2_9946055 S2_9946162 S2_9946200 S2_9946307
## M186 "G" "C" "G" "T" "C" "G"
## M033 "G" "C" "G" "T" "C" "G"
## M062 "G" "C" "G" "T" "C" "G"
## M085 "G" "C" "A" "T" "C" "G"
## M102 "G" "C" "G" "T" "C" "G"
## M099 "G" "C" "G" "T" "C" "G"
## S2_9946392 S2_9946785 S2_9947667 S2_9948354 S2_9948472 S2_9949173
## M186 "G" "C" "T" "G" "T" "A"
## M033 "G" "C" "T" "G" "T" "A"
## M062 "G" "C" "T" "G" "T" "A"
## M085 "G" "C" "T" "G" "T" "A"
## M102 "G" "C" "T" "G" "T" "A"
## M099 "G" "C" "T" "G" "T" "A"
## S2_9949856 S2_9950092 S2_9951368 S2_9953377 S2_9954657 S2_9954672
## M186 "A" "C" "G" "G" "C" "C"
## M033 "A" "C" "G" "G" "C" "C"
## M062 "A" "C" "G" "G" "C" "C"
## M085 "A" "C" "G" "G" "C" "C"
## M102 "A" "C" "G" "G" "C" "C"
## M099 "A" "C" "G" "G" "C" "C"
## S2_9955094 S2_9955104 S2_9955838 S2_9956316 S2_9956503 S2_9956593
## M186 "T" "G" "G" "C" "C" "A"
## M033 "T" "G" "G" "C" "C" "A"
## M062 "T" "G" "G" "C" "C" "A"
## M085 "T" "G" "G" "T" "G" "A"
## M102 "T" "G" "G" "C" "C" "A"
## M099 "T" "G" "G" "C" "C" "A"
## S2_9956648 S2_9956793 S2_9957010 S2_9957169 S2_9959219 S2_9959459
## M186 "T" "A" "G" "A" "T" "C"
## M033 "T" "A" "G" "A" "T" "C"
## M062 "T" "A" "G" "A" "T" "C"
## M085 "C" "A" "G" "A" "T" "C"
## M102 "T" "A" "G" "A" "T" "C"
## M099 "T" "A" "G" "A" "T" "C"
## S2_9959887 S2_9960828 S2_9961591 S2_9961706
## M186 "G" "C" "G" "C"
## M033 "G" "C" "G" "C"
## M062 "G" "C" "G" "C"
## M085 "G" "C" "G" "A"
## M102 "G" "C" "G" "C"
## M099 "G" "C" "G" "C"dim(locusT)## [1] 193 118dim(locus2)## [1] 118 193unique(locusT[,1])## [1] "G"if(length(unique(locusT[,2])) > 1){
good_locus <- locusT[,2]
}
good_locus## M186 M033 M062 M085 M102 M099 M056 M061 M282 M105 M106 M107 M108 M109 M112 M147
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M158 M159 M135 M273 M275 M199 M200 M059 M010 M013 M053 M065 M022 M072 M078 M070
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M036 M073 M097 M258 M261 M268 M269 M270 M005 M015 M017 M019 M020 M083 M028 M003
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "T" "G" "T" "G" "G" "G" "G"
## M021 M023 M024 M029 M030 M076 M075 M074 M093 M002 M088 M086 M103 M087 M037 M040
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M042 M001 M044 M045 M046 M048 M050 M051 M185 M160 M052 M245 M244 M241 M240 M239
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M238 M237 M236 M166 M140 M141 M142 M143 M144 M145 M006 M007 M008 M167 M009 M113
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M115 M116 M117 M118 M119 M120 M170 M171 M172 M174 M151 M155 M156 M184 M157 M067
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M058 M066 M057 M125 M127 M183 M133 M080 M094 M247 M248 M250 M251 M252 M253 M178
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M255 M256 M271 M284 M095 M096 M198 M176 M187 M189 M190 M191 M192 M196 M197 M221
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M223 M227 M229 M231 M234 M285 M286 M287 M288 M289 M290 M031 M291 M292 M293 M294
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M298 M301 M210 M211 M212 M217 M220 M214 M216 M011 M012 M014 M018 M034 M038 M054
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M055 M060 M064 M079 M082 M136 M138 M169 M179 M181 M193 M201 M207 M209 M225 M228
## "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G" "G"
## M264
## "G"ncol(locusT)## [1] 118counter <- 2
for(i in 3:118){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1
}
}
head(good_locus)## good_locus S2_9864600 S2_9864680 S2_9865743 S2_9866011 S2_9866538
## M186 "G" "T" "T" "G" "T" "A"
## M033 "G" "T" "T" "G" "T" "A"
## M062 "G" "T" "T" "G" "T" "A"
## M085 "G" "T" "T" "G" "T" "A"
## M102 "G" "T" "T" "G" "T" "A"
## M099 "G" "T" "T" "G" "T" "A"
## S2_9866617 S2_9867765 S2_9867877 S2_9868140 S2_9868905 S2_9869345
## M186 "G" "G" "G" "G" "C" "G"
## M033 "G" "G" "G" "G" "C" "G"
## M062 "G" "G" "G" "G" "C" "G"
## M085 "G" "G" "G" "G" "C" "G"
## M102 "G" "G" "G" "G" "C" "G"
## M099 "G" "G" "G" "G" "C" "G"
## S2_9869515 S2_9869516 S2_9870354 S2_9873014 S2_9873802 S2_9874502
## M186 "A" "T" "C" "C" "G" "G"
## M033 "A" "T" "C" "C" "G" "G"
## M062 "A" "T" "C" "C" "G" "G"
## M085 "A" "T" "C" "C" "G" "G"
## M102 "A" "T" "C" "C" "G" "G"
## M099 "A" "T" "C" "C" "G" "G"
## S2_9875049 S2_9877397 S2_9878194 S2_9881483 S2_9881625 S2_9881763
## M186 "C" "C" "C" "C" "T" "G"
## M033 "C" "C" "C" "C" "T" "G"
## M062 "C" "C" "C" "C" "T" "G"
## M085 "C" "C" "C" "C" "T" "G"
## M102 "C" "C" "C" "C" "T" "G"
## M099 "C" "C" "C" "C" "T" "G"
## S2_9883761 S2_9884714 S2_9884989 S2_9885370 S2_9885428 S2_9885670
## M186 "C" "A" "C" "C" "C" "C"
## M033 "C" "A" "C" "C" "C" "C"
## M062 "C" "A" "C" "C" "C" "C"
## M085 "T" "A" "C" "C" "C" "C"
## M102 "C" "A" "C" "C" "C" "C"
## M099 "C" "A" "C" "C" "C" "C"
## S2_9886011 S2_9886267 S2_9886302 S2_9910564 S2_9911113 S2_9911273
## M186 "C" "G" "A" "A" "T" "C"
## M033 "C" "G" "A" "A" "T" "C"
## M062 "C" "G" "A" "A" "T" "C"
## M085 "C" "G" "A" "A" "T" "C"
## M102 "C" "G" "A" "A" "T" "C"
## M099 "C" "G" "A" "A" "T" "C"
## S2_9911502 S2_9912134 S2_9912150 S2_9912327 S2_9912583 S2_9912725
## M186 "G" "G" "G" "T" "A" "C"
## M033 "G" "G" "G" "T" "A" "C"
## M062 "G" "G" "G" "T" "A" "C"
## M085 "G" "G" "G" "T" "A" "C"
## M102 "G" "G" "G" "T" "A" "C"
## M099 "G" "G" "G" "T" "A" "C"
## S2_9913207 S2_9913486 S2_9914307 S2_9916714 S2_9917030 S2_9917034
## M186 "C" "C" "G" "A" "G" "G"
## M033 "C" "C" "G" "A" "G" "G"
## M062 "C" "C" "G" "A" "G" "G"
## M085 "C" "C" "G" "A" "G" "G"
## M102 "C" "C" "G" "A" "G" "G"
## M099 "C" "C" "G" "A" "G" "G"
## S2_9917217 S2_9917416 S2_9926278 S2_9927586 S2_9927597 S2_9928443
## M186 "G" "C" "G" "C" "T" "T"
## M033 "G" "C" "G" "C" "T" "T"
## M062 "G" "C" "G" "C" "T" "T"
## M085 "G" "C" "G" "C" "T" "T"
## M102 "G" "C" "G" "C" "T" "T"
## M099 "G" "C" "G" "C" "T" "T"
## S2_9929195 S2_9929471 S2_9931061 S2_9931481 S2_9931701 S2_9932342
## M186 "G" "C" "C" "C" "C" "C"
## M033 "G" "C" "C" "C" "C" "C"
## M062 "G" "C" "C" "C" "C" "C"
## M085 "G" "C" "C" "C" "C" "C"
## M102 "G" "C" "C" "C" "C" "C"
## M099 "G" "C" "C" "C" "C" "C"
## S2_9932547 S2_9932910 S2_9933239 S2_9934742 S2_9934827 S2_9935103
## M186 "C" "T" "G" "G" "A" "A"
## M033 "C" "T" "G" "G" "A" "A"
## M062 "C" "T" "G" "G" "A" "A"
## M085 "C" "T" "G" "G" "A" "A"
## M102 "C" "T" "G" "G" "A" "A"
## M099 "C" "T" "G" "G" "A" "A"
## S2_9936032 S2_9936138 S2_9937415 S2_9937750 S2_9937827 S2_9937998
## M186 "G" "C" "C" "T" "G" "G"
## M033 "G" "C" "C" "T" "G" "G"
## M062 "G" "C" "C" "T" "G" "G"
## M085 "G" "C" "C" "G" "G" "G"
## M102 "G" "C" "C" "T" "G" "G"
## M099 "G" "C" "C" "T" "G" "G"
## S2_9938217 S2_9939747 S2_9941225 S2_9941965 S2_9942959 S2_9943748
## M186 "A" "G" "C" "C" "C" "C"
## M033 "A" "G" "C" "C" "C" "C"
## M062 "A" "G" "C" "C" "C" "C"
## M085 "A" "G" "C" "C" "C" "C"
## M102 "A" "G" "C" "C" "C" "C"
## M099 "A" "G" "C" "C" "C" "C"
## S2_9944395 S2_9945460 S2_9945758 S2_9946055 S2_9946200 S2_9946307
## M186 "G" "G" "C" "G" "C" "G"
## M033 "G" "G" "C" "G" "C" "G"
## M062 "G" "G" "C" "G" "C" "G"
## M085 "A" "G" "C" "A" "C" "G"
## M102 "G" "G" "C" "G" "C" "G"
## M099 "G" "G" "C" "G" "C" "G"
## S2_9946392 S2_9947667 S2_9948354 S2_9949173 S2_9949856 S2_9950092
## M186 "G" "T" "G" "A" "A" "C"
## M033 "G" "T" "G" "A" "A" "C"
## M062 "G" "T" "G" "A" "A" "C"
## M085 "G" "T" "G" "A" "A" "C"
## M102 "G" "T" "G" "A" "A" "C"
## M099 "G" "T" "G" "A" "A" "C"
## S2_9953377 S2_9954672 S2_9955094 S2_9955104 S2_9955838 S2_9956316
## M186 "G" "C" "T" "G" "G" "C"
## M033 "G" "C" "T" "G" "G" "C"
## M062 "G" "C" "T" "G" "G" "C"
## M085 "G" "C" "T" "G" "G" "T"
## M102 "G" "C" "T" "G" "G" "C"
## M099 "G" "C" "T" "G" "G" "C"
## S2_9956503 S2_9956648 S2_9956793 S2_9957010 S2_9957169 S2_9959219
## M186 "C" "T" "A" "G" "A" "T"
## M033 "C" "T" "A" "G" "A" "T"
## M062 "C" "T" "A" "G" "A" "T"
## M085 "G" "C" "A" "G" "A" "T"
## M102 "C" "T" "A" "G" "A" "T"
## M099 "C" "T" "A" "G" "A" "T"
## S2_9959459 S2_9959887 S2_9960828 S2_9961591 S2_9961706
## M186 "C" "G" "C" "G" "C"
## M033 "C" "G" "C" "G" "C"
## M062 "C" "G" "C" "G" "C"
## M085 "C" "G" "C" "G" "A"
## M102 "C" "G" "C" "G" "C"
## M099 "C" "G" "C" "G" "C"colnames(good_locus)[1] <- colnames(locusT)[1]
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`
}
}
colnames(mymatrix_D)## [1] "S2_9864039" "S2_9864600" "S2_9864680" "S2_9865743" "S2_9866011"
## [6] "S2_9866538" "S2_9866617" "S2_9867765" "S2_9867877" "S2_9868140"
## [11] "S2_9868905" "S2_9869345" "S2_9869515" "S2_9869516" "S2_9870354"
## [16] "S2_9873014" "S2_9873802" "S2_9874502" "S2_9875049" "S2_9877397"
## [21] "S2_9878194" "S2_9881483" "S2_9881625" "S2_9881763" "S2_9883761"
## [26] "S2_9884714" "S2_9884989" "S2_9885370" "S2_9885428" "S2_9885670"
## [31] "S2_9886011" "S2_9886267" "S2_9886302" "S2_9910564" "S2_9911113"
## [36] "S2_9911273" "S2_9911502" "S2_9912134" "S2_9912150" "S2_9912327"
## [41] "S2_9912583" "S2_9912725" "S2_9913207" "S2_9913486" "S2_9914307"
## [46] "S2_9916714" "S2_9917030" "S2_9917034" "S2_9917217" "S2_9917416"
## [51] "S2_9926278" "S2_9927586" "S2_9927597" "S2_9928443" "S2_9929195"
## [56] "S2_9929471" "S2_9931061" "S2_9931481" "S2_9931701" "S2_9932342"
## [61] "S2_9932547" "S2_9932910" "S2_9933239" "S2_9934742" "S2_9934827"
## [66] "S2_9935103" "S2_9936032" "S2_9936138" "S2_9937415" "S2_9937750"
## [71] "S2_9937827" "S2_9937998" "S2_9938217" "S2_9939747" "S2_9941225"
## [76] "S2_9941965" "S2_9942959" "S2_9943748" "S2_9944395" "S2_9945460"
## [81] "S2_9945758" "S2_9946055" "S2_9946200" "S2_9946307" "S2_9946392"
## [86] "S2_9947667" "S2_9948354" "S2_9949173" "S2_9949856" "S2_9950092"
## [91] "S2_9953377" "S2_9954672" "S2_9955094" "S2_9955104" "S2_9955838"
## [96] "S2_9956316" "S2_9956503" "S2_9956648" "S2_9956793" "S2_9957010"
## [101] "S2_9957169" "S2_9959219" "S2_9959459" "S2_9959887" "S2_9960828"
## [106] "S2_9961591" "S2_9961706"dim(mymatrix_D)## [1] 107 107dim(good_locus)## [1] 193 107head(good_locus)## S2_9864039 S2_9864600 S2_9864680 S2_9865743 S2_9866011 S2_9866538
## M186 "G" "T" "T" "G" "T" "A"
## M033 "G" "T" "T" "G" "T" "A"
## M062 "G" "T" "T" "G" "T" "A"
## M085 "G" "T" "T" "G" "T" "A"
## M102 "G" "T" "T" "G" "T" "A"
## M099 "G" "T" "T" "G" "T" "A"
## S2_9866617 S2_9867765 S2_9867877 S2_9868140 S2_9868905 S2_9869345
## M186 "G" "G" "G" "G" "C" "G"
## M033 "G" "G" "G" "G" "C" "G"
## M062 "G" "G" "G" "G" "C" "G"
## M085 "G" "G" "G" "G" "C" "G"
## M102 "G" "G" "G" "G" "C" "G"
## M099 "G" "G" "G" "G" "C" "G"
## S2_9869515 S2_9869516 S2_9870354 S2_9873014 S2_9873802 S2_9874502
## M186 "A" "T" "C" "C" "G" "G"
## M033 "A" "T" "C" "C" "G" "G"
## M062 "A" "T" "C" "C" "G" "G"
## M085 "A" "T" "C" "C" "G" "G"
## M102 "A" "T" "C" "C" "G" "G"
## M099 "A" "T" "C" "C" "G" "G"
## S2_9875049 S2_9877397 S2_9878194 S2_9881483 S2_9881625 S2_9881763
## M186 "C" "C" "C" "C" "T" "G"
## M033 "C" "C" "C" "C" "T" "G"
## M062 "C" "C" "C" "C" "T" "G"
## M085 "C" "C" "C" "C" "T" "G"
## M102 "C" "C" "C" "C" "T" "G"
## M099 "C" "C" "C" "C" "T" "G"
## S2_9883761 S2_9884714 S2_9884989 S2_9885370 S2_9885428 S2_9885670
## M186 "C" "A" "C" "C" "C" "C"
## M033 "C" "A" "C" "C" "C" "C"
## M062 "C" "A" "C" "C" "C" "C"
## M085 "T" "A" "C" "C" "C" "C"
## M102 "C" "A" "C" "C" "C" "C"
## M099 "C" "A" "C" "C" "C" "C"
## S2_9886011 S2_9886267 S2_9886302 S2_9910564 S2_9911113 S2_9911273
## M186 "C" "G" "A" "A" "T" "C"
## M033 "C" "G" "A" "A" "T" "C"
## M062 "C" "G" "A" "A" "T" "C"
## M085 "C" "G" "A" "A" "T" "C"
## M102 "C" "G" "A" "A" "T" "C"
## M099 "C" "G" "A" "A" "T" "C"
## S2_9911502 S2_9912134 S2_9912150 S2_9912327 S2_9912583 S2_9912725
## M186 "G" "G" "G" "T" "A" "C"
## M033 "G" "G" "G" "T" "A" "C"
## M062 "G" "G" "G" "T" "A" "C"
## M085 "G" "G" "G" "T" "A" "C"
## M102 "G" "G" "G" "T" "A" "C"
## M099 "G" "G" "G" "T" "A" "C"
## S2_9913207 S2_9913486 S2_9914307 S2_9916714 S2_9917030 S2_9917034
## M186 "C" "C" "G" "A" "G" "G"
## M033 "C" "C" "G" "A" "G" "G"
## M062 "C" "C" "G" "A" "G" "G"
## M085 "C" "C" "G" "A" "G" "G"
## M102 "C" "C" "G" "A" "G" "G"
## M099 "C" "C" "G" "A" "G" "G"
## S2_9917217 S2_9917416 S2_9926278 S2_9927586 S2_9927597 S2_9928443
## M186 "G" "C" "G" "C" "T" "T"
## M033 "G" "C" "G" "C" "T" "T"
## M062 "G" "C" "G" "C" "T" "T"
## M085 "G" "C" "G" "C" "T" "T"
## M102 "G" "C" "G" "C" "T" "T"
## M099 "G" "C" "G" "C" "T" "T"
## S2_9929195 S2_9929471 S2_9931061 S2_9931481 S2_9931701 S2_9932342
## M186 "G" "C" "C" "C" "C" "C"
## M033 "G" "C" "C" "C" "C" "C"
## M062 "G" "C" "C" "C" "C" "C"
## M085 "G" "C" "C" "C" "C" "C"
## M102 "G" "C" "C" "C" "C" "C"
## M099 "G" "C" "C" "C" "C" "C"
## S2_9932547 S2_9932910 S2_9933239 S2_9934742 S2_9934827 S2_9935103
## M186 "C" "T" "G" "G" "A" "A"
## M033 "C" "T" "G" "G" "A" "A"
## M062 "C" "T" "G" "G" "A" "A"
## M085 "C" "T" "G" "G" "A" "A"
## M102 "C" "T" "G" "G" "A" "A"
## M099 "C" "T" "G" "G" "A" "A"
## S2_9936032 S2_9936138 S2_9937415 S2_9937750 S2_9937827 S2_9937998
## M186 "G" "C" "C" "T" "G" "G"
## M033 "G" "C" "C" "T" "G" "G"
## M062 "G" "C" "C" "T" "G" "G"
## M085 "G" "C" "C" "G" "G" "G"
## M102 "G" "C" "C" "T" "G" "G"
## M099 "G" "C" "C" "T" "G" "G"
## S2_9938217 S2_9939747 S2_9941225 S2_9941965 S2_9942959 S2_9943748
## M186 "A" "G" "C" "C" "C" "C"
## M033 "A" "G" "C" "C" "C" "C"
## M062 "A" "G" "C" "C" "C" "C"
## M085 "A" "G" "C" "C" "C" "C"
## M102 "A" "G" "C" "C" "C" "C"
## M099 "A" "G" "C" "C" "C" "C"
## S2_9944395 S2_9945460 S2_9945758 S2_9946055 S2_9946200 S2_9946307
## M186 "G" "G" "C" "G" "C" "G"
## M033 "G" "G" "C" "G" "C" "G"
## M062 "G" "G" "C" "G" "C" "G"
## M085 "A" "G" "C" "A" "C" "G"
## M102 "G" "G" "C" "G" "C" "G"
## M099 "G" "G" "C" "G" "C" "G"
## S2_9946392 S2_9947667 S2_9948354 S2_9949173 S2_9949856 S2_9950092
## M186 "G" "T" "G" "A" "A" "C"
## M033 "G" "T" "G" "A" "A" "C"
## M062 "G" "T" "G" "A" "A" "C"
## M085 "G" "T" "G" "A" "A" "C"
## M102 "G" "T" "G" "A" "A" "C"
## M099 "G" "T" "G" "A" "A" "C"
## S2_9953377 S2_9954672 S2_9955094 S2_9955104 S2_9955838 S2_9956316
## M186 "G" "C" "T" "G" "G" "C"
## M033 "G" "C" "T" "G" "G" "C"
## M062 "G" "C" "T" "G" "G" "C"
## M085 "G" "C" "T" "G" "G" "T"
## M102 "G" "C" "T" "G" "G" "C"
## M099 "G" "C" "T" "G" "G" "C"
## S2_9956503 S2_9956648 S2_9956793 S2_9957010 S2_9957169 S2_9959219
## M186 "C" "T" "A" "G" "A" "T"
## M033 "C" "T" "A" "G" "A" "T"
## M062 "C" "T" "A" "G" "A" "T"
## M085 "G" "C" "A" "G" "A" "T"
## M102 "C" "T" "A" "G" "A" "T"
## M099 "C" "T" "A" "G" "A" "T"
## S2_9959459 S2_9959887 S2_9960828 S2_9961591 S2_9961706
## M186 "C" "G" "C" "G" "C"
## M033 "C" "G" "C" "G" "C"
## M062 "C" "G" "C" "G" "C"
## M085 "C" "G" "C" "G" "A"
## M102 "C" "G" "C" "G" "C"
## M099 "C" "G" "C" "G" "C"SNPdist <- gsub("S2_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
SNPdist## [1] 9864039 9864600 9864680 9865743 9866011 9866538 9866617 9867765 9867877
## [10] 9868140 9868905 9869345 9869515 9869516 9870354 9873014 9873802 9874502
## [19] 9875049 9877397 9878194 9881483 9881625 9881763 9883761 9884714 9884989
## [28] 9885370 9885428 9885670 9886011 9886267 9886302 9910564 9911113 9911273
## [37] 9911502 9912134 9912150 9912327 9912583 9912725 9913207 9913486 9914307
## [46] 9916714 9917030 9917034 9917217 9917416 9926278 9927586 9927597 9928443
## [55] 9929195 9929471 9931061 9931481 9931701 9932342 9932547 9932910 9933239
## [64] 9934742 9934827 9935103 9936032 9936138 9937415 9937750 9937827 9937998
## [73] 9938217 9939747 9941225 9941965 9942959 9943748 9944395 9945460 9945758
## [82] 9946055 9946200 9946307 9946392 9947667 9948354 9949173 9949856 9950092
## [91] 9953377 9954672 9955094 9955104 9955838 9956316 9956503 9956648 9956793
## [100] 9957010 9957169 9959219 9959459 9959887 9960828 9961591 9961706LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S2_9912134"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus_Na_Chr2.1_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S2_9912134"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2STI <- subset(results2, results2$Chr == 4)
STI <- subset(STI, STI$Pos > (17995110- 50000))
STI <- subset(STI, STI$Pos < (17995110 + 50000))
STINa_Chr4_1 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "TissueNa_me_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr4_1
pdf("Mitch_Locus_Na_Chr4.1_miniManhattan.pdf")
plot(Na_Chr4_1)
dev.off()## quartz_off_screen
## 2locus <- subset(PimpGeno2, PimpGeno2$chrom == "4")
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (17995110- 50000))
locus <- subset(locus, locus$pos < (17995110 + 50000))
locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
locusT <- t(locus2)
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1}}
colnames(good_locus)[1] <- colnames(locusT)[1]
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`}}
SNPdist <- gsub("S4_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_17995110"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus_Na_Chr4.1_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_17995110"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2STI <- subset(results2, results2$Chr == 4)
STI <- subset(STI, STI$Pos > (18007088- 50000))
STI <- subset(STI, STI$Pos < (18007088 + 50000))
STINa_Chr4_2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "TissueNa_me_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr4_2
pdf("Mitch_Locus_Na_Chr4.2_miniManhattan.pdf")
plot(Na_Chr4_2)
dev.off()## quartz_off_screen
## 2locus <- subset(PimpGeno2, PimpGeno2$chrom == "4")
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (18007088- 50000))
locus <- subset(locus, locus$pos < (18007088 + 50000))
locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
locusT <- t(locus2)
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1}}
colnames(good_locus)[1] <- colnames(locusT)[1]
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`}}
SNPdist <- gsub("S4_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_18007088"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus_Na_Chr4.2_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_18007088"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2STI <- subset(results2, results2$Chr == 4)
STI <- subset(STI, STI$Pos > (18486112- 50000))
STI <- subset(STI, STI$Pos < (18514133 + 50000))
STINa_Chr4_3 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "TissueNa_me_SALT", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr4_3
pdf("Mitch_Locus_Na_Chr4.3_miniManhattan.pdf")
plot(Na_Chr4_3)
dev.off()## quartz_off_screen
## 2locus <- subset(PimpGeno2, PimpGeno2$chrom == "4")
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (18486112- 50000))
locus <- subset(locus, locus$pos < (18486112 + 50000))
locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
locusT <- t(locus2)
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1}}
colnames(good_locus)[1] <- colnames(locusT)[1]
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`}}
SNPdist <- gsub("S4_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_18486112", "S4_18487801", "S4_18488327", "S4_18495923", "S4_18498875", "S4_18509635", "S4_18514133"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus_Na_Chr4.3_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S4_18007088"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2load(my_files2[30])
my_files2[30]## [1] "STI3.RGR_gwas2029.rda"results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)
results2_L7 <- subset(results2, results2$LOD > 7)
results2_L7Looks like we have only ONE locus of interest Chr. 3 Pos 1198094
STI <- subset(results2, results2$Chr == 10)
STI <- subset(STI, STI$Pos > (1198094- 50000))
STI <- subset(STI, STI$Pos < (1198094 + 50000))
STINa_Chr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "STI3.RGR", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr2
pdf("Mitch_Locus_STI3.RGR_Chr10.1_miniManhattan.pdf")
plot(Na_Chr2)
dev.off()## quartz_off_screen
## 2locus <- subset(PimpGeno2, PimpGeno2$chrom == "10")
locus$pos <- as.numeric(as.character(locus$pos))
locus <- subset(locus, locus$pos > (1198094- 50000))
locus <- subset(locus, locus$pos < (1198094 + 50000))
locus2 <- locus[,12:204]
rownames(locus2) <- locus$`rs#`
locusT <- t(locus2)
if(length(unique(locusT[,1])) > 1){
good_locus <- locusT[,1]
}
counter <- 2
for(i in 2:ncol(locusT)){
if(length(unique(locusT[,i])) > 1){
good_locus <- cbind(good_locus, locusT[,i])
colnames(good_locus)[counter] <- colnames(locusT)[i]
counter <- counter + 1}}
colnames(good_locus)[1] <- colnames(locusT)[1]
mymatrix_D <- matrix(, nrow = dim(good_locus)[2], ncol = dim(good_locus)[2])
colnames(mymatrix_D) <- colnames(good_locus)
rownames(mymatrix_D) <- colnames(good_locus)
for(c in 1:ncol(mymatrix_D)){
for(r in 1:nrow(mymatrix_D)){
geno1 <- genotype(good_locus[,c], good_locus[,c])
geno2 <- genotype(good_locus[,r], good_locus[,r])
LD_matrix <- LD(geno1, geno2)
mymatrix_D[r,c]<- LD_matrix$`D'`}}
SNPdist <- gsub("S10_", "", colnames(mymatrix_D))
SNPdist <- as.numeric(as.character(SNPdist))
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S10_1198094"),
color = heat.colors(20), add.key = TRUE)
pdf("Mitch_Locus_STI3.RGR_Chr10.1_AND_RGR.10to14.BLUE_SALT_LD_heatmap.pdf", width = 10, height = 10)
LDheatmap(mymatrix_D, SNPdist, LDmeasure = "D'", title= "Pairwise LD in D'", add.map = TRUE, name = "LD_locus1",
SNP.name = c("S10_1198094"),
color = heat.colors(20), add.key = TRUE)
dev.off()## quartz_off_screen
## 2load(my_files2[11])
my_files2[11]## [1] "Area.smooth.RGR.10to14.BLUE_SALT_gwas2029.rda"results2 <- merge(SNP, results)
results2$LOD <- -log10(results2$Pval)
results2_L7 <- subset(results2, results2$LOD > 7)
results2_L7Looks like we have only ONE locus of interest Chr. 3 Pos 1198094 - and it’s EXACTLY the same one as the previous locus! What are the chances!?!?
STI <- subset(results2, results2$Chr == 10)
STI <- subset(STI, STI$Pos > (1198094- 50000))
STI <- subset(STI, STI$Pos < (1198094 + 50000))
STINa_Chr2 <- plot(ggscatter(STI, x = "Pos", y = "LOD", size = "MAC", color = "darkorchid4", title = "STI3.RGR", ylim=c(0, 12))+theme(legend.position = "right") + labs(size="MAC"))
Na_Chr2
pdf("Mitch_Locus_RGR.10to14.BLUE_SALT_Chr10.1_miniManhattan.pdf")
plot(Na_Chr2)
dev.off()## quartz_off_screen
## 2