## ---- eval=FALSE---------------------------------------------------------
#  geno <- MatrixGenotypeReader(genotype = genotype, snpID = snpID, chromosome = chromosome,
#                               position = position, scanID = scanID)
#  genoData <- GenotypeData(geno)

## ---- eval=FALSE---------------------------------------------------------
#  geno <- GdsGenotypeReader(filename = "genotype.gds")
#  genoData <- GenotypeData(geno)

## ---- eval=FALSE---------------------------------------------------------
#  snpgdsBED2GDS(bed.fn = "genotype.bed", bim.fn = "genotype.bim", fam.fn = "genotype.fam",
#                out.gdsfn = "genotype.gds")

## ---- echo=FALSE, message=FALSE------------------------------------------
library(GENESIS)
library(GWASTools)

## ------------------------------------------------------------------------
# read in GDS data
gdsfile <- system.file("extdata", "HapMap_ASW_MXL_geno.gds", package="GENESIS")
HapMap_geno <- GdsGenotypeReader(filename = gdsfile)
# create a GenotypeData class object
HapMap_genoData <- GenotypeData(HapMap_geno)
HapMap_genoData

## ------------------------------------------------------------------------
# read individual IDs from GenotypeData object
iids <- getScanID(HapMap_genoData)
head(iids)
# create matrix of KING estimates
KINGmat <- king2mat(file.kin0 = system.file("extdata", "MXL_ASW.kin0", package="GENESIS"), 
    				file.kin = system.file("extdata", "MXL_ASW.kin", package="GENESIS"), 
					iids = iids)
KINGmat[1:5,1:5]

## ------------------------------------------------------------------------
# run PC-AiR
mypcair <- pcair(genoData = HapMap_genoData, kinMat = KINGmat, divMat = KINGmat)

## ---- eval=FALSE---------------------------------------------------------
#  mypcair <- pcair(genoData = HapMap_genoData, unrel.set = IDs)

## ---- eval=FALSE---------------------------------------------------------
#  mypcair <- pcair(genoData = HapMap_genoData, kinMat = KINGmat, divMat = KINGmat, unrel.set = IDs)

## ------------------------------------------------------------------------
part <- pcairPartition(kinMat = KINGmat, divMat = KINGmat)

## ------------------------------------------------------------------------
head(part$unrels); head(part$rels)

## ------------------------------------------------------------------------
summary(mypcair)

## ---- fig.show='hold', fig.width=3.4, fig.height=3.4, dev.args=list(pointsize = 10, bg='white')----
# plot top 2 PCs
plot(mypcair)
# plot PCs 3 and 4
plot(mypcair, vx = 3, vy = 4)

## ------------------------------------------------------------------------
# run PC-Relate
mypcrelate <- pcrelate(genoData = HapMap_genoData, pcMat = mypcair$vectors[,1:2], 
                       training.set = mypcair$unrels)

## ---- eval=FALSE---------------------------------------------------------
#  library(gdsfmt)
#  mypcrelate <- openfn.gds("tmp_pcrelate.gds")

## ------------------------------------------------------------------------
pcrelateReadKinship(pcrelObj = mypcrelate, scan.include = iids[1:40], kin.thresh = 2^(-9/2))
pcrelateReadInbreed(pcrelObj = mypcrelate, scan.include = iids[1:40], f.thresh = 2^(-11/2))
pcrelateMakeGRM(pcrelObj = mypcrelate, scan.include = iids[1:5], scaleKin = 2)

## ---- echo=FALSE---------------------------------------------------------
close(HapMap_genoData)