Hi, I am trying to integrate data with 3 conditions. Here is my code that I used for all three conditions to create the Seurat objects: I quantified all the three conditions using custom peak set, such that the features are consistent across,
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
set.seed(1234)
regulome_peaks = read.table("scATAC/regulome_clustering_V7.txt")
regulome_peaks = StringToGRanges(regulome_peaks$V1, sep = c(":", "-"))
Cond1_metadata <- read.csv(
file = "scATAC/Cond1/outs/singlecell.csv", eader = TRUE, ow.names = 1
)
Cond1_fragment.path <- 'scATAC/Cond1/outs/fragments.tsv.gz'
Cond1_fragment <- CreateFragmentObject(
path = Cond1_fragment.path, alidate.fragments = FALSE
)
Cond1_drug_FM <- FeatureMatrix(
fragments = Cond1_fragment,features = regulome_peaks,
)
rowcounts <- rowSums(Cond1_drug_FM)
Cond1_drug_FM <- Cond1_drug_FM[rowcounts>30,] #To eliminate features with less than 30 counts
Cond1_drug_assay <- CreateChromatinAssay(
counts = Cond1_drug_FM, ep = c("-", "-"), enome = 'hg19', in.features = 5000, in.cells = 0
)
Cond1_drug <- CreateSeuratObject(
counts = Cond1_drug_assay, ssay = "peaks",meta.data = Cond1_metadata
)
Fragments(Cond1_drug) <- CreateFragmentObject(
path = Cond1_fragment.path, cells = colnames(Cond1_drug), tolerance = 0.5
)
Cond1_drug <- NucleosomeSignal(object = Cond1_drug)
Cond1_drug$nucleosome_group <- ifelse(Cond1_drug$nucleosome_signal > 10, 'NS > 10', 'NS < 10')
Cond1_drug$pct_reads_in_peaks <- Cond1_drug$peak_region_fragments / Cond1_drug$passed_filters * 100
Cond1_drug$blacklist_ratio <- Cond1_drug$blacklist_region_fragments / Cond1_drug$peak_region_fragments
Cond1_drug <- subset(Cond1_drug, subset = peak_region_fragments > 3000 & peak_region_fragments < 25000 & pct_reads_in_peaks > 25 & blacklist_ratio < 0.05 & nucleosome_signal < 5)
Cond1_drug <- RunTFIDF(Cond1_drug)
Cond1_drug <- FindTopFeatures(Cond1_drug)
Cond1_drug <- RunSVD(object = Cond1_drug,assay = 'peaks', irlba.work=500
)
Cond1_drug <- RunUMAP(object = Cond1_drug, reduction = 'lsi', dims = 2:30)
Cond1_drug <- FindNeighbors(object = Cond1_drug, reduction = 'lsi', dims = 2:30)
Cond1_drug <- FindClusters(object = Cond1_drug, verbose = FALSE, algorithm = 3)
This was done for 3 conditions, which gives beautiful clusters that are cell-type specific. First I subset peaks that are common across all the conditions and then used them to subset my objects and tried integration.
This throws the following error:
Warning message in CheckDuplicateCellNames(object.list = object.list):
"Some cell names are duplicated across objects provided. Renaming to enforce unique cell names."Scaling features for provided objects
Warning message:
"UNRELIABLE VALUE: Future ('future_lapply-1') unexpectedly generated random numbers without specifying argument '[future.]seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify argument '[future.]seed', e.g. 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use [future].seed=NULL, or set option 'future.rng.onMisuse' to "ignore"."Warning message:
"UNRELIABLE VALUE: Future ('future_lapply-2') unexpectedly generated random numbers without specifying argument '[future.]seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify argument '[future.]seed', e.g. 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use [future].seed=NULL, or set option 'future.rng.onMisuse' to "ignore"."Warning message:
"UNRELIABLE VALUE: Future ('future_lapply-3') unexpectedly generated random numbers without specifying argument '[future.]seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify argument '[future.]seed', e.g. 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use [future].seed=NULL, or set option 'future.rng.onMisuse' to "ignore"."Finding anchors between all query and reference datasets
Running CCA
Merging objects
Error in CreateChromatinAssay(data = merged.data, min.cells = 0, min.features = 0, : Length of supplied genomic ranges does not match number
of rows in matrix
Traceback:
1. FindIntegrationAnchors(object.list = list(Cond2_drug_comm,
. Cond1_drug_comm,Cond3_drug_comm), anchor.features = peaks.use_common_all,
. reference = 1)
2. my.lapply(X = 1:nrow(x = combinations), FUN = function(row) {
. i <- combinations[row, 1]
. j <- combinations[row, 2]
. object.1 <- DietSeurat(object = object.list[[i]], assays = assay[i],
. features = anchor.features, counts = FALSE, scale.data = TRUE,
. dimreducs = reduction)
. object.2 <- DietSeurat(object = object.list[[j]], assays = assay[j],
. features = anchor.features, counts = FALSE, scale.data = TRUE,
. dimreducs = reduction)
. suppressWarnings(object.1[["ToIntegrate"]] <- object.1[[assay[i]]])
. DefaultAssay(object = object.1) <- "ToIntegrate"
. if (reduction %in% Reductions(object = object.1)) {
. slot(object = object.1[[reduction]], name = "assay.used") <- "ToIntegrate"
. }
. object.1 <- DietSeurat(object = object.1, assays = "ToIntegrate",
. scale.data = TRUE, dimreducs = reduction)
. suppressWarnings(object.2[["ToIntegrate"]] <- object.2[[assay[j]]])
. DefaultAssay(object = object.2) <- "ToIntegrate"
. if (reduction %in% Reductions(object = object.2)) {
. slot(object = object.2[[reduction]], name = "assay.used") <- "ToIntegrate"
. }
. object.2 <- DietSeurat(object = object.2, assays = "ToIntegrate",
. scale.data = TRUE, dimreducs = reduction)
. object.pair <- switch(EXPR = reduction, cca = {
. object.pair <- RunCCA(object1 = object.1, object2 = object.2,
. assay1 = "ToIntegrate", assay2 = "ToIntegrate", features = anchor.features,
. num.cc = max(dims), renormalize = FALSE, rescale = FALSE,
. verbose = verbose)
. if (l2.norm) {
. object.pair <- L2Dim(object = object.pair, reduction = reduction)
. reduction <- paste0(reduction, ".l2")
. nn.reduction <- reduction
. }
. reduction.2 <- character()
. object.pair
. }, pca = {
. common.features <- intersect(x = rownames(x = Loadings(object = object.1[["pca"]])),
. y = rownames(x = Loadings(object = object.2[["pca"]])))
. object.pair <- merge(x = object.1, y = object.2, merge.data = TRUE)
. projected.embeddings.1 <- t(x = GetAssayData(object = object.1,
. slot = "scale.data")[common.features, ]) %*% Loadings(object = object.2[["pca"]])[common.features,
. ]
. object.pair[["projectedpca.1"]] <- CreateDimReducObject(embeddings = rbind(projected.embeddings.1,
. Embeddings(object = object.2[["pca"]])), assay = DefaultAssay(object = object.1),
. key = "projectedpca1_")
. projected.embeddings.2 <- t(x = GetAssayData(object = object.2,
. slot = "scale.data")[common.features, ]) %*% Loadings(object = object.1[["pca"]])[common.features,
. ]
. object.pair[["projectedpca.2"]] <- CreateDimReducObject(embeddings = rbind(projected.embeddings.2,
. Embeddings(object = object.1[["pca"]])), assay = DefaultAssay(object = object.2),
. key = "projectedpca2_")
. object.pair[["pca"]] <- CreateDimReducObject(embeddings = rbind(Embeddings(object = object.1[["pca"]]),
. Embeddings(object = object.2[["pca"]])), assay = DefaultAssay(object = object.1),
. key = "pca_")
. reduction <- "projectedpca.1"
. reduction.2 <- "projectedpca.2"
. if (l2.norm) {
. slot(object = object.pair[["projectedpca.1"]], name = "cell.embeddings") <- Sweep(x = Embeddings(object = object.pair[["projectedpca.1"]]),
. MARGIN = 2, STATS = apply(X = Embeddings(object = object.pair[["projectedpca.1"]]),
. MARGIN = 2, FUN = sd), FUN = "/")
. slot(object = object.pair[["projectedpca.2"]], name = "cell.embeddings") <- Sweep(x = Embeddings(object = object.pair[["projectedpca.2"]]),
. MARGIN = 2, STATS = apply(X = Embeddings(object = object.pair[["projectedpca.2"]]),
. MARGIN = 2, FUN = sd), FUN = "/")
. object.pair <- L2Dim(object = object.pair, reduction = "projectedpca.1")
. object.pair <- L2Dim(object = object.pair, reduction = "projectedpca.2")
. reduction <- paste0(reduction, ".l2")
. reduction.2 <- paste0(reduction.2, ".l2")
. }
. object.pair
. }, stop("Invalid reduction parameter. Please choose either cca or rpca"))
. internal.neighbors <- internal.neighbors[c(i, j)]
. anchors <- FindAnchors(object.pair = object.pair, assay = c("ToIntegrate",
. "ToIntegrate"), slot = slot, cells1 = colnames(x = object.1),
. cells2 = colnames(x = object.2), internal.neighbors = internal.neighbors,
. reduction = reduction, reduction.2 = reduction.2, nn.reduction = nn.reduction,
. dims = dims, k.anchor = k.anchor, k.filter = k.filter,
. k.score = k.score, max.features = max.features, nn.method = nn.method,
. eps = eps, verbose = verbose)
. anchors[, 1] <- anchors[, 1] + offsets[i]
. anchors[, 2] <- anchors[, 2] + offsets[j]
. return(anchors)
. })
3. future_xapply(FUN = FUN, nX = nX, chunk_args = X, args = list(...),
. get_chunk = `[`, expr = expr, envir = envir, future.globals = future.globals,
. future.packages = future.packages, future.scheduling = future.scheduling,
. future.chunk.size = future.chunk.size, future.stdout = future.stdout,
. future.conditions = future.conditions, future.seed = future.seed,
. future.lazy = future.lazy, future.label = future.label, fcn_name = fcn_name,
. args_name = args_name, debug = debug)
4. value(fs)
5. value.list(fs)
6. resolve(y, result = TRUE, stdout = stdout, signal = signal, force = TRUE)
7. resolve.list(y, result = TRUE, stdout = stdout, signal = signal,
. force = TRUE)
8. signalConditionsASAP(obj, resignal = FALSE, pos = ii)
9. signalConditions(obj, exclude = getOption("future.relay.immediate",
. "immediateCondition"), resignal = resignal, ...)
This worked with my previous installation with 2 conditions, and now it doesnt work even with 2 conditions. I am trying to use the merge protocol, but I think integration gives better results.
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4