## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
suppressPackageStartupMessages({
  library(comapr)
  library(SummarizedExperiment)
  })

## -----------------------------------------------------------------------------
demo_path <-paste0(system.file("extdata",package = "comapr"),"/")

## -----------------------------------------------------------------------------
list.files(demo_path)

## -----------------------------------------------------------------------------
pcQC <- perCellChrQC(sampleName="s1",chroms=c("chr1"),
                     path=paste0(demo_path,"/"),
                     barcodeFile=NULL)

## -----------------------------------------------------------------------------
s1_rse_state <- readHapState("s1",chroms=c("chr1","chr2"),
                             path=demo_path,barcodeFile=NULL,
                             minSNP = 0,
                             minlogllRatio = 50,
                             bpDist = 100,
                             maxRawCO=10,
                             minCellSNP = 1)

s2_rse_state <- readHapState("s2",chroms=c("chr1","chr2"),
                             path=demo_path,
                             barcodeFile=NULL,
                             minSNP = 0,
                             minlogllRatio = 50,
                             bpDist = 100,
                             maxRawCO=10,
                             minCellSNP = 1)

## -----------------------------------------------------------------------------
s1_rse_state

## -----------------------------------------------------------------------------
assay(s1_rse_state)

## -----------------------------------------------------------------------------
rowRanges(s1_rse_state)

## -----------------------------------------------------------------------------
colData(s1_rse_state)

colData(s1_rse_state)$sampleGroup <- "s1"

colData(s2_rse_state)$sampleGroup <- "s2"


## -----------------------------------------------------------------------------
twoSample <- combineHapState(s1_rse_state,
                             s2_rse_state)

## -----------------------------------------------------------------------------
twoSample <- combineHapState(s1_rse_state,s2_rse_state)

## -----------------------------------------------------------------------------

assay(twoSample)

## -----------------------------------------------------------------------------
cocounts <- countCOs(twoSample)

## -----------------------------------------------------------------------------
rowRanges(cocounts)

## -----------------------------------------------------------------------------
colData(cocounts)

## -----------------------------------------------------------------------------
assay(cocounts)

## -----------------------------------------------------------------------------
dist_gr <- calGeneticDist(co_count = cocounts, mapping_fun = "k")
dist_gr

## -----------------------------------------------------------------------------
rowData(dist_gr)

## -----------------------------------------------------------------------------
## sampleGroup is a column in the colData slot
dist_gr <- calGeneticDist(co_count = cocounts,
                          group_by  = "sampleGroup",
                          mapping_fun = "k")


## -----------------------------------------------------------------------------
rowData(dist_gr)$kosambi

## -----------------------------------------------------------------------------
p_gr <- granges(dist_gr)
mcols(p_gr) <- rowData(dist_gr)$kosambi



## -----------------------------------------------------------------------------
plotWholeGenome(p_gr)


## -----------------------------------------------------------------------------
plotGeneticDist(p_gr,chr = "chr1",cumulative = TRUE)


## -----------------------------------------------------------------------------
set.seed(100)
bootsDiff <- bootstrapDist(co_gr = cocounts,B=10,
                           group_by = "sampleGroup")

## -----------------------------------------------------------------------------
hist(bootsDiff)

## -----------------------------------------------------------------------------
quantile(bootsDiff,c(0.025,0.975),na.rm =TRUE)

## -----------------------------------------------------------------------------
set.seed(100)
perms <- permuteDist(cocounts,B=1000,group_by = "sampleGroup")

## -----------------------------------------------------------------------------
statmod::permp(x = sum(perms$permutes>=perms$observed_diff,
                       na.rm = TRUE),
               nperm = sum(!is.na(perms$permutes)),
               n1 = perms$nSample[1],
               n2 = perms$nSample[2])


## -----------------------------------------------------------------------------
sessionInfo()