By default, Haplin estimates the relative risk (RR) of a phenotype
associated with a haplotype, based on triad or dyad genotype data. The
output of the genDataPreprocess
function (or
genDataLoad
) is used to run the analysis. Haplin analysis
is run by the single command:
haplin( my.prepared.gen.data )
CAUTION: this command will try to provide estimates based on all the markers in the data object! Therefore, if you have a large dataset, such as from GWAS analysis, first try running a scan over the markers with a small window size, to determine where to focus your subsequent more detailed analysis:
haplinSlide( my.prepared.gen.data, use.missing = TRUE, winlength = 3 )
This performs haplin analysis on the marker window of length given by
the winlegth
argument above in order to search for the most
significant regions in the dataset.
For more options and examples of how to run Haplin, see below or the haplin help file, obtained by writing in R:
?haplin
and
?haplinSlide
To test that Haplin runs properly, you can use the exemplary data provided with Haplin.
Here, the data includes only genotypes and the analysis is performed on all markers:
<- system.file( "extdata", package = "Haplin" )
dir.exmpl <- paste0( dir.exmpl, "/HAPLIN.trialdata.txt" )
exemplary.file1
<- genDataRead( file.in = exemplary.file1, file.out = "trial_data1",
trial.data1 dir.out = ".", format = "haplin", n.vars = 0 )
## Reading the data in chunks...
## -- chunk 1 --
## -- chunk 2 --
## ... done reading.
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Preparing data...
## ... done preparing
## Saving data...
## ... saved to files: ./trial_data1_gen.ffData, ./trial_data1_gen.RData
<- genDataPreprocess( data.in = trial.data1, design = "triad",
trial.data1.prep file.out = "trial_data1_prep", dir.out = "." )
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Recoding genetic data (no. of loci: 2)...
## ...running on only one CPU core! This may take some time...
## ...checking alleles per SNP...
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be8376f7174.ff
## ...done, all alleles: 1 2 3 4 C T
## ...recoding SNPs...
## ...done
## Saving data...
## ... saved to files: ./trial_data1_prep_gen.ffData , ./trial_data1_prep_gen.RData
haplin( trial.data1.prep, use.missing = TRUE, maternal = TRUE )
##
## ## HAPLIN, VERSION 7.3.2 ##
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be8694d1fee.ff
## There were 5 rows with missing data
## All rows retained in analysis
## No lines contained Mendelian inconsistencies
##
## Running EM for preliminary estimates of haplotype frequencies... Done
## Removing unused haplotypes... Done
##
## Using EM to estimate model with no effect:
## EM iter: 1 |GLM deviance: 7.77156e-15 |Coefficients: -2.78681e-17 -1.39629e-18 1.77335e-18 3.34327e-18 1.13512e-18
## EM iter: 2 |GLM deviance: 345.258 |Coefficients: -0.899013 0.0911016 -1.40653 0.684214 -1.50184
## EM iter: 3 |GLM deviance: 344.375 |Coefficients: -0.90182 0.0914246 -1.41264 0.685363 -1.50184
## EM iter: 4 |GLM deviance: 344.373 |Coefficients: -0.901826 0.0914253 -1.41265 0.685366 -1.50184
## EM iter: 5 |GLM deviance: 344.373 |Coefficients: -0.901826 0.0914253 -1.41265 0.685366 -1.50184
## EM iter: 6 |GLM deviance: 344.373 |Coefficients: -0.901826 0.0914253 -1.41265 0.685366 -1.50184
##
## Done
##
## Using EM to estimate full model:
## EM iter: 1 |GLM deviance: 5.55112e-15 |Coefficients: -1.02642e-17 1.20332e-17 3.19889e-18 3.72728e-18 4.71522e-18 7.37845e-20 -1.40935e-17 -4.03287e-18 -2.10359e-18 -3.09759e-17 -3.61426e-17 2.63283e-18 -8.86023e-18 -2.77969e-18 -3.27258e-18 -7.33247e-18 -2.91917e-18 2.35916e-19 -2.40023e-17 -4.22268e-17 1.49274e-17 -7.57695e-18 -4.42036e-18
## EM iter: 2 |GLM deviance: 324.308 |Coefficients: -0.786562 0.234966 -1.23777 0.627763 -1.46355 -0.603638 -0.37483 -0.782319 -0.137819 0.942851 -0.262305 0.700751 -0.274264 0.972414 0.30563 0.245949 0.28865 0.00556478 -0.468796 -0.353194 0.717524 0.521617 -14.3119
## EM iter: 3 |GLM deviance: 324.026 |Coefficients: -0.785755 0.238644 -1.24484 0.632153 -1.46081 -0.609407 -0.3795 -0.77606 -0.14077 0.953784 -0.258198 0.703041 -0.281416 0.974712 0.301628 0.239555 0.284818 0.00131739 -0.462487 -0.346801 0.728309 0.512107 -14.3091
## EM iter: 4 |GLM deviance: 324.025 |Coefficients: -0.785716 0.238678 -1.24486 0.632184 -1.46078 -0.609489 -0.379552 -0.776042 -0.140807 0.953882 -0.258161 0.70306 -0.28148 0.974733 0.301582 0.239504 0.284766 0.00128555 -0.46244 -0.346727 0.728394 0.512041 -14.3091
## EM iter: 5 |GLM deviance: 324.025 |Coefficients: -0.785716 0.238678 -1.24486 0.632185 -1.46078 -0.60949 -0.379553 -0.776042 -0.140807 0.953884 -0.258161 0.70306 -0.281481 0.974733 0.301581 0.239504 0.284766 0.0012854 -0.46244 -0.346726 0.728394 0.51204 -14.3091
## EM iter: 6 |GLM deviance: 324.025 |Coefficients: -0.785716 0.238678 -1.24486 0.632185 -1.46078 -0.60949 -0.379553 -0.776042 -0.140807 0.953884 -0.258161 0.70306 -0.281481 0.974733 0.301581 0.239504 0.284766 0.0012854 -0.46244 -0.346726 0.728394 0.51204 -14.3091
##
## Done
##
## Estimation finished, preparing output... Done
##
## Note: Some relative risk estimates fall outside the default plotting range.
## Consider replotting, with argument "ylim" set wider
##
## #################################
## ----Arguments supplied to haplin in this run:----
##
## filespecs: markers = 1:2
## model: design = "triad", use.missing = TRUE, xchrom = FALSE, comb.sex = "double", maternal = TRUE, poo = FALSE, test.maternal = FALSE, scoretest = "no"
## variables: ccvar = NULL, strata = NULL, sex = NULL
## haplos: reference = "reciprocal", response = "free", threshold = 0.01, max.haplos = NULL, haplo.file = NULL
## control: resampling = "no", max.EM.iter = 50, data.out = "no", verbose = TRUE, printout = TRUE
##
## ----Data summary:----
##
## Number of triads in original file: 249
##
## Accounting for possible loss of triads:
##
## Cause of loss Triads removed Triads remaining
## Missing data 0 249
## Mendelian incons. 0 249
## Unused haplotypes 5 244
##
## Triads remaining for analysis: 244
##
## NOTE: In the following, the most frequent allele
## is printed as upper-case, all others are lower-case
##
## Marker m1 (raw counts):
## Missing alleles: 6
## Allele Frequency Percent
## 1 149 10.0
## 2 412 27.7
## 3 87 5.8
## 4 840 56.5
## total 1488 100.0
## Chi-squared test for HWE, p-value: 0.6585
##
## Marker m2 (raw counts):
## Missing alleles: 4
## Allele Frequency Percent
## C 1397 93.8
## t 93 6.2
## total 1490 100.0
## Chi-squared test for HWE, p-value: 0.5721
##
## --------
## Haplotypes removed because of low frequencies:
## 1-t 2-t 3-t
##
## Haplotypes used in the analysis, with coding:
## 1-C 2-C 3-C 4-C 4-t
## 1 2 3 4 5
##
## ----Estimation results:----
##
## Date of call:
## Tue Aug 20 10:48:35 2024
##
## Number of triads: 244
##
## Number of haplotypes: 5
##
## Haplotype frequencies with 95% confidence intervals:
## Haplotype Frequency(%) lower upper
## 1-C 11.06 8.00 14.97
## 2-C 30.67 25.84 35.94
## 3-C 6.90 4.59 10.32
## 4-C 45.40 39.92 50.83
## 4-t 5.59 3.51 8.68
##
## Single- and double dose effects (Relative Risk) with 95% confidence intervals:
## Reference method: reciprocal
## Response model: free
##
## ----Child haplotypes----
## Haplotype Dose Relative Risk Lower CI Upper CI P-value
## 1-C Single 0.727 0.454 1.16 0.183
## 1-C Double 1.21 0.389 3.84 0.738
##
## 2-C Single 1.04 0.706 1.54 0.837
## 2-C Double 0.641 0.309 1.32 0.226
##
## 3-C Single 0.614 0.343 1.09 0.102
## 3-C Double 0.674 0.0822 5.26 0.713
##
## 4-C Single 1.62 1.03 2.55 0.0405
## 4-C Double 1.88 1.05 3.39 0.0354
##
## 4-t Single 1.25 0.712 2.26 0.442
## 4-t Double 3.53 0.709 18.2 0.121
##
## ----Maternal haplotypes----
## Haplotype Dose Relative Risk Lower CI Upper CI P-value
## 1-C Single 1.17 0.723 1.86 0.51
## 1-C Double 0.896 0.191 4.31 0.889
##
## 2-C Single 1.07 0.711 1.6 0.754
## 2-C Double 0.875 0.436 1.75 0.702
##
## 3-C Single 1.14 0.654 2.02 0.639
## 3-C Double 2.85 0.562 13.8 0.21
##
## 4-C Single 0.799 0.529 1.2 0.269
## 4-C Double 1.05 0.619 1.78 0.864
##
## 4-t Single 0.832 0.464 1.5 0.545
## 4-t Double 3.02e-05 0 Inf 1
##
##
## Overall test for difference between null model (no effects) and full model:
## ------------
## LIKELIHOOD RATIO TEST:
## Loglike null model: -1228.0728
## Loglike full model: -1217.7436
## df: 18.0000
## Likelihood ratio p-value: 0.2970
##
## (NOTE: The test may be sensitive to rare haplotypes)
First, the information about the calculation process is printed:
use.missing = F
);Next, the default is to print the summary of the results:
As you can see, a lot of information is printed out by default. One
can change it by setting options verbose
and/or
printout
to FALSE
. Moreover, it’s usually
quite useful to save the results into an object:
<- haplin( trial.data1.prep, use.missing = TRUE, maternal = TRUE,
my.results verbose = FALSE, printout = FALSE )
my.results
## This is the result of a haplin run.
## Number of data lines used: 244 | Number of haplotypes used: 5
## Please use the "summary", "plot", "haptable" or "output" functions to obtain
## more details.
To access the details of the results, we can use the following
functions with the saved object my.results
as the
argument:
summary
, which outputs all the results in a nicely
formatted text (as outputted above) or,haptable
, which gives the same data only in a matrix
format,plot
, which replots the figure with the results,output
, which saves the results to text and JPG
files.This example shows analysis of data where apart from genotypes, there
are two covariate columns (n.vars = 2
), one coding for the
case-control status (ccvar = 2
). In the analysis, we take
into account all the available data, imputing the parts that are missing
(use.missing = TRUE
). The estimation will be based on the
dose-response model (response = "mult"
).
<- paste0( dir.exmpl, "/HAPLIN.trialdata2.txt" )
exemplary.file2
<- genDataRead( file.in = exemplary.file2, file.out = "trial_data2",
trial.data2 dir.out = ".", format = "haplin", allele.sep = "", n.vars = 2 )
## The format of the file is 'haplin' with covariate data but no names of the covariate data is given. Will generate dummy names.
## Reading the data in chunks...
## -- chunk 1 --
## -- chunk 2 --
## ... done reading.
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Preparing data...
## ... done preparing
## Saving data...
## ... saved to files: ./trial_data2_gen.ffData, ./trial_data2_gen.RData
<- genDataPreprocess( data.in = trial.data2, design = "triad",
trial.data2.prep file.out = "trial_data2_prep", dir.out = "." )
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Recoding covariate data...
## ...done
## Recoding genetic data (no. of loci: 2)...
## ...running on only one CPU core! This may take some time...
## ...checking alleles per SNP...
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be81ace3f45.ff
## ...done, all alleles: A G B
## ...recoding SNPs...
## ...done
## Saving data...
## ... saved to files: ./trial_data2_prep_gen.ffData , ./trial_data2_prep_gen.RData
haplin( trial.data2.prep, use.missing = TRUE, ccvar = 2, design =
"cc.triad", reference = "ref.cat", response = "mult" )
##
## ## HAPLIN, VERSION 7.3.2 ##
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be832d73de0.ff
## There were 508 rows with missing data
## All rows retained in analysis
## The following 5 data lines were dropped due to Mendelian inconsistencies:
## 240 257 302 320 780
##
## Running EM for preliminary estimates of haplotype frequencies... Done
## Removing unused haplotypes... Done
##
## Note:
## The selected case/control variable is column 2: 'cov.2.c'
## The following case/control coding and frequencies have been found:
## controls(0): 762, cases(1): 377
##
## Using EM to estimate model with no effect:
## EM iter: 1 |GLM deviance: 0 |Coefficients: 1.22453e-17 -2.87212e-18 -1.60405e-18 6.10876e-19
## EM iter: 2 |GLM deviance: 359.799 |Coefficients: -0.713055 -1.18646 0.320706 1.27225
## EM iter: 3 |GLM deviance: 104.965 |Coefficients: -0.713055 -1.90555 0.224806 1.34845
## EM iter: 4 |GLM deviance: 101.348 |Coefficients: -0.713055 -2.02624 0.2036 1.3596
## EM iter: 5 |GLM deviance: 101.879 |Coefficients: -0.713055 -2.04019 0.19944 1.36137
## EM iter: 6 |GLM deviance: 101.978 |Coefficients: -0.713055 -2.04171 0.198655 1.36167
## EM iter: 7 |GLM deviance: 101.994 |Coefficients: -0.713055 -2.04188 0.198509 1.36172
## EM iter: 8 |GLM deviance: 101.997 |Coefficients: -0.713055 -2.04189 0.198483 1.36173
## EM iter: 9 |GLM deviance: 101.998 |Coefficients: -0.713055 -2.0419 0.198478 1.36173
## EM iter: 10 |GLM deviance: 101.998 |Coefficients: -0.713055 -2.0419 0.198477 1.36173
## EM iter: 11 |GLM deviance: 101.998 |Coefficients: -0.713055 -2.0419 0.198477 1.36173
## EM iter: 12 |GLM deviance: 101.998 |Coefficients: -0.713055 -2.0419 0.198477 1.36173
## EM iter: 13 |GLM deviance: 101.998 |Coefficients: -0.713055 -2.0419 0.198477 1.36173
##
## Done
##
## Using EM to estimate full model:
## EM iter: 1 |GLM deviance: 0 |Coefficients: 1.18534e-17 -2.87051e-18 -1.76678e-18 6.01548e-19 2.8841e-19 7.6097e-19
## EM iter: 2 |GLM deviance: 310.934 |Coefficients: -0.808547 -1.05734 0.26884 1.27994 -1.28603 0.328171
## EM iter: 3 |GLM deviance: 89.7213 |Coefficients: -0.893282 -1.79929 0.167796 1.36204 -0.501711 0.387945
## EM iter: 4 |GLM deviance: 86.3708 |Coefficients: -0.905895 -1.9526 0.145439 1.37468 -0.342099 0.402278
## EM iter: 5 |GLM deviance: 87.1614 |Coefficients: -0.908033 -1.97423 0.140962 1.37675 -0.319404 0.405527
## EM iter: 6 |GLM deviance: 87.3163 |Coefficients: -0.908428 -1.97705 0.140093 1.3771 -0.316378 0.40624
## EM iter: 7 |GLM deviance: 87.3418 |Coefficients: -0.908505 -1.97742 0.139926 1.37716 -0.315974 0.406393
## EM iter: 8 |GLM deviance: 87.346 |Coefficients: -0.90852 -1.97747 0.139895 1.37717 -0.315919 0.406426
## EM iter: 9 |GLM deviance: 87.3468 |Coefficients: -0.908524 -1.97747 0.139889 1.37718 -0.315911 0.406432
## EM iter: 10 |GLM deviance: 87.3469 |Coefficients: -0.908524 -1.97747 0.139888 1.37718 -0.31591 0.406434
## EM iter: 11 |GLM deviance: 87.3469 |Coefficients: -0.908524 -1.97747 0.139888 1.37718 -0.315909 0.406434
## EM iter: 12 |GLM deviance: 87.3469 |Coefficients: -0.908524 -1.97747 0.139887 1.37718 -0.315909 0.406434
## EM iter: 13 |GLM deviance: 87.3469 |Coefficients: -0.908524 -1.97747 0.139887 1.37718 -0.315909 0.406434
##
## Done
##
## Estimation finished, preparing output... Done
##
## #################################
## ----Arguments supplied to haplin in this run:----
##
## filespecs: markers = 1:2
## model: design = "cc.triad", use.missing = TRUE, xchrom = FALSE, comb.sex = "double", maternal = FALSE, poo = FALSE, test.maternal = FALSE, scoretest = "no"
## variables: ccvar = 2, strata = NULL, sex = NULL
## haplos: reference = "ref.cat", response = "mult", threshold = 0.01, max.haplos = NULL, haplo.file = NULL
## control: resampling = "no", max.EM.iter = 50, data.out = "no", verbose = TRUE, printout = TRUE
##
## ----Data summary:----
##
## Number of triads in original file: 1139
##
## Accounting for possible loss of triads:
##
## Cause of loss Triads removed Triads remaining
## Missing data 0 1139
## Mendelian incons. 5 1134
## Unused haplotypes 0 1134
##
## Triads remaining for analysis: 1134
##
## NOTE: In the following, the most frequent allele
## is printed as upper-case, all others are lower-case
##
## Marker m1 (raw counts):
## Missing alleles: 1424
## Allele Frequency Percent
## a 1290 23.8
## G 4120 76.2
## total 5410 100.0
## Chi-squared test for HWE, p-value: 0.6875
##
## Marker m2 (raw counts):
## Missing alleles: 954
## Allele Frequency Percent
## a 142 2.4
## B 5738 97.6
## total 5880 100.0
## Chi-squared test for HWE, p-value: 0.01014
##
## --------
## Haplotypes removed because of low frequencies:
## a-a
##
## Haplotypes used in the analysis, with coding:
## G-a a-B G-B
## 1 2 3
##
## ----Estimation results:----
##
## Date of call:
## Tue Aug 20 10:48:37 2024
##
## Number of triads: 1134
##
## Number of haplotypes: 3
##
## Haplotype frequencies with 95% confidence intervals:
## Haplotype Frequency(%) lower upper
## G-a 2.64 2.15 3.23
## a-B 21.88 20.47 23.40
## G-B 75.46 73.90 76.96
##
## Single- and double dose effects (Relative Risk) with 95% confidence intervals:
## Reference method: ref.cat
## Reference category: 3 (Haplotype G-B)
## Response model: mult
##
## ----Child haplotypes----
## Haplotype Dose Relative Risk Lower CI Upper CI P-value
## G-a Single 0.729 0.408 1.33 0.295
## G-a Double 0.532 0.166 1.76 0.295
##
## a-B Single 1.5 1.23 1.83 5.35e-05
## a-B Double 2.24 1.52 3.35 5.35e-05
##
## G-B Single REF
## G-B Double 1 1 1 NA
##
##
## Overall test for difference between null model (no effects) and full model:
## ------------
## LIKELIHOOD RATIO TEST:
## Loglike null model: -3117.5699752
## Loglike full model: -3108.4934459
## df: 2.0000000
## Likelihood ratio p-value: 0.0001143
##
## (NOTE: The test may be sensitive to rare haplotypes)
This is very similar to running haplin
, but the results
are a list of haptable. This type of analysis is usually performed when
we have much bigger data. Let’s read a proper .ped file:
<- paste0( dir.exmpl, "/exmpl_data.ped" )
exemplary.file3
<- genDataRead( exemplary.file3, file.out = "hapSlide_data",
hapSlide.data format = "ped" )
## 'n.vars' was not given explicitly and will be set to 6 based on the format given.
## Reading the data in chunks...
## -- chunk 1 --
## -- chunk 2 --
## ... done reading.
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Preparing data...
## ... done preparing
## Saving data...
## ... saved to files: ./hapSlide_data_gen.ffData, ./hapSlide_data_gen.RData
hapSlide.data
## This is raw genetic data read in through genDataRead.
##
## It contains the following parts:
## cov.data, gen.data, aux
##
## with following dimensions:
## - covariate variables = id.fam, id.c, id.f, id.m, sex, cc
## (total 6 covariate variables),
## - number of markers = 429 ,
## - number of data lines = 1659
The preprocessing of this file would take a while. Let’s focus then on first 100 markers:
<- genDataGetPart( hapSlide.data, design = "cc",
hapSlide.subset.data markers = 1:100, file.out = "hapSlide_subset_data" )
## Provided arguments:
## --- chosen markers: 1, 2, 3, 4, 5, 6 ... 95, 96, 97, 98, 99, 100
## INFO: Will select 1659 rows and 200 columns.
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be85bc869d7.ff
## Saving data...
## ... saved to files: ./hapSlide_subset_data_gen.ffData, ./hapSlide_subset_data_gen.RData
<- genDataPreprocess( hapSlide.subset.data,
hapSlide.subset.data.prep design = "cc.triad", file.out = "hapSlide_subset_prep" )
## Reading the marker names...
## Warning: No map file given, map file empty or the number of map file rows not
## equal to the number of markers in data; will generate dummy marker names.
## ...done
## Converting PED format to internal haplin...
## Creating unique IDs for individuals...
## ...done.
##
## Checking family and id variables...
## Warning: Found family size larger than 3! Will extract trios from general
## pedigree.
## Warning: 1100 mother code(s) and 1100 father code(s) refer to non-existing
## individuals and have been set to missing
## Sorting and re-coding families...
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be83e8c6699.ff
## Recoding covariate data...
## ...done
## Recoding genetic data (no. of loci: 100)...
## ...running on only one CPU core! This may take some time...
## ...checking alleles per SNP...
## ...done, all alleles: C G A T
## ...recoding SNPs...
## ...done
## Saving data...
## ... saved to files: ./hapSlide_subset_prep_gen.ffData , ./hapSlide_subset_prep_gen.RData
We run the haplinSlide analysis:
<- haplinSlide( hapSlide.subset.data.prep, use.missing = TRUE,
hapSlide.res1 ccvar = 10, design = "cc.triad", reference = "ref.cat", response = "mult" )
## INFO: No parallel environment selected. Will run in sequential mode.
## Running Haplin on Window 'm1' (1/100)...:
## opening ff C:/Users/hkgje/AppData/Local/Temp/Rtmpo5NE6q/ff/ff8be84d77eec.ff
## done
## Running Haplin on Window 'm2' (2/100)...: done
## Running Haplin on Window 'm3' (3/100)...: done
## Running Haplin on Window 'm4' (4/100)...: done
## Running Haplin on Window 'm5' (5/100)...: done
## Running Haplin on Window 'm6' (6/100)...: done
## Running Haplin on Window 'm7' (7/100)...: done
## Running Haplin on Window 'm8' (8/100)...: done
## Running Haplin on Window 'm9' (9/100)...: done
## Running Haplin on Window 'm10' (10/100)...: done
## Running Haplin on Window 'm11' (11/100)...: done
## Running Haplin on Window 'm12' (12/100)...: done
## Running Haplin on Window 'm13' (13/100)...: done
## Running Haplin on Window 'm14' (14/100)...: done
## Running Haplin on Window 'm15' (15/100)...: done
## Running Haplin on Window 'm16' (16/100)...: done
## Running Haplin on Window 'm17' (17/100)...: done
## Running Haplin on Window 'm18' (18/100)...: done
## Running Haplin on Window 'm19' (19/100)...: done
## Running Haplin on Window 'm20' (20/100)...: done
## Running Haplin on Window 'm21' (21/100)...: done
## Running Haplin on Window 'm22' (22/100)...: done
## Running Haplin on Window 'm23' (23/100)...: done
## Running Haplin on Window 'm24' (24/100)...: done
## Running Haplin on Window 'm25' (25/100)...: done
## Running Haplin on Window 'm26' (26/100)...: done
## Running Haplin on Window 'm27' (27/100)...: done
## Running Haplin on Window 'm28' (28/100)...: done
## Running Haplin on Window 'm29' (29/100)...: done
## Running Haplin on Window 'm30' (30/100)...: done
## Running Haplin on Window 'm31' (31/100)...: done
## Running Haplin on Window 'm32' (32/100)...: done
## Running Haplin on Window 'm33' (33/100)...: done
## Running Haplin on Window 'm34' (34/100)...: done
## Running Haplin on Window 'm35' (35/100)...: done
## Running Haplin on Window 'm36' (36/100)...: done
## Running Haplin on Window 'm37' (37/100)...: done
## Running Haplin on Window 'm38' (38/100)...: done
## Running Haplin on Window 'm39' (39/100)...: done
## Running Haplin on Window 'm40' (40/100)...: done
## Running Haplin on Window 'm41' (41/100)...: done
## Running Haplin on Window 'm42' (42/100)...: done
## Running Haplin on Window 'm43' (43/100)...: done
## Running Haplin on Window 'm44' (44/100)...: done
## Running Haplin on Window 'm45' (45/100)...: done
## Running Haplin on Window 'm46' (46/100)...: done
## Running Haplin on Window 'm47' (47/100)...: done
## Running Haplin on Window 'm48' (48/100)...: done
## Running Haplin on Window 'm49' (49/100)...: done
## Running Haplin on Window 'm50' (50/100)...: done
## Running Haplin on Window 'm51' (51/100)...: done
## Running Haplin on Window 'm52' (52/100)...: done
## Running Haplin on Window 'm53' (53/100)...: done
## Running Haplin on Window 'm54' (54/100)...: done
## Running Haplin on Window 'm55' (55/100)...: done
## Running Haplin on Window 'm56' (56/100)...: done
## Running Haplin on Window 'm57' (57/100)...: done
## Running Haplin on Window 'm58' (58/100)...: done
## Running Haplin on Window 'm59' (59/100)...: done
## Running Haplin on Window 'm60' (60/100)...: done
## Running Haplin on Window 'm61' (61/100)...: done
## Running Haplin on Window 'm62' (62/100)...: done
## Running Haplin on Window 'm63' (63/100)...: done
## Running Haplin on Window 'm64' (64/100)...: done
## Running Haplin on Window 'm65' (65/100)...: done
## Running Haplin on Window 'm66' (66/100)...: done
## Running Haplin on Window 'm67' (67/100)...: done
## Running Haplin on Window 'm68' (68/100)...: done
## Running Haplin on Window 'm69' (69/100)...: done
## Running Haplin on Window 'm70' (70/100)...: done
## Running Haplin on Window 'm71' (71/100)...: done
## Running Haplin on Window 'm72' (72/100)...: done
## Running Haplin on Window 'm73' (73/100)...: done
## Running Haplin on Window 'm74' (74/100)...: done
## Running Haplin on Window 'm75' (75/100)...: done
## Running Haplin on Window 'm76' (76/100)...: done
## Running Haplin on Window 'm77' (77/100)...: done
## Running Haplin on Window 'm78' (78/100)...: done
## Running Haplin on Window 'm79' (79/100)...: done
## Running Haplin on Window 'm80' (80/100)...: done
## Running Haplin on Window 'm81' (81/100)...: done
## Running Haplin on Window 'm82' (82/100)...: done
## Running Haplin on Window 'm83' (83/100)...: done
## Running Haplin on Window 'm84' (84/100)...: done
## Running Haplin on Window 'm85' (85/100)...: done
## Running Haplin on Window 'm86' (86/100)...: done
## Running Haplin on Window 'm87' (87/100)...: done
## Running Haplin on Window 'm88' (88/100)...: done
## Running Haplin on Window 'm89' (89/100)...: done
## Running Haplin on Window 'm90' (90/100)...: done
## Running Haplin on Window 'm91' (91/100)...: done
## Running Haplin on Window 'm92' (92/100)...: done
## Running Haplin on Window 'm93' (93/100)...: done
## Running Haplin on Window 'm94' (94/100)...: done
## Running Haplin on Window 'm95' (95/100)...: done
## Running Haplin on Window 'm96' (96/100)...: done
## Running Haplin on Window 'm97' (97/100)...: done
## Running Haplin on Window 'm98' (98/100)...: done
## Running Haplin on Window 'm99' (99/100)...: done
## Running Haplin on Window 'm100' (100/100)...: done
##
## --- haplinSlide has completed ---
Here, the output is stored as a list of tables, where each table
contains results of the haplin
run on a certain marker or a
set of markers, called a window. The length of the window can be
controlled by the winlength
argument (see the haplinSlide
help for details). Here, the default window length was used, 1
marker.
Due to the size of haplinSlide
results, it is best to
narrow down the visualization of the results. If you have installed the
ggplot2
package, you can use the functions plotPValues
and
plot.haplinSlide
. First, we advise to plot only the
p-values of the estimation we’re interested in, using the
plotPValues
function:
<- plotPValues( hapSlide.res1 ) all.p.values
By default, all the windows are plotted, but we can choose to plot results for the windows only with p-values below a certain threshold:
<- plotPValues( hapSlide.res1, plot.signif.only = TRUE,
all.p.values signif.thresh = 0.25 )
The object all.p.values
holds the table with p-values
for all the windows plotted.
head( all.p.values )
## marker haplos est.name est. lower upper p.value star
## m1 m1 m1 pv.overall 0.15404214 NA NA 0.15404214
## m2 m2 m2 pv.overall 0.05300921 NA NA 0.05300921
## m4 m4 m4 pv.overall 0.07280665 NA NA 0.07280665
## m5 m5 m5 pv.overall 0.03280638 NA NA 0.03280638 *
## m6 m6 m6 pv.overall 0.24185705 NA NA 0.24185705
## m7 m7 m7 pv.overall 0.24607721 NA NA 0.24607721
If we ran analysis with estimation of the maternal or parent-of-origin (PoO) effects, then we can choose to plot these p-values:
<- haplinSlide( hapSlide.subset.data.prep, use.missing = TRUE,
hapSlide.res2 ccvar = 10, design = "cc.triad", poo = TRUE, reference = "ref.cat",
response = "mult" )
## INFO: No parallel environment selected. Will run in sequential mode.
## Running Haplin on Window 'm1' (1/100)...: done
## Running Haplin on Window 'm2' (2/100)...: done
## Running Haplin on Window 'm3' (3/100)...: done
## Running Haplin on Window 'm4' (4/100)...: done
## Running Haplin on Window 'm5' (5/100)...: done
## Running Haplin on Window 'm6' (6/100)...: done
## Running Haplin on Window 'm7' (7/100)...: done
## Running Haplin on Window 'm8' (8/100)...: done
## Running Haplin on Window 'm9' (9/100)...: done
## Running Haplin on Window 'm10' (10/100)...: done
## Running Haplin on Window 'm11' (11/100)...: done
## Running Haplin on Window 'm12' (12/100)...: done
## Running Haplin on Window 'm13' (13/100)...: done
## Running Haplin on Window 'm14' (14/100)...: done
## Running Haplin on Window 'm15' (15/100)...: done
## Running Haplin on Window 'm16' (16/100)...: done
## Running Haplin on Window 'm17' (17/100)...: done
## Running Haplin on Window 'm18' (18/100)...: done
## Running Haplin on Window 'm19' (19/100)...: done
## Running Haplin on Window 'm20' (20/100)...: done
## Running Haplin on Window 'm21' (21/100)...: done
## Running Haplin on Window 'm22' (22/100)...: done
## Running Haplin on Window 'm23' (23/100)...: done
## Running Haplin on Window 'm24' (24/100)...: done
## Running Haplin on Window 'm25' (25/100)...: done
## Running Haplin on Window 'm26' (26/100)...: done
## Running Haplin on Window 'm27' (27/100)...: done
## Running Haplin on Window 'm28' (28/100)...: done
## Running Haplin on Window 'm29' (29/100)...: done
## Running Haplin on Window 'm30' (30/100)...: done
## Running Haplin on Window 'm31' (31/100)...: done
## Running Haplin on Window 'm32' (32/100)...: done
## Running Haplin on Window 'm33' (33/100)...: done
## Running Haplin on Window 'm34' (34/100)...: done
## Running Haplin on Window 'm35' (35/100)...: done
## Running Haplin on Window 'm36' (36/100)...: done
## Running Haplin on Window 'm37' (37/100)...: done
## Running Haplin on Window 'm38' (38/100)...: done
## Running Haplin on Window 'm39' (39/100)...: done
## Running Haplin on Window 'm40' (40/100)...: done
## Running Haplin on Window 'm41' (41/100)...: done
## Running Haplin on Window 'm42' (42/100)...: done
## Running Haplin on Window 'm43' (43/100)...: done
## Running Haplin on Window 'm44' (44/100)...: done
## Running Haplin on Window 'm45' (45/100)...: done
## Running Haplin on Window 'm46' (46/100)...: done
## Running Haplin on Window 'm47' (47/100)...: done
## Running Haplin on Window 'm48' (48/100)...: done
## Running Haplin on Window 'm49' (49/100)...: done
## Running Haplin on Window 'm50' (50/100)...: done
## Running Haplin on Window 'm51' (51/100)...: done
## Running Haplin on Window 'm52' (52/100)...: done
## Running Haplin on Window 'm53' (53/100)...: done
## Running Haplin on Window 'm54' (54/100)...: done
## Running Haplin on Window 'm55' (55/100)...: done
## Running Haplin on Window 'm56' (56/100)...: done
## Running Haplin on Window 'm57' (57/100)...: done
## Running Haplin on Window 'm58' (58/100)...: done
## Running Haplin on Window 'm59' (59/100)...: done
## Running Haplin on Window 'm60' (60/100)...: done
## Running Haplin on Window 'm61' (61/100)...: done
## Running Haplin on Window 'm62' (62/100)...: done
## Running Haplin on Window 'm63' (63/100)...: done
## Running Haplin on Window 'm64' (64/100)...: done
## Running Haplin on Window 'm65' (65/100)...: done
## Running Haplin on Window 'm66' (66/100)...: done
## Running Haplin on Window 'm67' (67/100)...: done
## Running Haplin on Window 'm68' (68/100)...: done
## Running Haplin on Window 'm69' (69/100)...: done
## Running Haplin on Window 'm70' (70/100)...: done
## Running Haplin on Window 'm71' (71/100)...: done
## Running Haplin on Window 'm72' (72/100)...: done
## Running Haplin on Window 'm73' (73/100)...: done
## Running Haplin on Window 'm74' (74/100)...: done
## Running Haplin on Window 'm75' (75/100)...: done
## Running Haplin on Window 'm76' (76/100)...: done
## Running Haplin on Window 'm77' (77/100)...: done
## Running Haplin on Window 'm78' (78/100)...: done
## Running Haplin on Window 'm79' (79/100)...: done
## Running Haplin on Window 'm80' (80/100)...: done
## Running Haplin on Window 'm81' (81/100)...: done
## Running Haplin on Window 'm82' (82/100)...: done
## Running Haplin on Window 'm83' (83/100)...: done
## Running Haplin on Window 'm84' (84/100)...: done
## Running Haplin on Window 'm85' (85/100)...: done
## Running Haplin on Window 'm86' (86/100)...: done
## Running Haplin on Window 'm87' (87/100)...: done
## Running Haplin on Window 'm88' (88/100)...: done
## Running Haplin on Window 'm89' (89/100)...: done
## Running Haplin on Window 'm90' (90/100)...: done
## Running Haplin on Window 'm91' (91/100)...: done
## Running Haplin on Window 'm92' (92/100)...: done
## Running Haplin on Window 'm93' (93/100)...: done
## Running Haplin on Window 'm94' (94/100)...: done
## Running Haplin on Window 'm95' (95/100)...: done
## Running Haplin on Window 'm96' (96/100)...: done
## Running Haplin on Window 'm97' (97/100)...: done
## Running Haplin on Window 'm98' (98/100)...: done
## Running Haplin on Window 'm99' (99/100)...: done
## Running Haplin on Window 'm100' (100/100)...: done
##
## --- haplinSlide has completed ---
<- plotPValues( hapSlide.res2, which.p.val = "poo",
all.p.values2 plot.signif.only = TRUE, signif.thresh = 0.2 )
head( all.p.values2 )
## marker haplos est.name est. lower
## 724 m31 m31:\nc (40.5%),\n ref: G (59.5%) RRcm_RRcf1 1.4903807 0.9919165
## 725 m32 m32:\na (35.9%),\n ref: T (64.1%) RRcm_RRcf1 1.5457860 1.0327751
## 729 m36 m36:\na (45.3%),\n ref: T (54.7%) RRcm_RRcf1 0.7436971 0.4825578
## 736 m46 m46:\nc (37.8%),\n ref: T (62.2%) RRcm_RRcf1 0.6638964 0.4280048
## 740 m50 m50:\na (36%),\n ref: T (64%) RRcm_RRcf1 0.6564842 0.4360486
## 745 m55 m55:\nc (44.2%),\n ref: G (55.8%) RRcm_RRcf1 1.4191366 0.9415219
## upper p.value star
## 724 2.217478 0.05305328
## 725 2.286814 0.03184453 *
## 729 1.128567 0.16509548
## 736 1.017188 0.06214700
## 740 1.002324 0.05056970
## 745 2.111067 0.09117895
NOTE: here, the star
column marks the
significant p-values:
*
marks p-value \(<
0.05\),**
marks p-value \(<
0.01\),***
marks p-value \(<
0.001\).And then, we can plot the relative risks:
plot( hapSlide.res2, plot.signif.only = TRUE )
Here’s also how the results look like when we estimate the maternal effects, with the window length set to two:
<- haplinSlide( hapSlide.subset.data.prep, use.missing = TRUE,
hapSlide.res3 ccvar = 10, design = "cc.triad", maternal = TRUE, reference = "ref.cat",
response = "mult", winlength = 2 )
##
## Important: Remember that SNPs must be in correct physical ordering
## for a sliding window approach to be meaningful!
## INFO: No parallel environment selected. Will run in sequential mode.
## Running Haplin on Window 'm1-m2' (1/99)...: done
## Running Haplin on Window 'm2-m3' (2/99)...: done
## Running Haplin on Window 'm3-m4' (3/99)...: done
## Running Haplin on Window 'm4-m5' (4/99)...: done
## Running Haplin on Window 'm5-m6' (5/99)...: done
## Running Haplin on Window 'm6-m7' (6/99)...: done
## Running Haplin on Window 'm7-m8' (7/99)...: done
## Running Haplin on Window 'm8-m9' (8/99)...: done
## Running Haplin on Window 'm9-m10' (9/99)...: done
## Running Haplin on Window 'm10-m11' (10/99)...: done
## Running Haplin on Window 'm11-m12' (11/99)...: done
## Running Haplin on Window 'm12-m13' (12/99)...: done
## Running Haplin on Window 'm13-m14' (13/99)...: done
## Running Haplin on Window 'm14-m15' (14/99)...: done
## Running Haplin on Window 'm15-m16' (15/99)...: done
## Running Haplin on Window 'm16-m17' (16/99)...: done
## Running Haplin on Window 'm17-m18' (17/99)...: done
## Running Haplin on Window 'm18-m19' (18/99)...: done
## Running Haplin on Window 'm19-m20' (19/99)...: done
## Running Haplin on Window 'm20-m21' (20/99)...: done
## Running Haplin on Window 'm21-m22' (21/99)...: done
## Running Haplin on Window 'm22-m23' (22/99)...: done
## Running Haplin on Window 'm23-m24' (23/99)...: done
## Running Haplin on Window 'm24-m25' (24/99)...: done
## Running Haplin on Window 'm25-m26' (25/99)...: done
## Running Haplin on Window 'm26-m27' (26/99)...: done
## Running Haplin on Window 'm27-m28' (27/99)...: done
## Running Haplin on Window 'm28-m29' (28/99)...: done
## Running Haplin on Window 'm29-m30' (29/99)...: done
## Running Haplin on Window 'm30-m31' (30/99)...: done
## Running Haplin on Window 'm31-m32' (31/99)...: done
## Running Haplin on Window 'm32-m33' (32/99)...: done
## Running Haplin on Window 'm33-m34' (33/99)...: done
## Running Haplin on Window 'm34-m35' (34/99)...: done
## Running Haplin on Window 'm35-m36' (35/99)...: done
## Running Haplin on Window 'm36-m37' (36/99)...: done
## Running Haplin on Window 'm37-m38' (37/99)...: done
## Running Haplin on Window 'm38-m39' (38/99)...: done
## Running Haplin on Window 'm39-m40' (39/99)...: done
## Running Haplin on Window 'm40-m41' (40/99)...: done
## Running Haplin on Window 'm41-m42' (41/99)...: done
## Running Haplin on Window 'm42-m43' (42/99)...: done
## Running Haplin on Window 'm43-m44' (43/99)...: done
## Running Haplin on Window 'm44-m45' (44/99)...: done
## Running Haplin on Window 'm45-m46' (45/99)...: done
## Running Haplin on Window 'm46-m47' (46/99)...: done
## Running Haplin on Window 'm47-m48' (47/99)...: done
## Running Haplin on Window 'm48-m49' (48/99)...: done
## Running Haplin on Window 'm49-m50' (49/99)...: done
## Running Haplin on Window 'm50-m51' (50/99)...: done
## Running Haplin on Window 'm51-m52' (51/99)...: done
## Running Haplin on Window 'm52-m53' (52/99)...: done
## Running Haplin on Window 'm53-m54' (53/99)...: done
## Running Haplin on Window 'm54-m55' (54/99)...: done
## Running Haplin on Window 'm55-m56' (55/99)...: done
## Running Haplin on Window 'm56-m57' (56/99)...: done
## Running Haplin on Window 'm57-m58' (57/99)...: done
## Running Haplin on Window 'm58-m59' (58/99)...: done
## Running Haplin on Window 'm59-m60' (59/99)...: done
## Running Haplin on Window 'm60-m61' (60/99)...: done
## Running Haplin on Window 'm61-m62' (61/99)...: done
## Running Haplin on Window 'm62-m63' (62/99)...: done
## Running Haplin on Window 'm63-m64' (63/99)...: done
## Running Haplin on Window 'm64-m65' (64/99)...: done
## Running Haplin on Window 'm65-m66' (65/99)...: done
## Running Haplin on Window 'm66-m67' (66/99)...: done
## Running Haplin on Window 'm67-m68' (67/99)...: done
## Running Haplin on Window 'm68-m69' (68/99)...: done
## Running Haplin on Window 'm69-m70' (69/99)...: done
## Running Haplin on Window 'm70-m71' (70/99)...: done
## Running Haplin on Window 'm71-m72' (71/99)...: done
## Running Haplin on Window 'm72-m73' (72/99)...: done
## Running Haplin on Window 'm73-m74' (73/99)...: done
## Running Haplin on Window 'm74-m75' (74/99)...: done
## Running Haplin on Window 'm75-m76' (75/99)...: done
## Running Haplin on Window 'm76-m77' (76/99)...: done
## Running Haplin on Window 'm77-m78' (77/99)...: done
## Running Haplin on Window 'm78-m79' (78/99)...: done
## Running Haplin on Window 'm79-m80' (79/99)...: done
## Running Haplin on Window 'm80-m81' (80/99)...: done
## Running Haplin on Window 'm81-m82' (81/99)...: done
## Running Haplin on Window 'm82-m83' (82/99)...: done
## Running Haplin on Window 'm83-m84' (83/99)...: done
## Running Haplin on Window 'm84-m85' (84/99)...: done
## Running Haplin on Window 'm85-m86' (85/99)...: done
## Running Haplin on Window 'm86-m87' (86/99)...: done
## Running Haplin on Window 'm87-m88' (87/99)...: done
## Running Haplin on Window 'm88-m89' (88/99)...: done
## Running Haplin on Window 'm89-m90' (89/99)...: done
## Running Haplin on Window 'm90-m91' (90/99)...: done
## Running Haplin on Window 'm91-m92' (91/99)...: done
## Running Haplin on Window 'm92-m93' (92/99)...: done
## Running Haplin on Window 'm93-m94' (93/99)...: done
## Running Haplin on Window 'm94-m95' (94/99)...: done
## Running Haplin on Window 'm95-m96' (95/99)...: done
## Running Haplin on Window 'm96-m97' (96/99)...: done
## Running Haplin on Window 'm97-m98' (97/99)...: done
## Running Haplin on Window 'm98-m99' (98/99)...: done
## Running Haplin on Window 'm99-m100' (99/99)...: done
##
## --- haplinSlide has completed ---
plot( hapSlide.res3, plot.signif.only = TRUE, signif.thresh = 0.01 )
If our dataset contains a column with environmental exposure of
(usually) the mother, we can analyse it with haplinStrat
,
that calculates relative risks for strata defined by the categories of
the environmental exposure variable. To do that, we use the
strata
argument which points to the column in the dataset
that contains the environmental exposure categories.
<- haplinStrat( data = trial.data2.prep, strata = 1, use.missing = TRUE,
hapStrat.res ccvar = 2, design = "cc.triad", reference = "ref.cat", response = "mult" )
##
## ## Running haplinStrat ##
##
## Selected stratification variable: cov.1.c
## Frequency distribution of stratification variable:
## 0 1 2 3 4
## 565 164 235 122 53
##
## Running Haplin on full data file...Done
##
## Running Haplin on stratum "0"...Done
##
## Running Haplin on stratum "1"...Done
##
## Running Haplin on stratum "2"...Done
##
## Running Haplin on stratum "3"...Done
##
## Running Haplin on stratum "4"...Done
The result of haplinStrat
run is a list of haplin
objects: one for the haplin run on the entire data, and then one for
every stratum.
lapply( hapStrat.res, haptable )
## $all
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 1290/4120 0.68753925 1139 1139 1134
## 2 m2 a/B 142/5738 0.01013796 1139 1139 1134
## 3 <NA> <NA> <NA> NA 1139 1139 1134
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 1134 0.0001143177 G-a 0.02637715 0.02147445
## 2 1134 0.0001143177 a-B 0.21878318 0.20469703
## 3 1134 0.0001143177 G-B 0.75462805 0.73895724
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value RRdd.est.
## 1 0.03228444 - 0.729373 0.4079393 1.327277 2.947870e-01 0.531985
## 2 0.23398201 - 1.498203 1.2322773 1.829178 5.347616e-05 2.244612
## 3 0.76964068 ref 1.000000 1.0000000 1.000000 NA 1.000000
## RRdd.lower RRdd.upper RRdd.p.value
## 1 0.1664145 1.761664 2.947870e-01
## 2 1.5185074 3.345894 5.347616e-05
## 3 1.0000000 1.000000 NA
##
## $`0`
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 611/2051 0.54762530 565 565 563
## 2 m2 a/B 83/2825 0.08595161 565 565 563
## 3 <NA> <NA> <NA> NA 565 565 563
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 563 0.05946876 G-a 0.03004278 0.0229362
## 2 563 0.05946876 a-B 0.21821984 0.1989562
## 3 563 0.05946876 G-B 0.75129396 0.7291091
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value RRdd.est.
## 1 0.03915257 - 0.591275 0.2358643 1.516845 0.26752685 0.3496062
## 2 0.23972144 - 1.372374 0.9963017 1.903200 0.05156035 1.8834096
## 3 0.77218341 ref 1.000000 1.0000000 1.000000 NA 1.0000000
## RRdd.lower RRdd.upper RRdd.p.value
## 1 0.05563195 2.300819 0.26752685
## 2 0.99261701 3.622170 0.05156035
## 3 1.00000000 1.000000 NA
##
## $`1`
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 196/574 0.428881321 164 164 163
## 2 m2 a/B 23/813 0.002084963 164 164 163
## 3 <NA> <NA> <NA> NA 164 164 163
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 163 0.04349637 G-a 0.02585087 0.01500796
## 2 163 0.04349637 a-B 0.21986259 0.18256901
## 3 163 0.04349637 G-B 0.75323598 0.70947240
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value RRdd.est.
## 1 0.0458127 - 1.638403 0.5106674 5.122134 0.40743557 2.684363
## 2 0.2614902 - 1.880685 1.1420259 3.009592 0.01171209 3.536976
## 3 0.7920373 ref 1.000000 1.0000000 1.000000 NA 1.000000
## RRdd.lower RRdd.upper RRdd.p.value
## 1 0.2607812 26.236259 0.40743557
## 2 1.3042231 9.057643 0.01171209
## 3 1.0000000 1.000000 NA
##
## $`2`
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 281/873 0.3411622 235 235 233
## 2 m2 a/B 20/1218 0.6828818 235 235 233
## 3 <NA> <NA> <NA> NA 235 235 233
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 233 0.1878774 G-a 0.01974803 0.01129351
## 2 233 0.1878774 a-B 0.22329923 0.19214328
## 3 233 0.1878774 G-B 0.75651115 0.72208306
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value RRdd.est.
## 1 0.03345271 - 0.6054025 0.1373635 2.731093 0.51436824 0.3665122
## 2 0.25682149 - 1.3965202 0.9498318 2.081453 0.09334688 1.9502688
## 3 0.78902968 ref 1.0000000 1.0000000 1.000000 NA 1.0000000
## RRdd.lower RRdd.upper RRdd.p.value
## 1 0.01886874 7.458870 0.51436824
## 2 0.90218040 4.332446 0.09334688
## 3 1.00000000 1.000000 NA
##
## $`3`
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 137/447 0.5291981 122 122 122
## 2 m2 a/B 12/610 0.7286501 122 122 122
## 3 <NA> <NA> <NA> NA 122 122 122
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 122 0.3866152 G-a 0.02092333 0.01001447
## 2 122 0.3866152 a-B 0.20799739 0.16674019
## 3 122 0.3866152 G-B 0.76898148 0.71823195
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value RRdd.est.
## 1 0.04399696 - 1.071522 0.2162447 5.159043 0.9391847 1.148160
## 2 0.25729337 - 1.464412 0.8531272 2.527055 0.1712284 2.144502
## 3 0.81300688 ref 1.000000 1.0000000 1.000000 NA 1.000000
## RRdd.lower RRdd.upper RRdd.p.value
## 1 0.04676178 26.615728 0.9391847
## 2 0.72782610 6.386005 0.1712284
## 3 1.00000000 1.000000 NA
##
## $`4`
## marker alleles counts HWE.pv Original After.rem.NA After.rem.Mend.inc.
## 1 m1 a/G 65/175 0.1952259 53 53 53
## 2 m2 a/B 4/272 0.8628440 53 53 53
## 3 <NA> <NA> <NA> NA 53 53 53
## After.rem.unused.haplos pv.overall haplos haplofreq haplofreq.lower
## 1 53 0.06436654 G-a 0.02845996 0.01077831
## 2 53 0.06436654 a-B 0.22419990 0.15893553
## 3 53 0.06436654 G-B 0.74378651 0.65902097
## haplofreq.upper reference RR.est. RR.lower RR.upper RR.p.value
## 1 0.07395674 - 0.000746149 0.000000 Inf 1.00000000
## 2 0.30457495 - 2.011376709 0.930743 4.329449 0.07208091
## 3 0.81404928 ref 1.000000000 1.000000 1.000000 NA
## RRdd.est. RRdd.lower RRdd.upper RRdd.p.value
## 1 6.419717e-07 0.0000000 Inf 1.00000000
## 2 4.045636e+00 0.8662826 18.74412 0.07208091
## 3 1.000000e+00 1.0000000 1.00000 NA
We can plot the results easily with the plot
function
(provided that the ggplot2 package is
installed on your system!):
plot( hapStrat.res )
To check whether there is any significant interaction between the
environmental exposure and genotypes, we use the gxe
function:
gxe( hapStrat.res )
## gxe.test chisq df pval
## 1 haplo.freq 2.74586596 8 0.9492800
## 2 child 3.45632944 8 0.9025524
## 3 haplo.freq.trend 0.08145869 2 0.9600889
## 4 child.trend 0.29423427 2 0.8631929