timoast
commented
2 years ago Can you show the code used to process the two objects (Vd0 and PD0)? Can you also edit the issue title to be more informative?
Closed nonafarbehi2018 closed 2 years ago
timoast
commented
2 years ago Can you show the code used to process the two objects (Vd0 and PD0)? Can you also edit the issue title to be more informative?
nonafarbehi2018
commented
2 years ago Thank you Tim for your response. Here is the code that I used to generate PD0 (it is the same for VD0):
#######
#PD0
#######
library(dplyr)
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(patchwork)
library(ggplot2)
library(cowplot)
set.seed(1234)
setwd (c("/10X_3pGEX_ATAC/210701/ATAC/PD0/outs"))
outputdir <- "/nonfar/analysis/multiome_2"
#load RNA and ATAC data
counts <- Read10X_h5("filtered_feature_bc_matrix.h5")
names(counts)
fragpath <- "atac_fragments.tsv.gz"
# get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
genome(annotation) <- "hg38"
# create a Seurat object containing the RNA adata
pbmc <- CreateSeuratObject(
counts = counts$`Gene Expression`,
assay = "RNA"
)
# create ATAC assay and add it to the object
pbmc[["ATAC"]] <- CreateChromatinAssay(
counts = counts$Peaks,
sep = c(":", "-"),
fragments = "atac_fragments.tsv.gz",
annotation = annotation
)
################
#Quality control
################
DefaultAssay(pbmc) <- "ATAC"
pbmc <- NucleosomeSignal(pbmc)
pbmc <- TSSEnrichment(pbmc)
VlnPlot(
object = pbmc,
features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal"),
ncol = 4,
pt.size = 0
)
# filter out low quality cells
pbmc <- subset(
x = pbmc,
subset = nCount_ATAC < 100000 &
nCount_RNA < 25000 &
nCount_ATAC > 1000 &
nCount_RNA > 1000 &
nucleosome_signal < 2 &
TSS.enrichment > 1
)
#Peak calling
# call peaks using MACS2
peaks <- CallPeaks(pbmc, macs2.path = "/share/ScratchGeneral/nonfar/.conda/envs/r_4/bin/macs2")
# remove peaks on nonstandard chromosomes and in genomic blacklist regions
peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")
peaks <- subsetByOverlaps(x = peaks, ranges = blacklist_hg38_unified, invert = TRUE)
# quantify counts in each peak
macs2_counts <- FeatureMatrix(
fragments = Fragments(pbmc),
features = peaks,
cells = colnames(pbmc)
)
# create a new assay using the MACS2 peak set and add it to the Seurat object
pbmc[["peaks"]] <- CreateChromatinAssay(
counts = macs2_counts,
fragments = "atac_fragments.tsv.gz",
annotation = annotation
)
#Gene expression data processing
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc)
pbmc <- RunPCA(pbmc)
#DNA accessibility data processing
DefaultAssay(pbmc) <- "peaks"
pbmc <- FindTopFeatures(pbmc, min.cutoff = 5)
pbmc <- RunTFIDF(pbmc)
pbmc <- RunSVD(pbmc)
DepthCor(pbmc)
#Non-linear dimension reduction and clustering
pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
DimPlot(object = pbmc, label = TRUE) + NoLegend()
#Create a gene activity matrix
gene.activities <- GeneActivity(pbmc)
# add the gene activity matrix to the Seurat object as a new assay and normalize it
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
object = pbmc,
assay = 'RNA',
normalization.method = 'LogNormalize',
scale.factor = median(pbmc$nCount_RNA)
)
saveRDS(pbmc, paste0(outdir,"_PD0_line151.rds"))
You need to use the same peak set across the two different ATAC assays that you want to merge. I'd suggest either combining the peak sets as shown here or choosing one set of peaks and quantifying the peaks in both datasets.
Here you're also overwriting the RNA assay (containing the actual gene expression measurements from the multiome experiment) with the gene activity values from the ATAC assay.
It's unclear what assay is being used for anchor-finding, but since you removed the gene expression assay it's either the gene activity or ATAC assay. I would recommend keeping the actual gene expression measurements, and using that RNA assay for integration.
Hello, I am trying to aggregate 10X multiome data (GEX+ATAC) using Seurat and Signac and got this error
This is my current R session info
and this is the code I used to aggregate the 2 datasets:
Many thanks Nona