A Monte Carlo permutation method for multiple test correction.

Share:

Description

Permutation tests exist for any test statistic, regardless of whether or not its distribution is known. Thus the permutation test is widely considered the gold standard for accurate multiple testing correction.

For example,for case/control association study for SNPs, the permutation test proceeds as follows:

1) Combine the observations from all the samples;

2) Shuffle them and rearrangements of the labels(case/control) on the observed data;

3) Record the genotype frequency of case and control samples, respectively;

4) Calculate the statistic of interest;

5) Repeat many times(at least 1000) to obtain the distribution of the statistic;

6) Determine how often the resampled statistic of interest is as extreme as the observed value of the same statistic.

Obviously, for multiple test correction in case/control association study for millions of SNPs, the traditional method—permutation test is very computationally impractical. Thus propose an accurate, rapid and efficient method for multiple testing correction in genome-wide association studies—MCPerm.

Method—MCPerm generates the genotype frequency for rearranged case and control data by twice generating random numbers for the hypergeometric distribution, based on the genotype statistic of original data, taking the place of the step 2) and step 3) of the traditional method. And the genotype frequency distribution generating by MCPerm is almost the same with permutation test, this simplified method greatly improves the efficiency of the permutation test and is faster. MCPerm method can be the perfect alternative to permutation test.

Details

Package: MCPerm
Type: Package
Version: 1.1.4
Date: 2013-06-12
License: GPL-2

Author(s)

Lanying Zhang and Yongshuai Jiang <jiangyongshuai@gmail.com>

References

William S Noble(Nat Biotechnol.2009): How does mutiple testing correction work?

Edgington. E.S.(1995): Randomization tests, 3rd ed.

Julian P.T.Higgins, Simon G.Thompson(Statistics in Medicine,2002): Quantifying heterogeneity in a meta-analysis.

See Also

Armitage, Armitage.TradPerm, Armitage.MCPerm, OR, OR.TradPerm, OR.MCPerm, permuteGenotype, permuteGenotypeCount, genotypeStat, chisq.TradPerm, chisq.MCPerm, fisher.TradPerm, fisher.MCPerm, rhyper, chisq.test, fisher.test, meta, meta.TradPerm, meta.MCPerm, VS.Genotype.Hist, VS.Allele.Hist, VS.Hist, PermMeta.LnOR.Hist, PermMeta.LnOR.CDC, PermMeta.LnOR.boxplot, PermMeta.boxplot, PermMeta.Hist, pearson_scatter, Q.TradPerm, Q.MCPerm, I2.TradPerm, I2.MCPerm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
## example1-----genotypeStat-----------
##  import example data
# data(genotypeData)
## get the first line: affection state for samples
# data1=genotypeData[1,,drop=FALSE]
## get the second line: genotype data for a sepcifed snp
# data2=genotypeData[2,,drop=FALSE]
## Statistical allele and genotype frequency of the specified snp for case-control data.
# result=genotypeStat(data2,data1,fromCol=2,naString="?_?",sep="_")
# genotypeCount=result$genotypeCount
# alleleCount=result$alleleCount

## example2-----permuteGenotype-----------
## a matrix with 1 row
# dataLine=matrix(c("rs12","1","2","3","4","5"),nrow=1)
## permutate the elements of the matrix
# newData=permuteGenotype(dataLine=dataLine,fromCol=2)

## example3-----permuteGenotypeCount-----------
# newMatrix=permuteGenotypeCount(case_11=24,case_12=34,case_22=45,
#   control_11=23,control_12=45,control_22=34,n=5)

## example4-----OR-----------
## OR(odd ratio) for the risk-allele
# ORvalue=OR(case_allele1=20,case_allele2=30,control_allele1=10,control_allele2=60)

## example5-----OR.TradPerm--------
## import example data
# data(genotypeData)
## get the first line: affection state for samples
# data1=genotypeData[1,,drop=FALSE]
## get the second line: genotype data for a sepcifed snp
# data2=genotypeData[2,,drop=FALSE]
# result1=OR.TradPerm(genotypeLine=data2,affectionLine=data1,fromCol=2,naString="?_?",
  # sep="_",repeatNum=5)
# risk_allele=result1$risk_allele
# p=result1$pValue
# obsOR=result$OR

## example6-----OR.MCPerm----------
# OR.MCPerm(case_allele1=34,case_allele2=23,control_allele1=27,control_allele2=45,repeatNum=5)

## example7-----Armitage-------
# Armitage(case_11=23,case_12=45,case_22=12,control_11=27,control_12=12,control_22=45)

## example8-----Armitage.TradPerm-----
## import example data
# data(genotypeData)
## get the first line: affection state for samples
# data1=genotypeData[1,,drop=FALSE]
## get the second line: genotype data for a sepcifed snp
# data2=genotypeData[2,,drop=FALSE]
# Armitage.TradPerm(genotypeLine=data2,affectionLine=data1,
   # fromCol=2,naString="?_?",sep="_",repeatNum=1000)

## example9----Armitage.MCPerm--------
# Armitage.MCPerm(case_11=23,case_12=45,case_22=12,
   # control_11=27,control_12=12,control_22=45,repeatNum=1000)

## example10---chisq.TradPerm---------
## import example data
data(genotypeData)
## get the first line: affection state for samples
# data1=genotypeData[1,,drop=FALSE]
## get the second line: genotype data for a sepcifed snp
# data2=genotypeData[2,,drop=FALSE]
# chisq.TradPerm(genotypeLine=data2,affectionLine=data1,
   # fromCol=2,naString="?_?",sep="_",repeatNum=1000)

## example11---chisq.MCPerm--------
# case_11=23
# case_12=45
# case_22=12
# control_11=27
# control_12=12
# control_22=45
# chisq.MCPerm(23,45,12,27,12,45,repeatNum=5)

## example12---fisher.TradPerm------
# import example data
# data(genotypeData)
## get the first line: affection state for samples
# data1=genotypeData[1,,drop=FALSE]
## get the second line: genotype data for a sepcifed snp
# data2=genotypeData[2,,drop=FALSE]
# fisher.TradPerm(genotypeLine=data2,affectionLine=data1,
   # fromCol=2,naString="?_?",sep="_",repeatNum=5)

## example13---fisher.MCPerm-------
#fisher.MCPerm(23,45,12,27,12,45,repeatNum=5)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.