## ----include=FALSE, cache=FALSE--------------------------------------------
library(RVS)
library(kinship2)

## --------------------------------------------------------------------------
data(samplePedigrees) # load sample pedigrees
kinship2::plot.pedigree(samplePedigrees$firstCousinPair) # plot pedigree

## --------------------------------------------------------------------------
RVsharing(samplePedigrees$firstCousinPair)

## --------------------------------------------------------------------------
defaultProbs <- sapply(samplePedigrees[1:4], function(p) 
suppressMessages(RVsharing(p)))

exactProbs <- list()
freq <- c(0.001,0.0025,0.005,0.01,0.025,0.05)
exactProbs$fistCousinPair <- sapply(freq, function(f) suppressMessages(
    RVsharing(samplePedigrees$firstCousinPair, alleleFreq=f)))
exactProbs$secondCousinPair <- sapply(freq, function(f) suppressMessages(
    RVsharing(samplePedigrees$secondCousinPair, alleleFreq=f)))
exactProbs$firstCousinTriple <- sapply(freq, function(f) suppressMessages(
    RVsharing(samplePedigrees$firstCousinTriple, alleleFreq=f)))
exactProbs$secondCousinTriple <- sapply(freq, function(f) suppressMessages(
    RVsharing(samplePedigrees$secondCousinTriple, alleleFreq=f)))

plot(NULL, xlim=c(0.001,0.05), ylim=c(0,0.12),log='x',xaxt='n', 
ylab='probability of sharing', xlab='variant frequency [%]')
axis(side=1, at=freq, labels=freq*100)
invisible(sapply(exactProbs, lines, x=freq))
invisible(sapply(exactProbs, points, x=freq, pch=19))
abline(h=defaultProbs,lty=2)

## --------------------------------------------------------------------------
ped <- samplePedigrees$firstCousinTriple
ped$affected[9] <- 0
plot(ped)

p <- RVsharing(ped)
p <- RVsharing(ped, useAffected=TRUE)
p <- RVsharing(ped, carriers=c(9,10))
p <- RVsharing(ped, carriers=c(10,11), useAffected=TRUE)

## --------------------------------------------------------------------------
p <- RVsharing(samplePedigrees$firstCousinPair)
p <- RVsharing(samplePedigrees$firstCousinPair, kinshipCoeff = 0.05)
kCoeff <- seq(0.01,0.2,0.01)
sharingProb <- sapply(kCoeff, function(c) suppressMessages(
    RVsharing(samplePedigrees$firstCousinPair, kinshipCoeff=c)))
plot(kCoeff, sharingProb, type='l')

## --------------------------------------------------------------------------
plot(samplePedigrees$twoGenerationsInbreeding)
ComputeKinshipPropCoef(samplePedigrees$twoGenerationsInbreeding)

## --------------------------------------------------------------------------
p <- RVsharing(samplePedigrees$firstCousinPair)
p <- RVsharing(samplePedigrees$firstCousinPair, nSim=1e4)

p <- RVsharing(samplePedigrees$firstCousinPair, alleleFreq=0.01)
p <- RVsharing(samplePedigrees$firstCousinPair, alleleFreq=0.01, nSim=1e4)

p <- RVsharing(samplePedigrees$firstCousinPair, kinshipCoeff=0.05)
p <- RVsharing(samplePedigrees$firstCousinPair, kinshipCoeff=0.05, nSim=1e4)

## --------------------------------------------------------------------------
# assumption that 1 founder introduces
fDist <- function(N) sample(c(rep(0,N-1), 1))
RVsharing(samplePedigrees$firstCousinPair, nSim=1e4, founderDist=fDist)
RVsharing(samplePedigrees$firstCousinPair)

## --------------------------------------------------------------------------
probs <- sapply(samplePedigrees, function(p) suppressMessages(RVsharing(p)))
observed <- c(rep(FALSE, 3), rep(TRUE, 2), rep(FALSE, 4))
multipleFamilyPValue(probs, observed)

## --------------------------------------------------------------------------
carriers = c(15,16,17)
carrier.sets = list()
for (i in length(carriers):1)
carrier.sets = c(carrier.sets, combn(carriers,i,simplify=FALSE))
fam15157.pattern.prob = sapply(carrier.sets,function (vec) RVsharing(samplePedigrees$secondCousinTriple,carriers=vec))
fam15157.N = sapply(carrier.sets,length)

## --------------------------------------------------------------------------
fam15157.pattern.prob = c(RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15,16,17)),
    RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15,16)),
    RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15)))
fam15157.N = 3:1
fam15157.nequiv = c(1,3,3)

## --------------------------------------------------------------------------
sum(fam15157.pattern.prob*fam15157.nequiv)

## --------------------------------------------------------------------------
fam15157 <- samplePedigrees$secondCousinTriple
fam15157$affected[3] = 1
plot(fam15157)

## --------------------------------------------------------------------------
carriers = c(3,15,16,17)
carrier.sets = list()
for (i in length(carriers):1)
carrier.sets = c(carrier.sets, combn(carriers,i,simplify=FALSE))
carrier.sets
carrier.sets = carrier.sets[-c(5,9,10)]
fam15157.pattern.prob = sapply(carrier.sets,function (vec) RVsharing(fam15157,carriers=vec,useAffected=TRUE))
fam15157.N = sapply(carrier.sets,length)

## --------------------------------------------------------------------------
sum(fam15157.pattern.prob)

## --------------------------------------------------------------------------
pobs = RVsharing(fam15157,carriers=c(3,16,17),useAffected=TRUE)

## --------------------------------------------------------------------------
sum(fam15157.pattern.prob[fam15157.pattern.prob<=pobs & fam15157.N >= 3])

## --------------------------------------------------------------------------
fam15157.pattern.prob = c(RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15,16,17)),
    RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15,16)),
    RVsharing(samplePedigrees$secondCousinTriple,carriers=c(15)))
fam15157.N = 3:1
fam15157.nequiv = c(1,3,3)

fam28003.pattern.prob = c(RVsharing(samplePedigrees$firstAndSecondCousinsTriple,carriers=c(36,104,110)),
RVsharing(samplePedigrees$firstAndSecondCousinsTriple,carriers=c(36,104)),
RVsharing(samplePedigrees$firstAndSecondCousinsTriple,carriers=c(104,110)),
RVsharing(samplePedigrees$firstAndSecondCousinsTriple,carriers=c(36)),
RVsharing(samplePedigrees$firstAndSecondCousinsTriple,carriers=c(104)))
fam28003.N = c(3,2,2,1,1)
fam28003.nequiv = c(1,2,1,1,2)

ex.pattern.prob.list = list("15157"=fam15157.pattern.prob,"28003"=fam28003.pattern.prob)
ex.nequiv.list = list("15157"=fam15157.nequiv,"28003"=fam28003.nequiv)
ex.N.list = list("15157"=fam15157.N,"28003"=fam28003.N)

## --------------------------------------------------------------------------
data(famVCF)
library(VariantAnnotation)
fam15157.snp = genotypeToSnpMatrix(fam15157.vcf)
fam28003.snp = genotypeToSnpMatrix(fam28003.vcf)

## --------------------------------------------------------------------------
ex.SnpMatrix.list = list(fam15157=fam15157.snp$genotypes,fam28003=fam28003.snp$genotypes)
ex.ped.obj = list(fam15157=samplePedigrees$secondCousinTriple,fam28003=samplePedigrees$firstAndSecondCousinsTriple)

## --------------------------------------------------------------------------
sites = c(92,119)
ex.SnpMatrix.list[["fam15157"]][,sites[1]]@.Data
ex.SnpMatrix.list[["fam28003"]][,sites[2]]@.Data

## --------------------------------------------------------------------------
RVgene(ex.SnpMatrix.list,ex.ped.obj,sites,pattern.prob.list=ex.pattern.prob.list,
nequiv.list=ex.nequiv.list,N.list=ex.N.list,type="count")