Replies: 2 comments 1 reply
-
|
My sessionInfo() |
Beta Was this translation helpful? Give feedback.
0 replies
-
|
Can you show the output of these lines: combined_merged
dim(combined_merged[['ATAC']][[]])
head(combined_merged[['ATAC']][[]])It seems like the feature metadata might be missing information for some peaks for some reason, or it has NA values (although this should show up when you run Can you also show the code used to create this object? By the way, here you are using all peaks with >0 counts for a group of cells in your motif enrichment test. More typically, we would do some differential testing to find significantly differentially accessible peaks with some greater stringency (minimum fold change and p-value thresholds). |
Beta Was this translation helpful? Give feedback.
1 reply
-
> combined_merged
An object of class Seurat
238720 features across 23684 samples within 5 assays
Active assay: ATAC (165621 features, 165621 variable features)
2 layers present: counts, data
4 other assays present: RNA, SCT, integrated, ACTIVITY
7 dimensional reductions calculated: pca, RNA_umap, tsne, lsi, RNA_integrated_umap, ATAC_integrated_umap, Multi_umap
> dim(combined_merged[['ATAC']][[]])
[1] 165621 20
> head(combined_merged[['ATAC']][[]])
count percentile AA AC AG AT CA CC
chr1-3094706-3095624 490 0.2142361 91 48 77 63 80 24
chr1-3113139-3114267 873 0.6265389 93 57 77 81 84 69
chr1-3119517-3120724 1373 0.8222327 99 44 81 99 72 56
chr1-3121251-3122029 493 0.2184989 77 48 77 48 80 40
chr1-3180850-3181790 550 0.3002759 85 43 69 81 53 31
chr1-3184713-3185642 468 0.1832678 39 31 83 46 61 81
CG CT GA GC GG GT TA TC TG TT
chr1-3094706-3095624 2 73 53 49 38 51 56 57 74 82
chr1-3113139-3114267 6 94 67 57 43 61 64 70 102 103
chr1-3119517-3120724 5 104 70 56 44 61 81 81 102 152
chr1-3121251-3122029 12 44 57 45 61 30 35 43 44 37
chr1-3180850-3181790 4 66 53 34 33 63 87 46 77 115
chr1-3184713-3185642 10 81 60 53 97 51 39 67 72 58
GC.percent sequence.length
chr1-3094706-3095624 40.26115 919
chr1-3113139-3114267 42.60407 1129
chr1-3119517-3120724 38.82450 1208
chr1-3121251-3122029 47.49679 779
chr1-3180850-3181790 35.81296 941
chr1-3184713-3185642 53.22581 930`To create the object, I had to combine 4 samples that has scRNAseq and acATACseq. I will only show 1 sample to have shorten version here. I tried Matchregionstat run on my combined_merged and got, For your comment about "all peaks with >0 counts", do you mean I should run FindMarkers first and feed those into features? (I did run FindMarkers and got DAs and saved as csv). # RNA part
# set sample ID
Samples <- c("PTZ_1hr", "PTZ_24hr", "SAL_1hr", "SAL_24hr")
# load h5 files and create RNA obejct
PTZ_1hr <- Read10X_h5(file.path(input_dir, Samples[1], "outs/filtered_feature_bc_matrix.h5"))
PTZ_1hr <- CreateSeuratObject(counts = PTZ_1hr$`Gene Expression`,
assay = "RNA",
project = "PTZ_1hr", # Sets orig.ident for all cells in this object
# min.cells = 3,
# min.features = 200,
meta.data = read.csv(file = file.path(input_dir, Samples[1], "outs/per_barcode_metrics.csv"), header = TRUE, row.names = 1))
PTZ_1hr <- PercentageFeatureSet(PTZ_1hr, pattern = "^mt-", col.name = "percent_mt_rna", assay = "RNA")
PTZ_1hr <-subset(PTZ_1hr, subset = nCount_RNA < 25000 & nCount_RNA > 1000 & percent_mt_rna < 5)
PTZ_1hr <-SCTransform(PTZ_1hr, method = "glmGamPoi", assay = "RNA", verbose = TRUE, vars.to.regress = c("percent_mt_rna"))
PTZ_1hr <- RunPCA(PTZ_1hr, assay = "SCT", verbose = TRUE)
PTZ_1hr <- RunUMAP(PTZ_1hr, dims = 1:30, reduction = "pca", assay = "SCT", verbose = TRUE)
PTZ_1hr <- RunTSNE(PTZ_1hr, dims = 1:30, reduction = "pca", assay = "SCT", verbose = TRUE)
identify <- function(seurat_obj, identity) {
seurat_obj$sample <- identity
seurat_obj$condition <- ifelse(grepl("PTZ", identity), "PTZ", "SAL")
seurat_obj$timepoint <- ifelse(grepl("1hr", identity), "1hr", "24hr")
return(seurat_obj)
}
PTZ_1hr <- identify(PTZ_1hr, Samples[1])
PTZ_1hr <- RenameCells(PTZ_1hr, add.cell.id = PTZ_1hr$sample)
# integrated RNA -
seurat_list <- list(PTZ_1hr = PTZ_1hr, PTZ_24hr = PTZ_24hr, SAL_1hr = SAL_1hr, SAL_24hr = SAL_24hr)
features <-SelectIntegrationFeatures(object.list = seurat_list, nfeatures = 3000)
seurat_list <-PrepSCTIntegration(object.list = seurat_list, anchor.features = features)
anchors <- FindIntegrationAnchors(object.list = seurat_list, normalization.method = "SCT", anchor.features = features)
integrated_rna <-IntegrateData(anchorset = anchors, normalization.method = "SCT")
integrated_rna <- RunPCA(integrated_rna, assay = "integrated", verbose = T)
integrated_rna <- FindNeighbors(integrated_rna, dims = 1:30, reduction = "pca", verbose = T, graph.name = c("RNA_nn", "RNA_snn"))
integrated_rna <- FindClusters(integrated_rna, graph.name = "RNA_snn", resolution = 0.8, verbose = T, cluster.name = "RNA_cluster_0.8")
integrated_rna <- PrepSCTFindMarkers(integrated_rna, verbose = T)
integrated_rna <- RunUMAP(integrated_rna, reduction = "pca", reduction.name = "RNA_umap", dims = 1:30, verbose = T)
integrated_rna <- RunTSNE(integrated_rna, reduction = "pca", dims = 1:30, verbose = T)
# for ATAC part
peak_beds <- c(
file.path(input_dir, Samples[1], "outs/atac_peaks.bed"),
file.path(input_dir, Samples[2], "outs/atac_peaks.bed"),
file.path(input_dir, Samples[3], "outs/atac_peaks.bed"),
file.path(input_dir, Samples[4], "outs/atac_peaks.bed"))
frag_paths <- c(
file.path(input_dir, Samples[1], "outs/atac_fragments.tsv.gz"),
file.path(input_dir, Samples[2], "outs/atac_fragments.tsv.gz"),
file.path(input_dir, Samples[3], "outs/atac_fragments.tsv.gz"),
file.path(input_dir, Samples[4], "outs/atac_fragments.tsv.gz"))merge_atac_id <-
function (peak_beds, min_peak_width = 20, max_peak_width = 10000, frag_paths, cell_ids, n_regions_simul = 2000, threads = 1){
message("Make sure that paths and cell_ids are in the same order")
message("peak_beds: ", paste(peak_beds, sep = " "))
message("frag_paths: ", paste(frag_paths, sep = " "))
message("Sample ID: ", paste(cell_ids, sep = " "))
# Reads peak BED files and reduces them using Signac::reduce to create a unified peak set.
reduced_peaks <- readr::read_tsv(peak_beds, comment = "#", col_names = c("chr", "start", "end"), show_col_types = FALSE) %>%
GenomicRanges::makeGRangesFromDataFrame() %>% Signac::reduce()
#Filters peaks by width (min_peak_width, max_peak_width) to remove noise or artifacts.
reduced_peaks <- reduced_peaks[reduced_peaks@ranges[1:3]@width >= min_peak_width & reduced_peaks@ranges[1:3]@width <= max_peak_width]
make_chrom_seurat <- function(i) {
# Creates a fragment object using Signac::CreateFragmentObject.
frags <- Signac::CreateFragmentObject(path = frag_paths[[i]])
# Computes a feature matrix using the reduced peaks and fragment files (Signac::FeatureMatrix).
counts <- Signac::FeatureMatrix(fragments = frags, features = reduced_peaks, process_n = n_regions_simul)
# Creates a chromatin assay and Seurat object.
obj <- Signac::CreateChromatinAssay(counts = counts, fragments = frags) %>% SeuratObject::CreateSeuratObject(assay = "ATAC", project = cell_ids[i])
return(obj)}
obj_list <- parallel::mclapply(seq_len(length(frag_paths)),
make_chrom_seurat, mc.cores = threads)
# Merges individual Seurat objects into a single Seurat object, labeling cells by their sample ID.
seurat_obj <- merge(x = obj_list[[1]], y = obj_list[-1], add.cell.ids = cell_ids)
return(seurat_obj)
}
merged_atac <- merge_atac_id(peak_beds = peak_beds, frag_paths = frag_paths, cell_ids = Samples, n_regions_simul = 100000, threads = 3)
mm10 <-import("/igm/apps/10X_chromium/refdata-cellranger-arc-mm10-2020-A-2.0.0/genes/genes.gtf.gz")
genome(mm10) <- "mm10"
seqlevelsStyle(mm10) <- "UCSC"
mm10$gene_biotype <- mm10$gene_type
annotations <- mm10
DefaultAssay(merged_atac) <- "ATAC"
Annotation(merged_atac) <- annotations
merged_atac <- NucleosomeSignal(merged_atac)
merged_atac <- TSSEnrichment(merged_atac, fast = FALSE)
merged_atac$high_tss <- ifelse(merged_atac$TSS.enrichment > 2, "High", "Low")
merged_atac$nucleosome_group <- ifelse(merged_atac$nucleosome_signal > 4,paste0("NS > ", 4), paste0("NS < ", 4))
merged_atac <- subset(merged_atac, subset = nCount_ATAC < 100000 & nCount_ATAC > 1000 & nucleosome_signal < 2 & TSS.enrichment > 2)
DefaultAssay(merged_atac) <- "ATAC"
merged_atac <- RunTFIDF(merged_atac) # Normalize data
merged_atac <- FindTopFeatures(merged_atac, min.cutoff = 10) # Select features
merged_atac <- RunSVD(merged_atac) # Generate 'lsi' reduction
merged_atac <- FindNeighbors(merged_atac, reduction = "lsi", dims = 2:30, graph.name = c("ATAC_nn", "ATAC_snn"), verbose = T)
merged_atac <- RunUMAP(merged_atac, reduction.name = "ATAC_umap", reduction = "lsi", dims = 1:30, verbose = T)
# Common cell find
common_cells <- intersect(colnames(integrated_rna), colnames(merged_atac))
integrated_rna <- integrated_rna[, common_cells]
merged_atac <- merged_atac[, common_cells]
integrated_rna[["ATAC"]] <- merged_atac[["ATAC"]]
atac_meta <- merged_atac@meta.data
rna_meta <- integrated_rna@meta.data
atac_specific_cols <-setdiff(colnames(atac_meta), colnames(rna_meta))
for (col in atac_specific_cols) {
integrated_rna[[col]] <- atac_meta[[col]]}
combined_merged <- integrated_rna
# combined process
DefaultAssay(combined_merged) <- "ATAC"
combined_merged <- RunTFIDF(combined_merged) # Normalize raw ATAC data
combined_merged <- FindTopFeatures(combined_merged, min.cutoff = 10) # Select variable peaks
combined_merged <- RunSVD(combined_merged) # Dimensionality reduction (creates "lsi")
DefaultAssay(combined_merged) <- "RNA"
combined_merged <- SCTransform(combined_merged, method = "glmGamPoi", assay = "RNA", vars.to.regress = c("percent_mt_rna"), vst.flavor = "v2", verbose = TRUE)
combined_merged <- RunPCA(combined_merged, verbose = TRUE) # Dimensionality reduction (creates "pca")
combined_merged <- FindNeighbors(combined_merged, dims = 1:30, assay = "SCT", reduction = "pca", graph.name = c("RNA_integrated_nn", "RNA_integrated_snn"), verbose = T)
combined_merged <- RunUMAP(combined_merged, reduction.name = "RNA_integrated_umap", assay = "SCT", reduction = "pca", dims = 1:30, verbose = T)
combined_merged <- FindNeighbors(combined_merged, dims = 2:30, assay = "ATAC", reduction = "lsi", graph.name = c("ATAC_integrated_nn", "ATAC_integrated_snn"), verbose = T)
combined_merged <- RunUMAP(combined_merged, reduction.name = "ATAC_integrated_umap", assay = "ATAC", reduction = "lsi", dims = 2:30, verbose = T)
# Multimodal
combined_merged <- FindMultiModalNeighbors(object = combined_merged,
reduction.list = list("pca", "lsi"), # reduction.list: Specifies RNA and ATAC reductions ("pca" for RNA, "lsi" for ATAC).
dims.list = list(1:30, 2:30), # dims.list: Defines the number of dimensions to use for each modality.
knn.graph.name = "Multimodal_wknn", # Multimodal_wknn: Weighted k-nearest neighbor graph.
snn.graph.name = "Multimodal_wsnn", # Multimodal_wsnn: Shared nearest neighbor (SNN) graph.
weighted.nn.name = "Multimodal_weighted.nn", # Multimodal_weighted.nn: Final multi-modal neighbor graph.
modality.weight.name = "Multimodal.weight", # modality.weight.name: Adds weights to indicate the relative contribution of RNA and ATAC.
verbose = TRUE)
combined_merged <- RunUMAP(object = combined_merged,
nn.name = "Multimodal_weighted.nn", # nn.name: Uses the joint neighbor graph.
verbose = TRUE,
reduction.name = "Multi_umap")
combined_merged <- FindClusters(combined_merged, graph.name = "Multimodal_wsnn", resolution = 0.8, verbose = T, cluster.name = "Multimodal_cluster_0.7")
# atac process
Annotation(combined_merged) <- keepStandardChromosomes(Annotation(combined_merged), pruning.mode = "coarse")
seq_lengths <- seqlengths(BSgenome.Mmusculus.UCSC.mm10)
seqlengths(Annotation(combined_merged)) <- seq_lengths[names(seqlengths(Annotation(combined_merged)))]
mcols(Annotation(combined_merged))$gene_id <- make.unique(as.character(mcols(Annotation(combined_merged))$gene_id))
combined_merged <- RegionStats(
object = combined_merged,
genome = BSgenome.Mmusculus.UCSC.mm10,
assay = "ATAC",
verbose = T
)
missing_gc <- rownames(combined_merged@assays$ATAC@meta.features)[
is.na(combined_merged@assays$ATAC@meta.features$GC.percent)
]
missing_seqnames <- sapply(strsplit(missing_gc, "-"), `[`, 1)
table(missing_seqnames)
all_seqnames <- sapply(strsplit(rownames(combined_merged@assays$ATAC@meta.features), "-"), `[`, 1)
valid_regions <- which(all_seqnames %in% seqlevels(Annotation(combined_merged)))
combined_merged@assays$ATAC <- subset(
combined_merged@assays$ATAC,
features = rownames(combined_merged@assays$ATAC@meta.features)[valid_regions]
)
# motif
pwm <- getMatrixSet(
x = JASPAR2020,
opts = list(collection = "CORE", tax_group = "vertebrates", species = 10090) # 10090 is the NCBI Taxonomy ID for mouse
)
combined_merged <- AddMotifs(
object = combined_merged,
genome = BSgenome.Mmusculus.UCSC.mm10, # Use the filtered genome
pfm = pwm
) |
Beta Was this translation helpful? Give feedback.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
It is a long post, sorry in advance! For context, I am doing data analysis from 10x Multiomics kit (scRNA and scATAC seq). I am new to this and this was hand down to me since the original person who did the exp left the lab. I am not savvy with coding nor r/python. I learnt all that while doing this analysis. (in R) I did scRNA integration/scATAC integration and all the fluff with scRNAseq part, at least i think so (up until GO analysis). Then I moved onto scATAC part. it was bit of struggle from the beginning tbh. the annotation was custom in our institution so i had to import gft files for annotation in the beginning and the peak alignment was a bit of process as well but I managed it to get all the process, integration and DAG so far. But when I tried to run Motif anlaysis i am having big issue that I can't fix for last 3 days... below is the code i am trying to run.
And this is basic info of my data.
The error I am getting is below.
And I did debugonce(FindMotifs) to see where is the issue comes from. Fisrt warning from here.
and because I have only 2 features in my background (i assume) this happens.
I think the problem is they don't have enough background features..? so i changed tried to use background.use to "all", default (gc content), and now using mannually putting high number (10000). but all not working..
I also, worried my GC.percent column is named not capitalized so tried to changed to capitalized one and that caused error on they can't find gc.percent so there's that.. I feel like i exhausted all options.. I sought out to get help within the institution but not many had done scATAC and also for those who did, they did not know my error (they never seen before). And now I am trying to feed DAG list for features as my last try.. But I would like to ask some other ideas on why it is like this..I am learning all tons of stuff but this is really frustrates me..(lol I enjoyed all the puzzles and this is fun ride for me but it is frustrating as well..)
Thank you for reading this really long post.. and please throw me any idea..like ANY! I will explore anything if I can fix this!
Beta Was this translation helpful? Give feedback.
All reactions