## ----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()