Gwa from GBLUP
AnS 656
Main Path: “C:/Users/jsteibel/OneDrive/Documents/job/ANS 656/hw/”
Packages
Current Package: DT
Current Package: MASS
Current Package: Matrix
Current Package: GGally
Current Package: plotly
Current Package: feather
Current Package: arrow
Current Package: pedigreemm
Current Package: BGLR
Current Package: regress
Current Package: gwaR
Current Package: lme4
Current Package: lmerTest
Current Package: car
Current Package: emmeans
Current Package: effects
Current Package: rsq
Current Package: performance
Current Package: coda
Current Package: pracma
Current Package: qvalue
Current Package: Matrix
Current Package: skimr
Current Package: data.table
Current Package: readxl
Current Package: lubridate
Current Package: broom
Current Package: stringr
Current Package: tidyverse initialize
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 2949478 157.6 4717335 252.0 4717335 252
Vcells 5496503 42.0 10146329 77.5 8387445 64Folder Path: C:/Users/jsteibel/OneDrive/Documents/job/ANS 656/hw/ Binary R Dataset: all_data.RData Data
Notice here I already read in many of the files you need.
##load data
load(data_file)
data.pheno <- data$Pheno
data.map <- data$Map
data.geno <- data$Geno
rm(data)
gc(verbose = FALSE) used (Mb) gc trigger (Mb) max used (Mb)
Ncells 3041251 162.5 4717674 252.0 4717674 252.0
Vcells 51832108 395.5 68496422 522.6 51885007 395.9# print first few rows and cols
data.geno[1:5, 1:5]# A tibble: 5 × 5
ID MARC0044150 ASGA0000014 H3GA0000026 ASGA0000021
<chr> <dbl> <dbl> <dbl> <dbl>
1 1001 1 1 1 1
2 1002 1 1 1 1
3 1003 1 1 1 1
4 1004 1 1 1 1
5 1006 0 1 2 2# M matrix
M <- Matrix::Matrix(sapply(data.geno[,-1], as.numeric))
# dim
dim(M)[1] 1015 45329# check dims on each file
dim(data.map)[1] 45329 3dim(data.geno)[1] 1015 45330dim(data.pheno)[1] 921 10id_geno<-data.geno$ID
id_pheno<-data.pheno$ID
skim(data.pheno)| Name | data.pheno |
| Number of rows | 921 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 1 |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| ID | 0 | 1 | 4 | 4 | 0 | 921 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Sex | 0 | 1 | FALSE | 2 | M: 490, F: 431 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Age | 0 | 1 | 166.17 | 9.61 | 148.00 | 157.00 | 164.00 | 176.00 | 182.00 | ▂▇▃▂▇ |
| PercDuroc | 0 | 1 | 0.50 | 0.06 | 0.31 | 0.47 | 0.50 | 0.53 | 0.84 | ▁▇▃▁▁ |
| BirthWt | 0 | 1 | 1.54 | 0.32 | 0.64 | 1.32 | 1.54 | 1.72 | 2.45 | ▁▅▇▃▁ |
| ADG | 0 | 1 | 881.11 | 102.55 | 536.06 | 812.25 | 880.50 | 950.38 | 1189.30 | ▁▃▇▅▁ |
| CarcassWt | 0 | 1 | 81.86 | 6.81 | 60.77 | 77.55 | 81.63 | 86.17 | 102.04 | ▁▃▇▅▁ |
| CarcassBF | 0 | 1 | 24.15 | 7.32 | 7.62 | 19.05 | 22.86 | 27.94 | 55.88 | ▂▇▃▁▁ |
| CarcassLMA | 0 | 1 | 40.60 | 4.73 | 28.06 | 37.42 | 40.65 | 43.87 | 58.71 | ▂▇▇▂▁ |
| Color | 0 | 1 | 3.16 | 0.82 | 1.00 | 3.00 | 3.00 | 4.00 | 6.00 | ▃▇▃▁▁ |
Create a G matrix using VanRaden method #1.
To perform this operation we need to center and then standardize the marker matrix. We center using the expected allele dosages. We standardize with the square root of the sum of 2pq.
Important checks: rownames of G, diagonals of G, columns of Z
p_i<-colMeans(M)/2
summary(p_i) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.01084 0.28916 0.51034 0.50804 0.73054 0.98916 Z<-sweep(M,MARGIN = 2,STATS = 2*p_i,FUN = "-")
max(abs((colMeans(Z)))) #check if all columns are centered[1] 2.404645e-15Z<-Z/sqrt(sum(2*p_i*(1-p_i)))
G<-tcrossprod(Z) #it takes a LONG time
summary(diag(G)) #Okish? Min. 1st Qu. Median Mean 3rd Qu. Max.
0.6220 0.9050 0.9429 0.9477 0.9844 1.4530 summary(as.vector(G[upper.tri(G,diag = FALSE)])) Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.5023635 -0.0709889 -0.0200363 -0.0009347 0.0506455 0.9278445 rownames(G)NULLdim(G)[1] 1015 1015dim(data.geno)[1] 1015 45330rownames(G)<-colnames(G)<-id_genoNotice: no NA in allele frequency, all MAF > 0.01. We can refine further and use MAF > 0.05 but I will leave that to a homework.
Let’s perform the same GWA as in the HW
Let’s use:
- ADG as the response variable (y)
- Sex as fixed categorical variable (class in R)
- Age as a covariate (fixed)
Important things to check: rownames of pheno should be ID. Rownames of G and geno matrix should also have animal ID
Convenient things: no need to match rows.
rownames(data.pheno)<-data.pheno$ID #notice the warningWarning: Setting row names on a tibble is deprecated.data.pheno<-as.data.frame(data.pheno) #neccesary due to not being tibble aware
class(data.pheno)[1] "data.frame"G<-as.matrix(G) #necesary due to not being sparse matrix aware
gb<-gblup(rsp = "ADG",design = c("~Sex+Age"),G = G,data = data.pheno)
gbgblup analysis of trait: ADG
fixed effects equation:
y ~ Sex + Age
random effects equation:
~G + In
log-likelihood: -249.6263 converged in: 5 iterations
estimated fixed effects:
Estimate StdError test.st p.value
(Intercept) 1988.597715 168.499939 11.8017712 0.000000e+00
SexM -7.929912 19.919779 -0.3980924 6.905621e-01
Age -6.644845 1.013622 -6.5555452 5.543899e-11
estimated variance components:
Estimate StdError prop.var
G 565.473 1447.148 0.1016316
In 4998.477 1460.166 0.8983684summary(gb,ehat=FALSE, sigma=TRUE, fe=TRUE)summary gblup analysis of trait: ADG
fixed effects equation:
y ~ Sex + Age
random effects equation:
~G + In
log-likelihood: -249.6263 converged in: 5 iterations
estimated fixed effects:
Estimate StdError test.st p.value
(Intercept) 1988.597715 168.499939 11.8017712 0.000000e+00
SexM -7.929912 19.919779 -0.3980924 6.905621e-01
Age -6.644845 1.013622 -6.5555452 5.543899e-11
estimated variance components:
Estimate StdError prop.var
G 565.473 1447.148 0.1016316
In 4998.477 1460.166 0.8983684
breeding value predicted for 55 first 5 animals shown:
[,1]
160 -14.158959
301 7.321146
303 -4.427691
311 1.106294
312 -1.810255
314 -1.056199varcomp(gb) Estimate StdError prop.var se
G 565.473 1447.148 0.1016316 0.2546509
In 4998.477 1460.166 0.8983684 0.3736658Now, let’s look into the gblup object
str(gb)List of 12
$ name : chr "ADG"
$ llik : num -250
$ cycle: int 5
$ Q : num [1:55, 1:55] 0.9344 -0.016 -0.0633 -0.016 0.028 ...
$ V : num [1:55, 1:55] 5451.3 25.8 56.8 -69.8 -25 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:55] "160" "301" "303" "311" ...
.. ..$ : chr [1:55] "160" "301" "303" "311" ...
$ Vinv : num [1:55, 1:55] 1.85e-04 -9.50e-07 -1.78e-06 1.89e-06 6.05e-07 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:55] "160" "301" "303" "311" ...
.. ..$ : NULL
$ sigma: Named num [1:2] 565 4998
..- attr(*, "names")= chr [1:2] "G" "In"
$ coefm: num [1:5, 1:2] 1988.6 -7.93 -6.64 565.47 4998.48 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "(Intercept)" "SexM" "Age" "G" ...
.. ..$ : chr [1:2] "Estimate" "StdError"
$ model:List of 7
..$ y : Named num [1:55] 879 794 896 731 811 ...
.. ..- attr(*, "names")= chr [1:55] "160" "301" "303" "311" ...
..$ Sex : Factor w/ 2 levels "F","M": 1 1 1 1 2 1 1 2 1 1 ...
..$ Age : num [1:55] 154 176 155 176 176 156 177 156 177 154 ...
..$ G : num [1:55, 1:55] 0.8008 0.0456 0.1004 -0.1235 -0.0443 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:55] "160" "301" "303" "311" ...
.. .. ..$ : chr [1:55] "160" "301" "303" "311" ...
..$ In : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ formula :Class 'formula' language y ~ Sex + Age
.. .. ..- attr(*, ".Environment")=<environment: 0x00000169bb69ba58>
..$ Vformula:Class 'formula' language ~G + In
.. .. ..- attr(*, ".Environment")=<environment: 0x00000169bb59f380>
$ V.CV : num [1:2, 1:2] 2094236 -1452785 -1452785 2132085
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "G" "In"
.. ..$ : chr [1:2] "G" "In"
$ ehat : num [1:55, 1] -86.102 -25.315 -62.807 -87.685 -0.375 ...
$ pos : logi [1:2] FALSE FALSE
- attr(*, "class")= chr "gblup"question to the class
Something unsual?
Rinse and repeat
rownames(data.pheno)<-data.pheno$ID #notice the warning
gb<-gblup(rsp = "ADG",design = c("~Sex+Age"),G = G,data = data.pheno)
gbgblup analysis of trait: ADG
fixed effects equation:
y ~ Sex + Age
random effects equation:
~G + In
log-likelihood: -4443.274 converged in: 6 iterations
estimated fixed effects:
Estimate StdError test.st p.value
(Intercept) 1760.17886 46.0572734 38.217175 0.000000e+00
SexM 37.50750 5.3626256 6.994242 2.666978e-12
Age -5.41047 0.2726298 -19.845482 0.000000e+00
estimated variance components:
Estimate StdError prop.var
G 2248.512 442.0670 0.3305127
In 4554.590 289.3812 0.6694873summary(gb,ehat=FALSE, sigma=TRUE, fe=TRUE)summary gblup analysis of trait: ADG
fixed effects equation:
y ~ Sex + Age
random effects equation:
~G + In
log-likelihood: -4443.274 converged in: 6 iterations
estimated fixed effects:
Estimate StdError test.st p.value
(Intercept) 1760.17886 46.0572734 38.217175 0.000000e+00
SexM 37.50750 5.3626256 6.994242 2.666978e-12
Age -5.41047 0.2726298 -19.845482 0.000000e+00
estimated variance components:
Estimate StdError prop.var
G 2248.512 442.0670 0.3305127
In 4554.590 289.3812 0.6694873
breeding value predicted for 921 first 5 animals shown:
[,1]
1001 58.719836
1002 7.837963
1003 42.809134
1004 58.956449
1006 -27.850185
1007 33.175948varcomp(gb) Estimate StdError prop.var se
G 2248.512 442.0670 0.3305127 0.05455006
In 4554.590 289.3812 0.6694873 0.06232086str(gb)List of 12
$ name : chr "ADG"
$ llik : num -4443
$ cycle: int 6
$ Q : num [1:921, 1:921] 0.99321 0.00214 -0.00148 0.00214 -0.00241 ...
$ V : num [1:921, 1:921] 6487 691 588 726 160 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:921] "1001" "1002" "1003" "1004" ...
.. ..$ : chr [1:921] "1001" "1002" "1003" "1004" ...
$ Vinv : num [1:921, 1:921] 1.82e-04 -5.86e-06 -6.00e-06 -7.46e-06 -1.17e-06 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:921] "1001" "1002" "1003" "1004" ...
.. ..$ : NULL
$ sigma: Named num [1:2] 2249 4555
..- attr(*, "names")= chr [1:2] "G" "In"
$ coefm: num [1:5, 1:2] 1760.18 37.51 -5.41 2248.51 4554.59 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:5] "(Intercept)" "SexM" "Age" "G" ...
.. ..$ : chr [1:2] "Estimate" "StdError"
$ model:List of 7
..$ y : Named num [1:921] 881 804 1041 887 680 ...
.. ..- attr(*, "names")= chr [1:921] "1001" "1002" "1003" "1004" ...
..$ Sex : Factor w/ 2 levels "F","M": 2 1 2 1 1 2 2 1 2 2 ...
..$ Age : num [1:921] 182 161 161 161 179 161 161 179 161 179 ...
..$ G : num [1:921, 1:921] 0.8595 0.3075 0.2617 0.3228 0.0714 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:921] "1001" "1002" "1003" "1004" ...
.. .. ..$ : chr [1:921] "1001" "1002" "1003" "1004" ...
..$ In : Factor w/ 921 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
..$ formula :Class 'formula' language y ~ Sex + Age
.. .. ..- attr(*, ".Environment")=<environment: 0x00000169c39247a8>
..$ Vformula:Class 'formula' language ~G + In
.. .. ..- attr(*, ".Environment")=<environment: 0x00000169c38d9ce0>
$ V.CV : num [1:2, 1:2] 195423 -62000 -62000 83741
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "G" "In"
.. ..$ : chr [1:2] "G" "In"
$ ehat : num [1:921, 1] 68.43 -84.99 114.6 -2.52 -111.31 ...
$ pos : logi [1:2] FALSE FALSE
- attr(*, "class")= chr "gblup"Not we can perform GWA.
We have two matrices: M and Z are them any different for generating a gwa?
class(Z)[1] "dgeMatrix"
attr(,"package")
[1] "Matrix"Z<-as.matrix(Z)
rownames(Z)<-id_geno
gw1<-gwas(gblup = gb,x = t(Z))
str(gw1)Classes 'gwas' and 'data.frame': 45329 obs. of 2 variables:
$ ghat : num 2.01 1.83 -2.82 -2.01 -2.82 ...
$ var_ghat: num 8.22 7.62 8.16 8.22 8.16 ...head(gw1) ghat var_ghat
MARC0044150 2.012046 8.219301
ASGA0000014 1.826943 7.624711
H3GA0000026 -2.824697 8.159150
ASGA0000021 -2.012046 8.219301
ALGA0000009 -2.824697 8.159150
ALGA0000014 -2.824697 8.159150summary(gw1) <0.001 <0.01 <0.025 <0.05 <0.1
signif.values 0 0 0 0 0
attr(,"class")
[1] "summary.gwas"summary(as.matrix(gw1)) ghat var_ghat
Min. :-13.88191 Min. : 0.1683
1st Qu.: -1.80213 1st Qu.: 5.8041
Median : -0.01486 Median : 8.5051
Mean : -0.02838 Mean : 8.3857
3rd Qu.: 1.77731 3rd Qu.:10.8911
Max. : 13.75541 Max. :34.7099 zvalues<-gw1$ghat/sqrt(gw1$var_ghat)
summary(zvalues) Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.323896 -0.674994 -0.005821 -0.010390 0.657597 4.278006 This is not very promising because the zvalues are too small
#test significance more formally obtain pvalues and qvalues
pv<-getpvalue(gwas = gw1,log.p = F)
summary(pv) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000153 0.2521190 0.5058183 0.5025033 0.7504822 0.9999919 qv<-qvalue(pv)
plot(qv)
summary(qv)
Call:
qvalue(p = pv)
pi0: 1
Cumulative number of significant calls:
<1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
p-value 12 75 584 1269 2372 4490 45329
q-value 0 0 0 0 0 0 45329
local FDR 0 0 0 0 0 0 1450qqgplot(pvector = pv)
head(data.map)# A tibble: 6 × 3
SNP Chr Pos
<chr> <dbl> <dbl>
1 MARC0044150 1 286933
2 ASGA0000014 1 342481
3 H3GA0000026 1 389876
4 ASGA0000021 1 489855
5 ALGA0000009 1 538161
6 ALGA0000014 1 565627map<-data.frame(data.map[,2:3])
rownames(map)<-data.map$SNP
colnames(map)<-c("chr","pos")
manhattan_plot(pvalues = pv,map = map,threshold = 0.000001)
manhattan_plot(pvalues = qv$qvalues,map = map,threshold = 0.1)
Not much at all. A temptation would be to relax the significance threshold and proceed as if there are significant results. Let’s discuss pro and cons of that approach.
Now, let’s also consider other traits to perform GWA on.
Carcass Weigth
LMA
Backfat
gb2<-gblup(rsp = "CarcassBF",design = c("~Sex+Age"),G = G,data = data.pheno)
summary(gb2,ehat=FALSE, sigma=TRUE, fe=TRUE)summary gblup analysis of trait: CarcassBF
fixed effects equation:
y ~ Sex + Age
random effects equation:
~G + In
log-likelihood: -2057.773 converged in: 5 iterations
estimated fixed effects:
Estimate StdError test.st p.value
(Intercept) 19.026668176 3.40114815 5.5941898 2.216544e-08
SexM 6.614743309 0.40138682 16.4797221 0.000000e+00
Age 0.009639276 0.02013143 0.4788174 6.320686e-01
estimated variance components:
Estimate StdError prop.var
G 17.68178 2.930463 0.4352325
In 22.94428 1.563773 0.5647675
breeding value predicted for 921 first 5 animals shown:
[,1]
1001 1.4102736
1002 0.6626867
1003 -0.1805185
1004 -5.1672973
1006 1.7765257
1007 -4.9536736varcomp(gb2) Estimate StdError prop.var se
G 17.68178 2.930463 0.4352325 0.05624899
In 22.94428 1.563773 0.5647675 0.05939821gw2<-gwas(gblup = gb2,x = t(Z))
pv2<-getpvalue(gwas = gw2,log.p = F)
qv2<-qvalue(pv2)
qqgplot(pvector = pv2)
manhattan_plot(pvalues = qv2$qvalues,map = map,threshold = 0.1)
let’s now try to use an alternative matrix of genotypes
M<-as.matrix(M)
rownames(M)<-rownames(Z)
gw3<-gwas(gblup = gb2,x = t(M))
head(gw2) ghat var_ghat
MARC0044150 0.01397086 0.08391781
ASGA0000014 0.06736581 0.08012698
H3GA0000026 -0.04771612 0.08319594
ASGA0000021 -0.01397086 0.08391781
ALGA0000009 -0.04771612 0.08319594
ALGA0000014 -0.04771612 0.08319594head(gw3) ghat var_ghat
MARC0044150 1.777401 1358.247
ASGA0000014 8.570415 1296.891
H3GA0000026 -6.070542 1346.563
ASGA0000021 -1.777401 1358.247
ALGA0000009 -6.070542 1346.563
ALGA0000014 -6.070542 1346.563pv3<-getpvalue(gwas = gw3,log.p = F)
head(pv3)MARC0044150 ASGA0000014 H3GA0000026 ASGA0000021 ALGA0000009 ALGA0000014
0.9615348 0.8118925 0.8686056 0.9615348 0.8686056 0.8686056 head(pv2)MARC0044150 ASGA0000014 H3GA0000026 ASGA0000021 ALGA0000009 ALGA0000014
0.9615348 0.8118925 0.8686056 0.9615348 0.8686056 0.8686056 Let’s now build an IC against the highest peak
find the peak in p-value and look at the corresponding region
for instance: +/- 10 MB
which.min(pv2)ALGA0104402
18326 mk_peak<-names(which.min(pv2))
map[mk_peak,] chr pos
ALGA0104402 6 136084448#select a region
around_peak<-filter(map,chr==6,pos<map[mk_peak,"pos"]+5000000,pos>map[mk_peak,"pos"]-5000000)
dim(around_peak)[1] 184 2manhattan_plot(pvalues = qv2$qvalues[rownames(around_peak)],map = around_peak)
head(around_peak) chr pos
ALGA0122487 6 131184737
ALGA0117684 6 131254700
ALGA0109995 6 131322492
ALGA0118586 6 131345131
MARC0082363 6 131496812
ALGA0117317 6 131527274tail(around_peak) chr pos
ASGA0029822 6 140724273
ALGA0037308 6 140754146
ALGA0113744 6 140776171
H3GA0019050 6 140804651
ALGA0037315 6 140829602
ALGA0037327 6 141025042IC_90<-IC_gwa(gb2,Z = Z,map = map,
rng =c(6,first(around_peak)[,"pos"],last(around_peak)[,"pos"]),
peak=map[mk_peak,"pos"],
iter=100,alpha = 0.9
)
IC_90 chr pos exact
ASGA0029653 6 134141272 134112676
ALGA0104402 6 136084448 136084448
DIAS0000053 6 137985648 138056220first(around_peak) chr pos
1 6 131184737last(around_peak) chr pos
1 6 141025042ICN_90<-IC_gwa(gb2,Z = Z,map = map,
rng =c(6,first(around_peak)[,"pos"],last(around_peak)[,"pos"]),
peak=map[mk_peak,"pos"],
iter=100,alpha = 0.9,NP = TRUE
)
ICN_90 chr pos exact
ALGA0037046 6 132322578 132322578
ALGA0104402 6 136084448 136084448
ALGA0104402.1 6 136084448 136084448map90<-filter(map,chr==6,pos>=ICN_90[1,3],pos<=ICN_90[3,3])
map90 chr pos
ALGA0037046 6 132322578
ALGA0037050 6 132420958
DRGA0006914 6 132505693
ASGA0029635 6 132543418
DRGA0006919 6 132565019
DRGA0006922 6 132577848
ALGA0037060 6 132598457
DRGA0006924 6 132658334
H3GA0018927 6 132687367
DRGA0006936 6 132729420
DRGA0006937 6 132741221
DRGA0006938 6 132753633
DRGA0006939 6 132769111
ALGA0037068 6 132794256
MARC0068815 6 132846088
ASGA0029645 6 132874648
ASGA0029640 6 132912255
ALGA0110589 6 133030240
MARC0064854 6 133060727
ALGA0115255 6 133084549
ALGA0117544 6 133090986
ALGA0115484 6 133181146
ALGA0103118 6 133185463
ALGA0115991 6 133315760
ASGA0029647 6 133682098
ALGA0037081 6 133743176
ALGA0037083 6 133775504
ASGA0029648 6 133801308
ASGA0029650 6 133839515
M1GA0008917 6 133885460
ASGA0029651 6 133929215
ASGA0029653 6 134141272
INRA0022476 6 134169998
ALGA0037096 6 134197060
ALGA0037101 6 134240349
ASGA0029661 6 134274826
ALGA0037105 6 134307926
DRGA0006951 6 134411824
ASGA0029677 6 134444944
MARC0051254 6 134484418
ASGA0029673 6 134552576
ALGA0112986 6 134616728
ALGA0109354 6 134685885
ASGA0106351 6 134691896
ALGA0107498 6 134736061
ASGA0096926 6 134750885
MARC0083918 6 134895009
MARC0101013 6 134910792
H3GA0053642 6 134996322
INRA0022506 6 135074196
ASGA0029689 6 135089645
ALGA0037123 6 135116578
MARC0025443 6 135143210
MARC0026985 6 135167453
ALGA0037119 6 135188378
ASGA0029681 6 135208598
ASGA0089937 6 135325751
ASGA0093565 6 135424176
MARC0043232 6 135447118
ASGA0086975 6 135523169
ASGA0029693 6 135703082
ALGA0037129 6 135815108
ALGA0037131 6 135878484
H3GA0053839 6 136045327
ALGA0122657 6 136078566
ALGA0104402 6 136084448position_L<-match(rownames(IC_90)[1],rownames(around_peak))
position_p<-match(rownames(IC_90)[2],rownames(around_peak))
position_R<-match(rownames(IC_90)[3],rownames(around_peak))
manhattan_plot(pvalues = qv2$qvalues[rownames(around_peak)],map = around_peak)+
abline(v=position_L,lty="dashed")+abline(v=position_R,lty="dashed")+abline(v=position_p,col="red")
integer(0)position_L<-match(rownames(ICN_90)[1],rownames(around_peak))
position_p<-match(rownames(ICN_90)[2],rownames(around_peak))
position_R<-match(rownames(ICN_90)[3],rownames(around_peak))
manhattan_plot(pvalues = qv2$qvalues[rownames(around_peak)],map = around_peak)+
abline(v=position_L,lty="dashed")+abline(v=position_p,col="red")+abline(v=position_R,lty="dashed")
integer(0)Some questions to ask yourself before using these functions: 1) what factors affect compute time? Sample size? Number of SNP? Search interval for IC?