Gwa from GBLUP

AnS 656

Author

Austin Putz and Juan Steibel

Published

March 19, 2024

Main Path: “C:/Users/jsteibel/OneDrive/Documents/job/ANS 656/hw/”


Packages

Current Package: DT 
Current Package: MASS 
Current Package: Matrix 
Current Package: GGally 
Loading required package: ggplot2
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
Current Package: plotly 

Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':

    last_plot
The following object is masked from 'package:MASS':

    select
The following object is masked from 'package:stats':

    filter
The following object is masked from 'package:graphics':

    layout
Current Package: feather 
Current Package: arrow 

Attaching package: 'arrow'
The following objects are masked from 'package:feather':

    read_feather, write_feather
The following object is masked from 'package:plotly':

    schema
The following object is masked from 'package:utils':

    timestamp
Current Package: pedigreemm 
Loading required package: lme4
Current Package: BGLR 
Current Package: regress 
Current Package: gwaR 
Current Package: lme4 
Current Package: lmerTest 

Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':

    lmer
The following object is masked from 'package:stats':

    step
Current Package: car 
Loading required package: carData
Current Package: emmeans 

Attaching package: 'emmeans'
The following object is masked from 'package:GGally':

    pigs
Current Package: effects 
lattice theme set by effectsTheme()
See ?effectsTheme for details.
Current Package: rsq 
Current Package: performance 
Current Package: coda 
Current Package: pracma 

Attaching package: 'pracma'
The following object is masked from 'package:car':

    logit
The following object is masked from 'package:lmerTest':

    rand
The following object is masked from 'package:arrow':

    null
The following objects are masked from 'package:Matrix':

    expm, lu, tril, triu
Current Package: qvalue 
Current Package: Matrix 
Current Package: skimr 
Current Package: data.table 
Current Package: readxl 
Current Package: lubridate 

Attaching package: 'lubridate'
The following objects are masked from 'package:data.table':

    hour, isoweek, mday, minute, month, quarter, second, wday, week,
    yday, year
The following object is masked from 'package:arrow':

    duration
The following objects are masked from 'package:base':

    date, intersect, setdiff, union
Current Package: broom 
Current Package: stringr 
Current Package: tidyverse 
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr   1.1.4     ✔ readr   2.1.5
✔ forcats 1.0.0     ✔ tibble  3.2.1
✔ purrr   1.0.2     ✔ tidyr   1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between()      masks data.table::between()
✖ purrr::cross()        masks pracma::cross()
✖ lubridate::duration() masks arrow::duration()
✖ tidyr::expand()       masks Matrix::expand()
✖ dplyr::filter()       masks plotly::filter(), stats::filter()
✖ dplyr::first()        masks data.table::first()
✖ lubridate::hour()     masks data.table::hour()
✖ lubridate::isoweek()  masks data.table::isoweek()
✖ dplyr::lag()          masks stats::lag()
✖ dplyr::last()         masks data.table::last()
✖ lubridate::mday()     masks data.table::mday()
✖ lubridate::minute()   masks data.table::minute()
✖ lubridate::month()    masks data.table::month()
✖ tidyr::pack()         masks Matrix::pack()
✖ lubridate::quarter()  masks data.table::quarter()
✖ dplyr::recode()       masks car::recode()
✖ lubridate::second()   masks data.table::second()
✖ dplyr::select()       masks plotly::select(), MASS::select()
✖ purrr::some()         masks car::some()
✖ purrr::transpose()    masks data.table::transpose()
✖ tidyr::unpack()       masks Matrix::unpack()
✖ lubridate::wday()     masks data.table::wday()
✖ lubridate::week()     masks data.table::week()
✖ lubridate::yday()     masks data.table::yday()
✖ lubridate::year()     masks data.table::year()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

initialize

          used  (Mb) gc trigger  (Mb) max used (Mb)
Ncells 2956191 157.9    4718254 252.0  4718254  252
Vcells 5510610  42.1   10146329  77.5  8388451   64
Folder 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  3050113 162.9    4718254 252.0  4718254  252
Vcells 51851656 395.6   76116327 580.8 51903453  396
# 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     3
dim(data.geno)
[1]  1015 45330
dim(data.pheno)
[1] 921  10
id_geno<-data.geno$ID
id_pheno<-data.pheno$ID
skim(data.pheno)
Data summary
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] 0.000000000000002404645
Z<-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)
NULL
dim(G)
[1] 1015 1015
dim(data.geno)
[1]  1015 45330
rownames(G)<-colnames(G)<-id_geno

Notice: 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 warning
Warning: 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)
gb
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.00000000000000000
SexM          -7.929912  19.919779 -0.3980924 0.69056210087545322
Age           -6.644845   1.013622 -6.5555452 0.00000000005543899

estimated variance components:
   Estimate StdError  prop.var
G   565.473 1447.148 0.1016316
In 4998.477 1460.166 0.8983684
summary(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.00000000000000000
SexM          -7.929912  19.919779 -0.3980924 0.69056210087545322
Age           -6.644845   1.013622 -6.5555452 0.00000000005543899

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.056199
varcomp(gb)
   Estimate StdError  prop.var        se
G   565.473 1447.148 0.1016316 0.2546509
In 4998.477 1460.166 0.8983684 0.3736658

Now, 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] 0.000184632 -0.00000095 -0.000001785 0.000001887 0.000000605 ...
  ..- 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: 0x00000258dc245bd8> 
  ..$ Vformula:Class 'formula'  language ~G + In
  .. .. ..- attr(*, ".Environment")=<environment: 0x00000258dc18f0e8> 
 $ 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)
gb
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.000000000000000000
SexM          37.50750  5.3626256   6.994242 0.000000000002666978
Age           -5.41047  0.2726298 -19.845482 0.000000000000000000

estimated variance components:
   Estimate StdError  prop.var
G  2248.512 442.0670 0.3305127
In 4554.590 289.3812 0.6694873
summary(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.000000000000000000
SexM          37.50750  5.3626256   6.994242 0.000000000002666978
Age           -5.41047  0.2726298 -19.845482 0.000000000000000000

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.175948
varcomp(gb)
   Estimate StdError  prop.var         se
G  2248.512 442.0670 0.3305127 0.05455006
In 4554.590 289.3812 0.6694873 0.06232086
str(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] 0.00018168 -0.00000586 -0.000006 -0.00000746 -0.00000117 ...
  ..- 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: 0x00000258be7834c0> 
  ..$ Vformula:Class 'formula'  language ~G + In
  .. .. ..- attr(*, ".Environment")=<environment: 0x00000258be741f70> 
 $ 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.159150
summary(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:

          <0.0001 <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  1450
qqgplot(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 565627
map<-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

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.08319594
head(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.563
pv3<-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   2
manhattan_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 131527274
tail(around_peak)
            chr       pos
ASGA0029822   6 140724273
ALGA0037308   6 140754146
ALGA0113744   6 140776171
H3GA0019050   6 140804651
ALGA0037315   6 140829602
ALGA0037327   6 141025042
IC_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 134096745
ALGA0104402   6 136084448 136084448
DIAS0000379   6 138133920 138072151
first(around_peak)
  chr       pos
1   6 131184737
last(around_peak)
  chr       pos
1   6 141025042
ICN_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 136084448
map90<-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 136084448
position_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?