On August 25th 2017:
The files from this analysis are saved in this folder: /uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/ExperimentalComparison/march17_finalAnalysis/effectsizes/stdAlleleFreqFiles
-- expComp_stdaf* = files for GLA
-- sla_expComp_stdaf* = files for SLA
-- gla_expComp.txt, sla_expComp.txt = files with combined data from the files above
Rcode for all this is saved in this folder:
/uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/ExperimentalComparison/march17_finalAnalysis/
--
ZG asked me to check the signs of the top snps shared between experimental results and baypass. So to that I first went and used the same code with glm for getting the top 0.01% SNPs (0.99 probs) for all glaMs. I just found the shared SNPs as follows:
#all glams survival top 0.01%
##note that mstop1 and bftop1 are in the code here:
top0.01=which(mstop1 %in% bftop1)
top0.01 #126 137 219 223 290 299 313 624 685
As these are just row numbers for the SNPs shared, I found the details of these rows from the main survival dataframe as follows:
survival[126,] #
scaf pos glaAc glaMs slaAc slaMs all_bf west_bf
126 0 540695 5.497935e-05 0.0002109367 NA NA -8.757504 -8.617567
east_bf all_af west_af east_af gla_af sla_af
126 -3.719619 0.00994 0.01451765 0.0002125 0.0398 1e-04
To check for the bayPass effect result for the reference alleles, I used the posterior mean of the standardized allele frequencies across each population for each locus.(THIS IS VERY CONFUSING!!) So here are the steps to do this, I have to be extra careful to make sure I am cross referencing the correct SNPs.
1). The output file from bayPass which contains the standardized allele frequencies for the analysis is : *_summary_pij.out. The allele frequency are in the column M_Pstd.
2). Since we need the allele frequencies for a particular SNP which is common in the experimental population: we need to find the SNP from the main list which I created for the final analysis with scaffold, position and bayesfactors. Once I find the SNP based on the position in the bfmeansScaffoldPosition_west file, I will look for the row number of this SNP as the locuslist with scaffold, position as the locuslist file and the output from BayPass are in the same order so the row number will correspond to the marker (MRK) number in the bayPass output file. So first look for the row number marker in the file with bayesfactors, then look for that marker in the standard allele freq out file from bayPass.
3). We ran 4 runs, so take averages of the allele frequency posterior means across the 4 runs for east and west.
4). After that, find out which population (1->25) are associated with alfalfa or which are not. (Populations are numbered in the bayPass output file).
5). Then take allele frequency means for that marker for the alfalfa pop and the other pops.
Here is the link to the google sheet for results: https://docs.google.com/spreadsheets/d/1xVghPDVE8-zIvyMbGv6rNav8WSfjuc88v3K85E04xlI/edit#gid=0
Here is the link to the google doc for allele freqs for each population for shared SNPs between experiment and bayPass:
https://docs.google.com/document/d/17fW3kdeARLEt9vVAjfF3RZX35nZ-ntre7-Ce848_pmM/edit
#######################################################################
Here is the code:
##all the correlation between west std allele freqs from bayPass output were 0.9813
###taking averages of the allele frequencies and creating a final file
##west
getwd()
[1] "/uufs/chpc.utah.edu/common/home/gompert-group1/projects/lmelissa_hostAdaptation/ExperimentalComparison/march17_finalAnalysis"
westc0 = read.table("./effectsizes/westc0_summary_pij.out", header=TRUE)
westc1 = read.table("./effectsizes/westc1_summary_pij.out", header=TRUE)
westc2 = read.table("./effectsizes/westc2_summary_pij.out", header=TRUE)
westc3 = read.table("./effectsizes/westc3_summary_pij.out", header=TRUE)
head(westc0)
POP MRK M_P SD_P M_Pstd SD_Pstd DELTA_P ACC_P
1 1 1 0.00056266 0.00609419 -1.0266206 0.6565197 0.0625 0.29738
2 2 1 0.00022675 0.00233418 -0.4678065 0.8712842 0.0625 0.29738
3 3 1 0.00013610 0.00148083 -0.1788254 0.9749945 0.0625 0.29738
4 4 1 0.00060143 0.00443222 -0.2106406 0.9417812 0.0625 0.29738
5 5 1 0.00051683 0.00374661 -0.1987589 0.9800291 0.0625 0.29738
6 6 1 0.00144003 0.00808444 -0.1278349 0.9791656 0.0625 0.29738
#do correlations
cor.test(westc0$M_Pstd, westc1$M_Pstd) #0.98
cor.test(eastc0$M_Pstd, eastc0$M_Pstd) #0.95
#taking averages of all the allele freqs from the 4 runs
avg_af_w = (westc0$M_Pstd + westc1$M_Pstd + westc2$M_Pstd + westc3$M_Pstd) /4
avg_af_e = (eastc0$M_Pstd + eastc1$M_Pstd + eastc2$M_Pstd + eastc3$M_Pstd) /4
head(avg_af)
[1] -0.9961686 -0.4486819 -0.1716657 -0.1474564 -0.1815082 -0.1166823
#create the final data file with population and snp numbers
af_datfile = cbind(westc0$POP, westc0$MRK, avg_af)
af_datfile_e = cbind(eastc0$POP, eastc0$MRK, avg_af_e)
head(af_datfile)
avg_af
[1,] 1 1 -0.9961686
[2,] 2 1 -0.4486819
[3,] 3 1 -0.1716657
[4,] 4 1 -0.1474564
[5,] 5 1 -0.1815082
[6,] 6 1 -0.1166823
colnames(af_datfile_e) <- c("pop","mrk","std_af")
af_datfile_e <- data.frame(af_datfile_e)
head(af_datfile)
pop mrk std_af
[1,] 1 1 -0.9961686
[2,] 2 1 -0.4486819
[3,] 3 1 -0.1716657
[4,] 4 1 -0.1474564
[5,] 5 1 -0.1815082
[6,] 6 1 -0.1166823
bfwest = read.table("bfmeansWithScafPosition_west", header=TRUE)
bfeast = read.table("bfmeansWithScafPosition_east", header=TRUE)
#to get the row number of the SNP from the main file
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 263802),] #115
#find the rows of the populations with this marker number in the allele freq file
snp46 = af_datfile[which(af_datfile[,2] == 115),]
#snp88/206
survival[88,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 417352),] #206
snp206 = af_datfile[which(af_datfile[,2] == 206),]
(snp206$std_af[2] + snp206$std_af[3] + snp206$std_af[5] + snp206$std_af[7])/4
(snp206$std_af[1] + snp206$std_af[4] + snp206$std_af[6] + snp206$std_af[8])/4
#snp104/240
survival[104,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 463964),] #240
snp240 = af_datfile[which(af_datfile[,2] == 240),]
(snp240$std_af[2] + snp240$std_af[3] + snp240$std_af[5] + snp240$std_af[7])/4
(snp240$std_af[1] + snp240$std_af[4] + snp240$std_af[6] + snp240$std_af[8])/4
#snp202/474
survival[202,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 817487),] #474
snp474 = af_datfile[which(af_datfile[,2] == 474),]
(snp474$std_af[2] + snp474$std_af[3] + snp474$std_af[5] + snp474$std_af[7])/4
(snp474$std_af[1] + snp474$std_af[4] + snp474$std_af[6] + snp474$std_af[8])/4
#snp219/519
survival[219,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 876848),] #519
snp519 = af_datfile[which(af_datfile[,2] == 519),]
(snp519$std_af[2] + snp519$std_af[3] + snp519$std_af[5] + snp519$std_af[7])/4
(snp519$std_af[1] + snp519$std_af[4] + snp519$std_af[6] + snp519$std_af[8])/4
#snp249
survival[249,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1053825),] #622
snp622 = af_datfile[which(af_datfile[,2] == 622),]
(snp622$std_af[2] + snp622$std_af[3] + snp622$std_af[5] + snp622$std_af[7])/4
(snp622$std_af[1] + snp622$std_af[4] + snp622$std_af[6] + snp622$std_af[8])/4
#snp260
survival[260,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1058723),] #643
snp643 = af_datfile[which(af_datfile[,2] == 643),]
(snp643$std_af[2] + snp643$std_af[3] + snp643$std_af[5] + snp643$std_af[7])/4
(snp643$std_af[1] + snp643$std_af[4] + snp643$std_af[6] + snp643$std_af[8])/4
#snp274
survival[274,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1076430),] #685
snp685 = af_datfile[which(af_datfile[,2] == 685),]
(snp685$std_af[2] + snp685$std_af[3] + snp685$std_af[5] + snp685$std_af[7])/4
(snp685$std_af[1] + snp685$std_af[4] + snp685$std_af[6] + snp685$std_af[8])/4
#snp300
survival[300,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1192252),] #757
snp757 = af_datfile[which(af_datfile[,2] == 757),]
(snp757$std_af[2] + snp757$std_af[3] + snp757$std_af[5] + snp757$std_af[7])/4
(snp757$std_af[1] + snp757$std_af[4] + snp757$std_af[6] + snp757$std_af[8])/4
#snp313
survival[313,]
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1243005),] #788
snp788 = af_datfile[which(af_datfile[,2] == 788),]
(snp788$std_af[2] + snp788$std_af[3] + snp788$std_af[5] + snp788$std_af[7])/4
(snp788$std_af[1] + snp788$std_af[4] + snp788$std_af[6] + snp788$std_af[8])/4
#snp345
survival[345,]
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 78251),] #854
snp854 = af_datfile[which(af_datfile[,2] == 854),]
(snp854$std_af[2] + snp854$std_af[3] + snp854$std_af[5] + snp854$std_af[7])/4
(snp854$std_af[1] + snp854$std_af[4] + snp854$std_af[6] + snp854$std_af[8])/4
#snp397
survival[397,]
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 270783),] #944
snp944 = af_datfile[which(af_datfile[,2] == 944),]
(snp944$std_af[2] + snp944$std_af[3] + snp944$std_af[5] + snp944$std_af[7])/4
(snp944$std_af[1] + snp944$std_af[4] + snp944$std_af[6] + snp944$std_af[8])/4
#snp464
survival[464,]
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 447931),] #1096
snp1096 = af_datfile[which(af_datfile[,2] == 1096),]
(snp1096$std_af[2] + snp1096$std_af[3] + snp1096$std_af[5] + snp1096$std_af[7])/4
(snp1096$std_af[1] + snp1096$std_af[4] + snp1096$std_af[6] + snp1096$std_af[8])/4
#snp519
survival[519,]
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 663897),] #1254
snp1254 = af_datfile[which(af_datfile[,2] == 1254),]
(snp1254$std_af[2] + snp1254$std_af[3] + snp1254$std_af[5] + snp1254$std_af[7])/4
(snp1254$std_af[1] + snp1254$std_af[4] + snp1254$std_af[6] + snp1254$std_af[8])/4
#snp623
survival[623,]
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 148081),] #1495
snp1495 = af_datfile[which(af_datfile[,2] == 1495),]
(snp1495$std_af[2] + snp1495$std_af[3] + snp1495$std_af[5] + snp1495$std_af[7])/4
(snp1495$std_af[1] + snp1495$std_af[4] + snp1495$std_af[6] + snp1495$std_af[8])/4
#snp624
survival[624,]
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 148082),] #1496
snp1496 = af_datfile[which(af_datfile[,2] == 1496),]
(snp1496$std_af[2] + snp1496$std_af[3] + snp1496$std_af[5] + snp1496$std_af[7])/4
(snp1496$std_af[1] + snp1496$std_af[4] + snp1496$std_af[6] + snp1496$std_af[8])/4
#snp685
survival[685,]
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 340419),] #1652
snp1652 = af_datfile[which(af_datfile[,2] == 1652),]
(snp1652$std_af[2] + snp1652$std_af[3] + snp1652$std_af[5] + snp1652$std_af[7])/4
(snp1652$std_af[1] + snp1652$std_af[4] + snp1652$std_af[6] + snp1652$std_af[8])/4
############# west glaAc survival ##########################
20 156 191 237 322 327 453 512 584 674 734 735 750 751 801
survival[20,] #0 215322
survival[156,] #0 639311
survival[191,] #0 804401
survival[237,] #0 947971
survival[322,] #0 1254761
survival[327,] #1 21202
survival[453,] #1 438855
survival[512,] #1 641713
survival[584,] #1 891331
survival[674,] #2 326859
survival[734,] #2 521619
survival[735,] #2 521627
survival[750,] #2 579506
survival[751,] #2 579510
survival[801,] #2 790051
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 215322),] #55
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 639311),] #355
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 804401),] #453
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 947971),] #586
bfwest[which(bfwest[,1] == 0 & bfwest[,2] == 1254761),] #804
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 21202),] #816
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 438855),] #1073
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 641713),] #1236
bfwest[which(bfwest[,1] == 1 & bfwest[,2] == 891331),] #1398
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 326859),] #1620
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 521619),] #1747
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 521627),] #1748
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 579506),] #1782
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 579510),] #1783
bfwest[which(bfwest[,1] == 2 & bfwest[,2] == 790051),] #1882
snp55 = af_datfile[which(af_datfile[,2] == 55),]
snp55
(snp55$std_af[2] + snp55$std_af[3] + snp55$std_af[5] + snp55$std_af[7])/4
(snp55$std_af[1] + snp55$std_af[4] + snp55$std_af[6] + snp55$std_af[8])/4
snp355 = af_datfile[which(af_datfile[,2] == 355),]
snp355
(snp355$std_af[2] + snp355$std_af[3] + snp355$std_af[5] + snp355$std_af[7])/4
(snp355$std_af[1] + snp355$std_af[4] + snp355$std_af[6] + snp355$std_af[8])/4
snp453 = af_datfile[which(af_datfile[,2] == 453),]
snp453
(snp453$std_af[2] + snp453$std_af[3] + snp453$std_af[5] + snp453$std_af[7])/4
(snp453$std_af[1] + snp453$std_af[4] + snp453$std_af[6] + snp453$std_af[8])/4
snp586 = af_datfile[which(af_datfile[,2] == 586),]
snp586
(snp586$std_af[2] + snp586$std_af[3] + snp586$std_af[5] + snp586$std_af[7])/4
(snp586$std_af[1] + snp586$std_af[4] + snp586$std_af[6] + snp586$std_af[8])/4
snp804 = af_datfile[which(af_datfile[,2] == 804),]
snp804
(snp804$std_af[2] + snp804$std_af[3] + snp804$std_af[5] + snp804$std_af[7])/4
(snp804$std_af[1] + snp804$std_af[4] + snp804$std_af[6] + snp804$std_af[8])/4
snp816 = af_datfile[which(af_datfile[,2] == 816),]
snp816
(snp816$std_af[2] + snp816$std_af[3] + snp816$std_af[5] + snp816$std_af[7])/4
(snp816$std_af[1] + snp816$std_af[4] + snp816$std_af[6] + snp816$std_af[8])/4
snp1073 = af_datfile[which(af_datfile[,2] == 1073),]
snp1073
(snp1073$std_af[2] + snp1073$std_af[3] + snp1073$std_af[5] + snp1073$std_af[7])/4
(snp1073$std_af[1] + snp1073$std_af[4] + snp1073$std_af[6] + snp1073$std_af[8])/4
snp1236 = af_datfile[which(af_datfile[,2] == 1236),]
snp1236
(snp1236$std_af[2] + snp1236$std_af[3] + snp1236$std_af[5] + snp1236$std_af[7])/4
(snp1236$std_af[1] + snp1236$std_af[4] + snp1236$std_af[6] + snp1236$std_af[8])/4
snp1398 = af_datfile[which(af_datfile[,2] == 1398),]
snp1398
(snp1398$std_af[2] + snp1398$std_af[3] + snp1398$std_af[5] + snp1398$std_af[7])/4
(snp1398$std_af[1] + snp1398$std_af[4] + snp1398$std_af[6] + snp1398$std_af[8])/4
snp1620 = af_datfile[which(af_datfile[,2] == 1620),]
snp1620
(snp1620$std_af[2] + snp1620$std_af[3] + snp1620$std_af[5] + snp1620$std_af[7])/4
(snp1620$std_af[1] + snp1620$std_af[4] + snp1620$std_af[6] + snp1620$std_af[8])/4
snp1747 = af_datfile[which(af_datfile[,2] == 1747),]
snp1747
(snp1747$std_af[2] + snp1747$std_af[3] + snp1747$std_af[5] + snp1747$std_af[7])/4
(snp1747$std_af[1] + snp1747$std_af[4] + snp1747$std_af[6] + snp1747$std_af[8])/4
snp1748 = af_datfile[which(af_datfile[,2] == 1748),]
snp1748
(snp1748$std_af[2] + snp1748$std_af[3] + snp1748$std_af[5] + snp1748$std_af[7])/4
(snp1748$std_af[1] + snp1748$std_af[4] + snp1748$std_af[6] + snp1748$std_af[8])/4
snp1782 = af_datfile[which(af_datfile[,2] == 1782),]
snp1782
(snp1782$std_af[2] + snp1782$std_af[3] + snp1782$std_af[5] + snp1782$std_af[7])/4
(snp1782$std_af[1] + snp1782$std_af[4] + snp1782$std_af[6] + snp1782$std_af[8])/4
snp1783 = af_datfile[which(af_datfile[,2] == 1783),]
snp1783
(snp1783$std_af[2] + snp1783$std_af[3] + snp1783$std_af[5] + snp1783$std_af[7])/4
(snp1783$std_af[1] + snp1783$std_af[4] + snp1783$std_af[6] + snp1783$std_af[8])/4
snp1882 = af_datfile[which(af_datfile[,2] == 1882),]
snp1882
(snp1882$std_af[2] + snp1882$std_af[3] + snp1882$std_af[5] + snp1882$std_af[7])/4
(snp1882$std_af[1] + snp1882$std_af[4] + snp1882$std_af[6] + snp1882$std_af[8])/4
#########east glaMs survival ############################
9 70 80 115 137 146 147 199 201 313 402 443 567 568 684 719 722
#snp9/25
survival[9,]
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 129970),] #25
snp25 = af_datfile_e[which(af_datfile_e[,2] == 25),]
(snp25$std_af[2] + snp25$std_af[3] + snp25$std_af[4] + snp25$std_af[5]+ snp25$std_af[7]+ snp25$std_af[8]+ snp25$std_af[10]+ snp25$std_af[11]+ snp25$std_af[14]+ snp25$std_af[15] + snp25$std_af[16])/11
(snp25$std_af[1] + snp25$std_af[6] + snp25$std_af[9]+ snp25$std_af[12]+ snp25$std_af[13]+ snp25$std_af[17]) / 6
#snp70
survival[70,] #0 332623
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 332623),] #169
snp169 = af_datfile_e[which(af_datfile_e[,2] == 169),]
(snp169$std_af[2] + snp169$std_af[3] + snp169$std_af[4] + snp169$std_af[5]+ snp169$std_af[7]+ snp169$std_af[8]+ snp169$std_af[10]+ snp169$std_af[11]+ snp169$std_af[14]+ snp169$std_af[15] + snp169$std_af[16])/11
(snp169$std_af[1] + snp169$std_af[6] + snp169$std_af[9]+ snp169$std_af[12]+ snp169$std_af[13]+ snp169$std_af[17]) / 6
#snp80
survival[80,] #0 337573
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 337573),] #192
snp192 = af_datfile_e[which(af_datfile_e[,2] == 192),]
(snp192$std_af[2] + snp192$std_af[3] + snp192$std_af[4] + snp192$std_af[5]+ snp192$std_af[7]+ snp192$std_af[8]+ snp192$std_af[10]+ snp192$std_af[11]+ snp192$std_af[14]+ snp192$std_af[15] + snp192$std_af[16])/11
(snp192$std_af[1] + snp192$std_af[6] + snp192$std_af[9]+ snp192$std_af[12]+ snp192$std_af[13]+ snp192$std_af[17]) / 6
#snp115
survival[115,] #0 506838
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 506838),] #257
snp257 = af_datfile_e[which(af_datfile_e[,2] == 257),]
snp257
(snp257$std_af[2] + snp257$std_af[3] + snp257$std_af[4] + snp257$std_af[5]+ snp257$std_af[7]+ snp257$std_af[8]+ snp257$std_af[10]+ snp257$std_af[11]+ snp257$std_af[14]+ snp257$std_af[15] + snp257$std_af[16])/11
(snp257$std_af[1] + snp257$std_af[6] + snp257$std_af[9]+ snp257$std_af[12]+ snp257$std_af[13]+ snp257$std_af[17]) / 6
#snp137
survival[137,] #0 570407
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 570407),] #311
snp311 = af_datfile_e[which(af_datfile_e[,2] == 311),]
snp311
(snp311$std_af[2] + snp311$std_af[3] + snp311$std_af[4] + snp311$std_af[5]+ snp311$std_af[7]+ snp311$std_af[8]+ snp311$std_af[10]+ snp311$std_af[11]+ snp311$std_af[14]+ snp311$std_af[15] + snp311$std_af[16])/11
(snp311$std_af[1] + snp311$std_af[6] + snp311$std_af[9]+ snp311$std_af[12]+ snp311$std_af[13]+ snp311$std_af[17]) / 6
#snp146
survival[146,] #0 584433
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 584433),] #323
snp323 = af_datfile_e[which(af_datfile_e[,2] == 323),]
snp323
(snp323$std_af[2] + snp323$std_af[3] + snp323$std_af[4] + snp323$std_af[5]+ snp323$std_af[7]+ snp323$std_af[8]+ snp323$std_af[10]+ snp323$std_af[11]+ snp323$std_af[14]+ snp323$std_af[15] + snp323$std_af[16])/11
(snp323$std_af[1] + snp323$std_af[6] + snp323$std_af[9]+ snp323$std_af[12]+ snp323$std_af[13]+ snp323$std_af[17]) / 6
#snp147
survival[147,] # 0 584441
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 584441),] #325
snp325 = af_datfile_e[which(af_datfile_e[,2] == 325),]
snp325
(snp325$std_af[2] + snp325$std_af[3] + snp325$std_af[4] + snp325$std_af[5]+ snp325$std_af[7]+ snp325$std_af[8]+ snp325$std_af[10]+ snp325$std_af[11]+ snp325$std_af[14]+ snp325$std_af[15] + snp325$std_af[16])/11
(snp325$std_af[1] + snp325$std_af[6] + snp325$std_af[9]+ snp325$std_af[12]+ snp325$std_af[13]+ snp325$std_af[17]) / 6
#snp199
survival[199,] #0 817432
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 817432),] #466
snp466 = af_datfile_e[which(af_datfile_e[,2] == 466),]
snp466
(snp466$std_af[2] + snp466$std_af[3] + snp466$std_af[4] + snp466$std_af[5]+ snp466$std_af[7]+ snp466$std_af[8]+ snp466$std_af[10]+ snp466$std_af[11]+ snp466$std_af[14]+ snp466$std_af[15] + snp466$std_af[16])/11
(snp466$std_af[1] + snp466$std_af[6] + snp466$std_af[9]+ snp466$std_af[12]+ snp466$std_af[13]+ snp466$std_af[17]) / 6
#snp201
survival[201,] #0 817480
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 817480),] #471
snp471 = af_datfile_e[which(af_datfile_e[,2] == 471),]
snp471
(snp471$std_af[2] + snp471$std_af[3] + snp471$std_af[4] + snp471$std_af[5]+ snp471$std_af[7]+ snp471$std_af[8]+ snp471$std_af[10]+ snp471$std_af[11]+ snp471$std_af[14]+ snp471$std_af[15] + snp471$std_af[16])/11
(snp471$std_af[1] + snp471$std_af[6] + snp471$std_af[9]+ snp471$std_af[12]+ snp471$std_af[13]+ snp471$std_af[17]) / 6
#snp313
survival[313,] #0 1243005
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 1243005),] #788
snp788 = af_datfile_e[which(af_datfile_e[,2] == 788),]
snp788
(snp788$std_af[2] + snp788$std_af[3] + snp788$std_af[4] + snp788$std_af[5]+ snp788$std_af[7]+ snp788$std_af[8]+ snp788$std_af[10]+ snp788$std_af[11]+ snp788$std_af[14]+ snp788$std_af[15] + snp788$std_af[16])/11
(snp788$std_af[1] + snp788$std_af[6] + snp788$std_af[9]+ snp788$std_af[12]+ snp788$std_af[13]+ snp788$std_af[17]) / 6
#snp402
survival[402,] #1 280954
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 280954),] #960
snp960 = af_datfile_e[which(af_datfile_e[,2] == 960),]
snp960
(snp960$std_af[2] + snp960$std_af[3] + snp960$std_af[4] + snp960$std_af[5]+ snp960$std_af[7]+ snp960$std_af[8]+ snp960$std_af[10]+ snp960$std_af[11]+ snp960$std_af[14]+ snp960$std_af[15] + snp960$std_af[16])/11
(snp960$std_af[1] + snp960$std_af[6] + snp960$std_af[9]+ snp960$std_af[12]+ snp960$std_af[13]+ snp960$std_af[17]) / 6
#snp443
survival[443,] #1 438516
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 438516),] #1052
snp1052 = af_datfile_e[which(af_datfile_e[,2] == 1052),]
snp1052
(snp1052$std_af[2] + snp1052$std_af[3] + snp1052$std_af[4] + snp1052$std_af[5]+ snp1052$std_af[7]+ snp1052$std_af[8]+ snp1052$std_af[10]+ snp1052$std_af[11]+ snp1052$std_af[14]+ snp1052$std_af[15] + snp1052$std_af[16])/11
(snp1052$std_af[1] + snp1052$std_af[6] + snp1052$std_af[9]+ snp1052$std_af[12]+ snp1052$std_af[13]+ snp1052$std_af[17]) / 6
#snp567
survival[567,] #1 799067
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 799067),] #1353
snp1353 = af_datfile_e[which(af_datfile_e[,2] == 1353),]
snp1353
(snp1353$std_af[2] + snp1353$std_af[3] + snp1353$std_af[4] + snp1353$std_af[5]+ snp1353$std_af[7]+ snp1353$std_af[8]+ snp1353$std_af[10]+ snp1353$std_af[11]+ snp1353$std_af[14]+ snp1353$std_af[15] + snp1353$std_af[16])/11
(snp1353$std_af[1] + snp1353$std_af[6] + snp1353$std_af[9]+ snp1353$std_af[12]+ snp1353$std_af[13]+ snp1353$std_af[17]) / 6
#snp568
survival[568,] #1 799071
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 799071),] #1354
snp1354 = af_datfile_e[which(af_datfile_e[,2] == 1354),]
snp1354
(snp1354$std_af[2] + snp1354$std_af[3] + snp1354$std_af[4] + snp1354$std_af[5]+ snp1354$std_af[7]+ snp1354$std_af[8]+ snp1354$std_af[10]+ snp1354$std_af[11]+ snp1354$std_af[14]+ snp1354$std_af[15] + snp1354$std_af[16])/11
(snp1354$std_af[1] + snp1354$std_af[6] + snp1354$std_af[9]+ snp1354$std_af[12]+ snp1354$std_af[13]+ snp1354$std_af[17]) / 6
#snp684
survival[684,] #2 340399
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 340399),] #1649
snp1649 = af_datfile_e[which(af_datfile_e[,2] == 1649),]
snp1649
(snp1649$std_af[2] + snp1649$std_af[3] + snp1649$std_af[4] + snp1649$std_af[5]+ snp1649$std_af[7]+ snp1649$std_af[8]+ snp1649$std_af[10]+ snp1649$std_af[11]+ snp1649$std_af[14]+ snp1649$std_af[15] + snp1649$std_af[16])/11
(snp1649$std_af[1] + snp1649$std_af[6] + snp1649$std_af[9]+ snp1649$std_af[12]+ snp1649$std_af[13]+ snp1649$std_af[17]) / 6
#snp719
survival[719,] #2 478731
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 478731),] #1726
snp1726 = af_datfile_e[which(af_datfile_e[,2] == 1726),]
snp1726
(snp1726$std_af[2] + snp1726$std_af[3] + snp1726$std_af[4] + snp1726$std_af[5]+ snp1726$std_af[7]+ snp1726$std_af[8]+ snp1726$std_af[10]+ snp1726$std_af[11]+ snp1726$std_af[14]+ snp1726$std_af[15] + snp1726$std_af[16])/11
(snp1726$std_af[1] + snp1726$std_af[6] + snp1726$std_af[9]+ snp1726$std_af[12]+ snp1726$std_af[13]+ snp1726$std_af[17]) / 6
#snp722
survival[722,] #2 491577
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 491577),] #1733
snp1733 = af_datfile_e[which(af_datfile_e[,2] == 1733),]
snp1733
(snp1733$std_af[2] + snp1733$std_af[3] + snp1733$std_af[4] + snp1733$std_af[5]+ snp1733$std_af[7]+ snp1733$std_af[8]+ snp1733$std_af[10]+ snp1733$std_af[11]+ snp1733$std_af[14]+ snp1733$std_af[15] + snp1733$std_af[16])/11
(snp1733$std_af[1] + snp1733$std_af[6] + snp1733$std_af[9]+ snp1733$std_af[12]+ snp1733$std_af[13]+ snp1733$std_af[17]) / 6
### east glaAc
5 55 121 156 187 217 374 393 425 522 587 594 678 741 762 772 820
#snp5
survival[5,] #0 68873
survival[55,] #0 292200
survival[121,] #0 529889
survival[156,] #0 639311
survival[187,] #0 789797
survival[217,] #0 865694
survival[374,] #1 150055
survival[393,] #1 248589
survival[425,] #1 357023
survival[522,] #1 663909
survival[587,] #1 893092
survival[594,] #1 954412
survival[678,] # 2 326895
survival[741,] # 2 524370
survival[762,] # 2 612668
survival[772,] # 2 637116
survival[820,] # 2 932707
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 68873),] #8
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 292200),] #143
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 529889),] #264
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 639311),] #355
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 789797),] #447
bfeast[which(bfeast[,1] == 0 & bfeast[,2] == 865694),] #511
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 150055),] #901
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 248589),] #936
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 357023),] #1008
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 663909),] #1258
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 893092),] #1405
bfeast[which(bfeast[,1] == 1 & bfeast[,2] == 954412),] #1414
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 326895),] #1626
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 524370),] #1757
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 612668),] #1806
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 637116),] #1825
bfeast[which(bfeast[,1] == 2 & bfeast[,2] == 932707),] #1924
snp8 = af_datfile_e[which(af_datfile_e[,2] == 8),]
snp8
(snp8$std_af[2] + snp8$std_af[3] + snp8$std_af[4] + snp8$std_af[5]+ snp8$std_af[7]+ snp8$std_af[8]+ snp8$std_af[10]+ snp8$std_af[11]+ snp8$std_af[14]+ snp8$std_af[15] + snp8$std_af[16])/11
(snp8$std_af[1] + snp8$std_af[6] + snp8$std_af[9]+ snp8$std_af[12]+ snp8$std_af[13]+ snp8$std_af[17]) / 6
snp143 = af_datfile_e[which(af_datfile_e[,2] == 143),]
snp143
(snp143$std_af[2] + snp143$std_af[3] + snp143$std_af[4] + snp143$std_af[5]+ snp143$std_af[7]+ snp143$std_af[8]+ snp143$std_af[10]+ snp143$std_af[11]+ snp143$std_af[14]+ snp143$std_af[15] + snp143$std_af[16])/11
(snp143$std_af[1] + snp143$std_af[6] + snp143$std_af[9]+ snp143$std_af[12]+ snp143$std_af[13]+ snp143$std_af[17]) / 6
snp264 = af_datfile_e[which(af_datfile_e[,2] == 264),]
snp264
(snp264$std_af[2] + snp264$std_af[3] + snp264$std_af[4] + snp264$std_af[5]+ snp264$std_af[7]+ snp264$std_af[8]+ snp264$std_af[10]+ snp264$std_af[11]+ snp264$std_af[14]+ snp264$std_af[15] + snp264$std_af[16])/11
(snp264$std_af[1] + snp264$std_af[6] + snp264$std_af[9]+ snp264$std_af[12]+ snp264$std_af[13]+ snp264$std_af[17]) / 6
snp355 = af_datfile_e[which(af_datfile_e[,2] == 355),]
snp355
(snp355$std_af[2] + snp355$std_af[3] + snp355$std_af[4] + snp355$std_af[5]+ snp355$std_af[7]+ snp355$std_af[8]+ snp355$std_af[10]+ snp355$std_af[11]+ snp355$std_af[14]+ snp355$std_af[15] + snp355$std_af[16])/11
(snp355$std_af[1] + snp355$std_af[6] + snp355$std_af[9]+ snp355$std_af[12]+ snp355$std_af[13]+ snp355$std_af[17]) / 6
snp447 = af_datfile_e[which(af_datfile_e[,2] == 447),]
snp447
(snp447$std_af[2] + snp447$std_af[3] + snp447$std_af[4] + snp447$std_af[5]+ snp447$std_af[7]+ snp447$std_af[8]+ snp447$std_af[10]+ snp447$std_af[11]+ snp447$std_af[14]+ snp447$std_af[15] + snp447$std_af[16])/11
(snp447$std_af[1] + snp447$std_af[6] + snp447$std_af[9]+ snp447$std_af[12]+ snp447$std_af[13]+ snp447$std_af[17]) / 6
snp511 = af_datfile_e[which(af_datfile_e[,2] == 511),]
snp511
(snp511$std_af[2] + snp511$std_af[3] + snp511$std_af[4] + snp511$std_af[5]+ snp511$std_af[7]+ snp511$std_af[8]+ snp511$std_af[10]+ snp511$std_af[11]+ snp511$std_af[14]+ snp511$std_af[15] + snp511$std_af[16])/11
(snp511$std_af[1] + snp511$std_af[6] + snp511$std_af[9]+ snp511$std_af[12]+ snp511$std_af[13]+ snp511$std_af[17]) / 6
snp901 = af_datfile_e[which(af_datfile_e[,2] == 901),]
snp901
(snp901$std_af[2] + snp901$std_af[3] + snp901$std_af[4] + snp901$std_af[5]+ snp901$std_af[7]+ snp901$std_af[8]+ snp901$std_af[10]+ snp901$std_af[11]+ snp901$std_af[14]+ snp901$std_af[15] + snp901$std_af[16])/11
(snp901$std_af[1] + snp901$std_af[6] + snp901$std_af[9]+ snp901$std_af[12]+ snp901$std_af[13]+ snp901$std_af[17]) / 6
snp936 = af_datfile_e[which(af_datfile_e[,2] == 936),]
snp936
(snp936$std_af[2] + snp936$std_af[3] + snp936$std_af[4] + snp936$std_af[5]+ snp936$std_af[7]+ snp936$std_af[8]+ snp936$std_af[10]+ snp936$std_af[11]+ snp936$std_af[14]+ snp936$std_af[15] + snp936$std_af[16])/11
(snp936$std_af[1] + snp936$std_af[6] + snp936$std_af[9]+ snp936$std_af[12]+ snp936$std_af[13]+ snp936$std_af[17]) / 6
snp1008 = af_datfile_e[which(af_datfile_e[,2] == 1008),]
snp1008
(snp1008$std_af[2] + snp1008$std_af[3] + snp1008$std_af[4] + snp1008$std_af[5]+ snp1008$std_af[7]+ snp1008$std_af[8]+ snp1008$std_af[10]+ snp1008$std_af[11]+ snp1008$std_af[14]+ snp1008$std_af[15] + snp1008$std_af[16])/11
(snp1008$std_af[1] + snp1008$std_af[6] + snp1008$std_af[9]+ snp1008$std_af[12]+ snp1008$std_af[13]+ snp1008$std_af[17]) / 6
snp1258 = af_datfile_e[which(af_datfile_e[,2] == 1258),]
snp1258
(snp1258$std_af[2] + snp1258$std_af[3] + snp1258$std_af[4] + snp1258$std_af[5]+ snp1258$std_af[7]+ snp1258$std_af[8]+ snp1258$std_af[10]+ snp1258$std_af[11]+ snp1258$std_af[14]+ snp1258$std_af[15] + snp1258$std_af[16])/11
(snp1258$std_af[1] + snp1258$std_af[6] + snp1258$std_af[9]+ snp1258$std_af[12]+ snp1258$std_af[13]+ snp1258$std_af[17]) / 6
snp1405 = af_datfile_e[which(af_datfile_e[,2] == 1405),]
snp1405
(snp1405$std_af[2] + snp1405$std_af[3] + snp1405$std_af[4] + snp1405$std_af[5]+ snp1405$std_af[7]+ snp1405$std_af[8]+ snp1405$std_af[10]+ snp1405$std_af[11]+ snp1405$std_af[14]+ snp1405$std_af[15] + snp1405$std_af[16])/11
(snp1405$std_af[1] + snp1405$std_af[6] + snp1405$std_af[9]+ snp1405$std_af[12]+ snp1405$std_af[13]+ snp1405$std_af[17]) / 6
snp1414 = af_datfile_e[which(af_datfile_e[,2] == 1414),]
snp1414
(snp1414$std_af[2] + snp1414$std_af[3] + snp1414$std_af[4] + snp1414$std_af[5]+ snp1414$std_af[7]+ snp1414$std_af[8]+ snp1414$std_af[10]+ snp1414$std_af[11]+ snp1414$std_af[14]+ snp1414$std_af[15] + snp1414$std_af[16])/11
(snp1414$std_af[1] + snp1414$std_af[6] + snp1414$std_af[9]+ snp1414$std_af[12]+ snp1414$std_af[13]+ snp1414$std_af[17]) / 6
snp1626 = af_datfile_e[which(af_datfile_e[,2] == 1626),]
snp1626
(snp1626$std_af[2] + snp1626$std_af[3] + snp1626$std_af[4] + snp1626$std_af[5]+ snp1626$std_af[7]+ snp1626$std_af[8]+ snp1626$std_af[10]+ snp1626$std_af[11]+ snp1626$std_af[14]+ snp1626$std_af[15] + snp1626$std_af[16])/11
(snp1626$std_af[1] + snp1626$std_af[6] + snp1626$std_af[9]+ snp1626$std_af[12]+ snp1626$std_af[13]+ snp1626$std_af[17]) / 6
snp1806 = af_datfile_e[which(af_datfile_e[,2] == 1806),]
snp1806
(snp1806$std_af[2] + snp1806$std_af[3] + snp1806$std_af[4] + snp1806$std_af[5]+ snp1806$std_af[7]+ snp1806$std_af[8]+ snp1806$std_af[10]+ snp1806$std_af[11]+ snp1806$std_af[14]+ snp1806$std_af[15] + snp1806$std_af[16])/11
(snp1806$std_af[1] + snp1806$std_af[6] + snp1806$std_af[9]+ snp1806$std_af[12]+ snp1806$std_af[13]+ snp1806$std_af[17]) / 6
snp1825 = af_datfile_e[which(af_datfile_e[,2] == 1825),]
snp1825
(snp1825$std_af[2] + snp1825$std_af[3] + snp1825$std_af[4] + snp1825$std_af[5]+ snp1825$std_af[7]+ snp1825$std_af[8]+ snp1825$std_af[10]+ snp1825$std_af[11]+ snp1825$std_af[14]+ snp1825$std_af[15] + snp1825$std_af[16])/11
(snp1825$std_af[1] + snp1825$std_af[6] + snp1825$std_af[9]+ snp1825$std_af[12]+ snp1825$std_af[13]+ snp1825$std_af[17]) / 6
snp1924 = af_datfile_e[which(af_datfile_e[,2] == 1924),]
snp1924
(snp1924$std_af[2] + snp1924$std_af[3] + snp1924$std_af[4] + snp1924$std_af[5]+ snp1924$std_af[7]+ snp1924$std_af[8]+ snp1924$std_af[10]+ snp1924$std_af[11]+ snp1924$std_af[14]+ snp1924$std_af[15] + snp1924$std_af[16])/11
(snp1924$std_af[1] + snp1924$std_af[6] + snp1924$std_af[9]+ snp1924$std_af[12]+ snp1924$std_af[13]+ snp1924$std_af[17]) / 6
Take the results from west and
******* survival *****************************
#east glaMs
##east glaMs
#glm and filtering
fil_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]
fit_glams = glm(fil_glams$glaMs ~ fil_glams$gla_af, family = gaussian)
fit_east = glm(fil_glams$east_bf ~ fil_glams$east_af , family = gaussian)
#doing correlations and randomization ####
ms = fit_glams$residuals
bf = fit_east$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsm = quantile(ms,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
mstop1 = which(ms > quantile(ms,probs=0.99,na.rm=T))
obs = sum(mstop1 %in% bftop1) #17
obs
snps_glams_e = which(mstop1 %in% bftop1)
snps_glams_e
#9 70 80 115 137 146 147 199 201 313 402 443 567 568 684 719 722
snps_glams_dat_e = rbind(survival[9,],survival[70,], survival[80,],survival[115,],survival[137,],survival[146,],survival[147,], survival[199,],survival[201,],survival[313,],survival[402,],survival[443,],survival[567,],survival[568,],survival[684,],survival[719,],survival[722,])
write.table(snps_glams_dat_e, "glaMs_surv_e.txt", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
##east glaAc
#glm and filtering east ac
fil_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]
fit_glaac = glm(fil_glaac$glaAc ~ fil_glaac$gla_af, family = gaussian)
fit_east = glm(fil_glaac$east_bf ~ fil_glaac$east_af , family = gaussian)
#doing correlations and randomization ####
ac = fit_glaac$residuals
bf = fit_east$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsa = quantile(ac,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
actop1 = which(ac > quantile(ac,probs=0.99,na.rm=T))
obs = sum(actop1 %in% bftop1) #17
snps_glaac_e = which(actop1 %in% bftop1)
snps_glaac_e
##west glaMs
#glm and filtering
fil_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]
fit_glams = glm(fil_glams$glaMs ~ fil_glams$gla_af, family = gaussian)
fit_west = glm(fil_glams$west_bf ~ fil_glams$west_af , family = gaussian)
#doing correlations and randomization ####
ms = fit_glams$residuals
bf = fit_west$residuals
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
mstop1 = which(ms > quantile(ms,probs=0.99,na.rm=T))
obs = sum(mstop1 %in% bftop1) #18
snps_gla_ms = which(mstop1 %in% bftop1)
snps_gla_ms
[1] 5 46 88 104 202 219 249 260 274 300 313 345 397 464 519 623 624 685
snps_glams_data=rbind(survival[5,],survival[46,], survival[88,],survival[104,],survival[202,],survival[219,],survival[249,], survival[260,],survival[274,],survival[300,],survival[313,],survival[5,],survival[397,],survival[464,],survival[519,],survival[623,],survival[624,],survival[685,])
write.table(snps_glams_data, "glaMs_surv.txt", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
#west glaAc
#glm and filtering west ac
fil_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]
fit_glaac = glm(fil_glaac$glaAc ~ fil_glaac$gla_af, family = gaussian)
fit_west = glm(fil_glaac$west_bf ~ fil_glaac$west_af , family = gaussian)
#doing correlations and randomization ####
ac = fit_glaac$residuals
bf = fit_west$residuals
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
actop1 = which(ac > quantile(ac,probs=0.99,na.rm=T))
obs = sum(actop1 %in% bftop1) #15
snps_gla_ac = which(actop1 %in% bftop1)
snps_gla_ac