satijalab / seurat

R toolkit for single cell genomics
http://www.satijalab.org/seurat
Other
2.29k stars 913 forks source link

Issues with Integrating scATAC and scRNA seq Datasets #5809

Closed carversh closed 2 years ago

carversh commented 2 years ago

Hello! I am trying to integrate two sc datasets and am having issues with the integration. The datasets obtained are both pre-processed. Here is my code:

# making chromatin assay to ATAC seurat object
# get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
# reading in processed atac matrix
m <- readMM(file='/n/scratch3/users/s/shc989/matrix.mtx')
m <- m*1

# creating barcode and feature labels
barcodes <- fread(file='barcodes_atac.csv', header=F)[[1]]
features <- fread(file='peaks.csv', header=F)[[1]]

# adding these feature and barcode names to the matrix
rownames(m) <- features
colnames(m) <- barcodes

# Define GRanges object using the features
granges <- StringToGRanges(features, sep = c(":", "-"))
granges <- granges[as.vector(seqnames(granges) %in% standardChromosomes(granges)),]

# Create Seurat Chromatin Assay
chrom_assay <- CreateChromatinAssay(
  counts = m,
  sep = c(":", "-"),
  ranges = granges,
  genome = 'hg38',
  fragments = 'fragments.tsv.gz',
  min.cells = 0,
  min.features = 0,
  annotation = annotation
)

# read in metadata
metadata <- read.csv(
  file = "snATAC_metadta.csv",
  header = TRUE,
  row.names = 1
)

morab.atac <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

# process atac data

# We exclude the first dimension as this is typically correlated with sequencing depth
morab.atac <- RunTFIDF(morab.atac)
morab.atac <- FindTopFeatures(morab.atac, min.cutoff = "q0")
morab.atac <- RunSVD(morab.atac)
morab.atac <- RunUMAP(morab.atac, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key = "atacUMAP_")
# making rna seq object
##loading in 10x processed data, scaling and saving as RDS file
data_dir <- '/home/shc989/Gusev_Lab/Real_Data/data.dir'
list.files(data_dir) # Should show barcodes.tsv, genes.tsv, and matrix.mtx
expression_matrix <- Read10X(data.dir = data_dir, gene.column=1, strip.suffix = TRUE)

# read in metadata
meta = read.csv(file = '/n/scratch3/users/s/shc989/snRNA_metadta.csv')

# create seurat object
morab.rna = CreateSeuratObject(counts = expression_matrix)

morab.rna <- NormalizeData(morab.rna)
morab.rna <- FindVariableFeatures(morab.rna)
morab.rna <- ScaleData(morab.rna)
morab.rna <- RunPCA(morab.rna)
morab.rna <- RunUMAP(morab.rna, dims = 1:30)

# add metadata to seurat object
meta = read.csv(file = 'snRNA_metadta.csv', header = TRUE, row.names = 1)
morab.rna <- AddMetaData(object =  morab.rna, metadata = meta)

Now, when I go to integrate the data I get the following error and am not sure what to do: gene.activities <- GeneActivity(morab.atac, features = VariableFeatures(morab.rna))

`Error in names(cell.convert) <- cells: attempt to set an attribute on NULL Traceback:

  1. GeneActivity(morab.atac, features = VariableFeatures(morab.rna))
  2. FeatureMatrix(fragments = frags, features = transcripts, cells = cells, . verbose = verbose, ...)
  3. sapply(X = obj.use, FUN = function(x) { . SingleFeatureMatrix(fragment = fragments[[x]], features = features, . cells = cells, sep = sep, verbose = verbose, process_n = process_n) . })
  4. lapply(X = X, FUN = FUN, ...)
  5. FUN(X[[i]], ...)
  6. SingleFeatureMatrix(fragment = fragments[[x]], features = features, . cells = cells, sep = sep, verbose = verbose, process_n = process_n)`

I would really appreciate help on this! Thank you in advance.

timoast commented 2 years ago

Since this is an issue with running Signac:: GeneActivity(), please open an issue in the Signac repository and include the output of sessionInfo()