Description Usage Arguments Details Value Author(s) References See Also Examples
Calculate the log likelihood values from maximum likelihood estimation.
1 2 3 | snpgdsPairIBDMLELogLik(geno1, geno2, allele.freq, k0=NaN, k1=NaN,
relatedness=c("", "self", "fullsib", "offspring", "halfsib",
"cousin", "unrelated"), verbose=TRUE)
|
geno1 |
the SNP genotypes for the first individual, 0 – BB, 1 – AB, 2 – AA, other values – missing |
geno2 |
the SNP genotypes for the second individual, 0 – BB, 1 – AB, 2 – AA, other values – missing |
allele.freq |
the allele frequencies |
k0 |
specified IBD coefficient |
k1 |
specified IBD coefficient |
relatedness |
specify a relatedness, otherwise use the values of k0 and k1 |
verbose |
if TRUE, show information |
If (relatedness
== "") and (k0 == NaN or k1 == NaN), then return
the log likelihood values for each (k0, k1) stored in ibdobj.
If (relatedness
== "") and (k0 != NaN) and (k1 != NaN), then return
the log likelihood values for a specific IBD coefficient (k0, k1).
If relatedness
is: "self", then k0 = 0, k1 = 0; "fullsib", then
k0 = 0.25, k1 = 0.5; "offspring", then k0 = 0, k1 = 1; "halfsib", then
k0 = 0.5, k1 = 0.5; "cousin", then k0 = 0.75, k1 = 0.25; "unrelated", then
k0 = 1, k1 = 0.
The value of log likelihood.
Xiuwen Zheng
Milligan BG. 2003. Maximum-likelihood estimation of relatedness. Genetics 163:1153-1167.
Weir BS, Anderson AD, Hepler AB. 2006. Genetic relatedness analysis: modern data and new challenges. Nat Rev Genet. 7(10):771-80.
Choi Y, Wijsman EM, Weir BS. 2009. Case-control association testing in the presence of unknown relationships. Genet Epidemiol 33(8):668-78.
snpgdsPairIBD
, snpgdsIBDMLE
,
snpgdsIBDMLELogLik
, snpgdsIBDMoM
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 | # open an example dataset (HapMap)
genofile <- snpgdsOpen(snpgdsExampleFileName())
YRI.id <- read.gdsn(index.gdsn(genofile, "sample.id"))[
read.gdsn(index.gdsn(genofile, "sample.annot/pop.group"))=="YRI"]
# SNP pruning
set.seed(10)
snpset <- snpgdsLDpruning(genofile, sample.id=YRI.id, maf=0.05,
missing.rate=0.05)
snpset <- unname(sample(unlist(snpset), 250))
# the number of samples
n <- 25
# specify allele frequencies
RF <- snpgdsSNPRateFreq(genofile, sample.id=YRI.id, snp.id=snpset,
with.id=TRUE)
summary(RF$AlleleFreq)
subMLE <- snpgdsIBDMLE(genofile, sample.id=YRI.id[1:n], snp.id=RF$snp.id,
allele.freq=RF$AlleleFreq)
subMoM <- snpgdsIBDMoM(genofile, sample.id=YRI.id[1:n], snp.id=RF$snp.id,
allele.freq=RF$AlleleFreq)
# genotype matrix
mat <- snpgdsGetGeno(genofile, sample.id=YRI.id[1:n], snp.id=snpset,
snpfirstdim=TRUE)
########################
rv <- NULL
for (i in 2:n)
{
rv <- rbind(rv, snpgdsPairIBD(mat[,1], mat[,i], RF$AlleleFreq, "EM"))
print(snpgdsPairIBDMLELogLik(mat[,1], mat[,i], RF$AlleleFreq,
relatedness="unrelated", verbose=TRUE))
}
rv
summary(rv$k0 - subMLE$k0[1, 2:n])
summary(rv$k1 - subMLE$k1[1, 2:n])
# ZERO
rv <- NULL
for (i in 2:n)
rv <- rbind(rv, snpgdsPairIBD(mat[,1], mat[,i], RF$AlleleFreq, "MoM"))
rv
summary(rv$k0 - subMoM$k0[1, 2:n])
summary(rv$k1 - subMoM$k1[1, 2:n])
# ZERO
# close the genotype file
snpgdsClose(genofile)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.