Skip to content

Commit

Permalink
styler
Browse files Browse the repository at this point in the history
  • Loading branch information
@mjz1
mjz1 committed Dec 28, 2022
1 parent 174d932 commit ce2cfe7
Show file tree
Hide file tree
Showing 6 changed files with 109 additions and 104 deletions.
2 changes: 1 addition & 1 deletion R/bin_utils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ bin_atac_frags <- function(ArrowFiles,
bpparams = BiocParallel::MulticoreParam(
workers = ncores,
progressbar = TRUE
),
),
overwrite = FALSE,
return_mat = FALSE,
...) {
Expand Down
4 changes: 3 additions & 1 deletion R/calc_allelic_imbalance.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@
#'
#' @export
#'
calc_allelic <- function(snp, ncores = 1,
calc_allelic <- function(snp,
ncores = 1,
bins = NULL,
group_var = NULL,
FUN = calc_ai,
Expand Down Expand Up @@ -109,6 +110,7 @@ calc_ai <- function(ref_counts, alt_counts) {
x <- pmax(ref_counts, alt_counts)
y <- pmin(ref_counts, alt_counts)
ai <- (x - y) / x
# ai <- y / (x + y) # Deviation from 0.5

return(ai)
}
148 changes: 74 additions & 74 deletions R/cnv_tools.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,135 +259,135 @@ calc_ratios <- function(sce, assay_name, fun = c("mean", "median"), new_assay =

#' @export
#' @importFrom stats ansari.test wilcox.test
mergeLevels <- function (vecObs, vecPred, pv.thres = 1e-04, ansari.sign = 0.05,
thresMin = 0.05, thresMax = 0.5, verbose = 1, scale = TRUE)
{
mergeLevels <- function(vecObs, vecPred, pv.thres = 1e-04, ansari.sign = 0.05,
thresMin = 0.05, thresMax = 0.5, verbose = 1, scale = TRUE) {
if (thresMin > thresMax) {
cat("Error, thresMax should be equal to or larger than thresMin\n")
return()
}
thresAbs = thresMin
thresAbs <- thresMin
sq <- numeric()
j = 0
ansari = numeric()
lv = numeric()
flag = 0
j <- 0
ansari <- numeric()
lv <- numeric()
flag <- 0
if (thresMin == thresMax) {
flag = 2
}
else {
l.step <- signif((thresMax - thresMin)/10, 1)
s.step <- signif((thresMax - thresMin)/200, 1)
flag <- 2
} else {
l.step <- signif((thresMax - thresMin) / 10, 1)
s.step <- signif((thresMax - thresMin) / 200, 1)
}
while (1) {
if (verbose >= 1) {
cat("\nCurrent thresAbs: ", thresAbs, "\n")
}
j = j + 1
j <- j + 1
sq[j] <- thresAbs
vecPredNow = vecPred
mnNow = unique(vecPred)
mnNow = mnNow[!is.na(mnNow)]
cont = 0
vecPredNow <- vecPred
mnNow <- unique(vecPred)
mnNow <- mnNow[!is.na(mnNow)]
cont <- 0
while (cont == 0 & length(mnNow) > 1) {
mnNow = sort(mnNow)
mnNow <- sort(mnNow)
n <- length(mnNow)
if (verbose >= 2) {
cat("\r", n, ":", length(unique(vecPred)), "\t")
}
if (scale) {
d <- (2 * 2^mnNow)[-n] - (2 * 2^mnNow)[-1]
}
else {
} else {
d <- (mnNow)[-n] - (mnNow)[-1]
}
dst <- cbind(abs(d)[order(abs(d))], (2:n)[order(abs(d))],
(1:(n - 1))[order(abs(d))])
dst <- cbind(
abs(d)[order(abs(d))], (2:n)[order(abs(d))],
(1:(n - 1))[order(abs(d))]
)
for (i in 1:nrow(dst)) {
cont = 1
out = combine.func(diff = dst[i, 1], vecObs,
vecPredNow, mnNow, mn1 = mnNow[dst[i, 2]],
mn2 = mnNow[dst[i, 3]], pv.thres = pv.thres,
thresAbs = if (scale) {
2 * 2^thresAbs - 2
}
else {
thresAbs
})
cont <- 1
out <- combine.func(
diff = dst[i, 1], vecObs,
vecPredNow, mnNow, mn1 = mnNow[dst[i, 2]],
mn2 = mnNow[dst[i, 3]], pv.thres = pv.thres,
thresAbs = if (scale) {
2 * 2^thresAbs - 2
} else {
thresAbs
}
)
if (out$pv > pv.thres) {
cont = 0
vecPredNow = out$vecPredNow
mnNow = out$mnNow
cont <- 0
vecPredNow <- out$vecPredNow
mnNow <- out$mnNow
break
}
}
}
ansari[j] = ansari.test(sort(vecObs - vecPredNow), sort(vecObs -
vecPred))$p.value
ansari[j] <- ansari.test(sort(vecObs - vecPredNow), sort(vecObs -
vecPred))$p.value
if (is.na(ansari[j])) {
ansari[j] = 0
ansari[j] <- 0
}
lv[j] = length(mnNow)
lv[j] <- length(mnNow)
if (flag == 2) {
break
}
if (ansari[j] < ansari.sign) {
flag = 1
flag <- 1
}
if (flag) {
if (ansari[j] > ansari.sign | thresAbs == thresMin) {
break
}
else {
thresAbs = signif(thresAbs - s.step, 3)
} else {
thresAbs <- signif(thresAbs - s.step, 3)
if (thresAbs <= thresMin) {
thresAbs = thresMin
thresAbs <- thresMin
}
}
}
else {
thresAbs = thresAbs + l.step
} else {
thresAbs <- thresAbs + l.step
}
if (thresAbs >= thresMax) {
thresAbs = thresMax
flag = 2
thresAbs <- thresMax
flag <- 2
}
}
return(list(vecMerged = vecPredNow, mnNow = mnNow, sq = sq,
ansari = ansari))
return(list(
vecMerged = vecPredNow, mnNow = mnNow, sq = sq,
ansari = ansari
))
}


#' @export
combine.func <- function (diff, vecObs, vecPredNow, mnNow, mn1, mn2, pv.thres = 1e-04,
thresAbs = 0)
{
vec1 = vecObs[which(vecPredNow == mn1)]
vec2 = vecObs[which(vecPredNow == mn2)]
combine.func <- function(diff, vecObs, vecPredNow, mnNow, mn1, mn2, pv.thres = 1e-04,
thresAbs = 0) {
vec1 <- vecObs[which(vecPredNow == mn1)]
vec2 <- vecObs[which(vecPredNow == mn2)]
if (diff <= thresAbs) {
pv = 1
}
else {
if ((length(vec1) > 10 & length(vec2) > 10) | sum(length(vec1),
length(vec2)) > 100) {
pv = wilcox.test(vec1, vec2)$p.value
}
else {
pv = wilcox.test(vec1, vec2, exact = TRUE)$p.value
pv <- 1
} else {
if ((length(vec1) > 10 & length(vec2) > 10) | sum(
length(vec1),
length(vec2)
) > 100) {
pv <- wilcox.test(vec1, vec2)$p.value
} else {
pv <- wilcox.test(vec1, vec2, exact = TRUE)$p.value
}
if (length(vec1) <= 3 | length(vec2) <= 3) {
pv = 0
pv <- 0
}
}
index.merged <- numeric()
if (pv > pv.thres) {
vec = c(vec1, vec2)
index.merged = which((vecPredNow == mn1) | (vecPredNow ==
mn2))
vecPredNow[index.merged] = median(vec, na.rm = TRUE)
mnNow[which((mnNow == mn1) | (mnNow == mn2))] = median(vec,
na.rm = TRUE)
mnNow = unique(mnNow)
vec <- c(vec1, vec2)
index.merged <- which((vecPredNow == mn1) | (vecPredNow ==
mn2))
vecPredNow[index.merged] <- median(vec, na.rm = TRUE)
mnNow[which((mnNow == mn1) | (mnNow == mn2))] <- median(vec,
na.rm = TRUE
)
mnNow <- unique(mnNow)
}
list(mnNow = mnNow, vecPredNow = vecPredNow, pv = pv)
}
Expand Down
19 changes: 9 additions & 10 deletions R/run_scatools.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ run_scatools <- function(sample_id,
verbose = TRUE,
save_h5ad = TRUE,
ncores = 1) {

# TODO INPUT VALIDATION

dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
Expand Down Expand Up @@ -59,12 +58,14 @@ run_scatools <- function(sample_id,
ArrowFile <- glue::glue("{archr_dir}/{sample_id}.arrow")
}

proj <- create_arrow_file(archr_dir = archr_dir,
ArrowFile = ArrowFile,
fragments = fragments,
sample_id = sample_id,
ncores = ncores,
force = force_arrow)
proj <- create_arrow_file(
archr_dir = archr_dir,
ArrowFile = ArrowFile,
fragments = fragments,
sample_id = sample_id,
ncores = ncores,
force = force_arrow
)

# Binning and CNV processing
dir.create(bins_out, recursive = TRUE, showWarnings = FALSE)
Expand Down Expand Up @@ -99,7 +100,7 @@ run_scatools <- function(sample_id,
save_to(object = sce, save_to = file.path(processed, glue::glue("{bin_name}_raw.sce")))

# Remove any small leftover bins less than 10% of the full bin length
sce_processed <- sce[rowRanges(sce)$binwidth > 0.1*getmode(width(bins)), ]
sce_processed <- sce[rowRanges(sce)$binwidth > 0.1 * getmode(width(bins)), ]

# Length normalize so the tail end bins are in the same scale
sce_processed <- length_normalize(sce_processed, assay_name = "counts", assay_to = "counts")
Expand Down Expand Up @@ -136,7 +137,6 @@ run_scatools <- function(sample_id,


create_arrow_file <- function(fragments = NULL, ArrowFile, archr_dir, sample_id = "sample", genome = "hg38", ncores = 1, force = FALSE) {

if (!require("ArchR", quietly = TRUE)) {
logger::log_error("'ArchR' package is required to create ArrowFiles from fragments files")
stop()
Expand Down Expand Up @@ -204,4 +204,3 @@ create_arrow_file <- function(fragments = NULL, ArrowFile, archr_dir, sample_id

return(proj)
}

20 changes: 11 additions & 9 deletions data-raw/test_sce.R
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
## code to prepare `test_sce` dataset goes here
library(scatools)

ncores = parallelly::availableCores()
ncores <- parallelly::availableCores()

fragment_file <- system.file("extdata", "fragments.bed.gz", package = "scatools")
names(fragment_file) <- "test_sample"
Expand All @@ -13,13 +13,15 @@ bin_name <- prettyMb(getmode(width(bins)))
outdir <- tempdir()


test_sce <- run_scatools(sample_id = names(fragment_file),
fragments = fragment_file,
outdir = outdir,
verbose = TRUE,
overwrite = TRUE,
force_arrow = TRUE,
ncores = ncores,
bins = bins)
test_sce <- run_scatools(
sample_id = names(fragment_file),
fragments = fragment_file,
outdir = outdir,
verbose = TRUE,
overwrite = TRUE,
force_arrow = TRUE,
ncores = ncores,
bins = bins
)

usethis::use_data(test_sce, overwrite = TRUE)
20 changes: 11 additions & 9 deletions vignettes/scatools.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,14 +62,16 @@ bins <- get_tiled_bins(bs_genome = BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapien

A convenience wrapper is provided to perform the analysis
```{r, eval=FALSE}
sce <- run_scatools(sample_id = names(fragment_file),
fragments = fragment_file,
outdir = "./example/",
verbose = TRUE,
overwrite = FALSE,
force_arrow = FALSE,
ncores = ncores,
bins = bins)
sce <- run_scatools(
sample_id = names(fragment_file),
fragments = fragment_file,
outdir = "./example/",
verbose = TRUE,
overwrite = FALSE,
force_arrow = FALSE,
ncores = ncores,
bins = bins
)
```


Expand All @@ -96,7 +98,7 @@ We can visualize the results as a heatmap.
```{r cnaheatmap_plot, message=FALSE}
library(ComplexHeatmap)
sce <- sce[,order(sce$tumor_cell, sce$clusters)]
sce <- sce[, order(sce$tumor_cell, sce$clusters)]
col_clones <- dittoColors()
col_clones <- col_clones[1:length(unique(sce[["clusters"]]))]
Expand Down

0 comments on commit ce2cfe7

Please sign in to comment.