timoast
commented
10 months ago Hi, can you show the full code you're running, and the output of head(VariableFeatures(data.rna))?
Could you also try updating to the latest Signac release?
Closed YaoLi3 closed 10 months ago
timoast
commented
10 months ago Hi, can you show the full code you're running, and the output of head(VariableFeatures(data.rna))?
Could you also try updating to the latest Signac release?
YaoLi3
commented
10 months ago the output of head(VariableFeatures(data.rna)) is:
> head(VariableFeatures(data.rna))
[1] "nbisL1-mrna-12678" "MSTRG.16655" "nbisL1-mrna-2064"
[4] "nbisL1-mrna-16480" "MSTRG.15682" "nbisL1-mrna-4839" the full code I'm running:
library(Seurat)
library(Signac)
library(Matrix)
# load scRNA data
fn <- 'liver.rds'
data.rna <- readRDS(fn)
# load scATAC data
load_atac <- function(directory_path){
barcode_fn <- file.path(directory_path, "barcodes.tsv")
peaks_fn <- file.path(directory_path, "peak.bed")
matrix_fn<- file.path(directory_path, "matrix.mtx")
# Load the data files
barcodes <- read.table(barcode_fn, sep = "\t", header = FALSE)
peaks <- read.table(peaks_fn, sep = "\t", check.names=FALSE)
frag_matrix <- readMM(matrix_fn)
# Set the row names for the matrix from the peaks object
rownames(frag_matrix) <- peaks$V1
# Set the column names for the matrix
colnames(frag_matrix) <- barcodes$V1
return (frag_matrix)
}
dir2 <- 'liver-2_1'
peaks <- load_atac(dir2)
# 2. chromain assay
metadata <- read.table('Metadata.tsv',sep='\t',header =TRUE)
frag_fn <- 'liver-2_1.fragments.tsv.gz'
chrom_assay <- CreateChromatinAssay(
counts = peaks,
sep = c(":", "-","_"),
fragments = frag_fn,
min.cells = 10,
min.features = 200
)
data.atac <- CreateSeuratObject(
counts = chrom_assay,
assay = "peaks",
meta.data = metadata
)
data.atac$tech <- 'atac'
# preprocessing
data.atac <- FindVariableFeatures(data.atac)#, assay = "ACTIVITY", selection.method = "vst", nfeatures = 2000)
data.atac <- NormalizeData(data.atac)
data.atac <- ScaleData(data.atac)
# ATAC analysis add gene annotation information
gff_fn = 'final.gff3'
gff = rtracklayer::import.gff(gff_fn)
gene_gff = gff[gff$type=='gene']
Annotation(data.atac) <- gene_gff
# Dimension Reduction of the peak matrix via Latent Semantic Indexing
# To study the internal structure of the scATAC-seq data, making anchor detection possible
VariableFeatures(data.atac) <- names(which(Matrix::rowSums(data.atac) > 100))
# We exclude the first dimension as this is typically correlated with sequencing depth
data.atac <- RunTFIDF(data.atac)
data.atac <- FindTopFeatures(data.atac, min.cutoff = "q0")
data.atac <- RunSVD(data.atac)
data.atac <- RunUMAP(data.atac, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key = "atacUMAP_")
# Identifying anchors between scRNA-seq and scATAC-seq datasets
# transfer label
# quantify gene activity
gene.activities <- GeneActivity(data.atac, features = VariableFeatures(data.rna))As for updating the Signac package, I'm afraid this is going to be difficult. I had a hard time installing the package on the Linux server my institute provided, so I'm currently using a R miniconda environment my colleague kindly gave to me. It's not very likely to get her to update it for me. Do I need to do this?
timoast
commented
10 months ago I think the issue here is you don't have biotypes information in the gene annotations that you're using. If you add a column like this:
gene_gff$biotype <- "protein_coding"can you see if that fixes the problem?
YaoLi3
commented
10 months ago > gene_gff$biotype <- "protein_coding"
> Annotation(data.atac) = gene_gff
> gene.activities <- GeneActivity(data.atac, features = VariableFeatures(data.rna))
Error in `[[<-`(`*tmp*`, name, value = logical(0)) :
0 elements in value to replace 22094 elementstried this, still getting the same error
timoast
commented
10 months ago sorry it should be gene_biotype not biotype. Try:
gene_gff$gene_biotype <- "protein_coding"sorry but still... Really don't understand where this error is from
Error in `[[<-`(`*tmp*`, name, value = logical(0)) :
0 elements in value to replace 22094 elementsYou also need to add gene_name information to your annotations since you're setting the features parameter in GeneActivity().
gene_gff$gene_name <- gene_gff$gene_idI think this should resolve the issue, but if not please take a look at the example gene annotations included in the atac_small dataset in the Signac package (Annotation(atac_small)), and ensure that the annotations you're using matches the format of those annotations. If they do match and you think this is a bug in Signac, please reopen this issue.
YaoLi3
commented
10 months ago yes, you are right, adding gene_name column into my annotations solved the issue. I falsely assumed GeneActivity() only uses gene_id information.
Thanks again your help!
Dear Signac team,
Hi, I'm trying to perform label transferring between ATAC-seq data and annotated scRNA-seq data. However, I keep getting this error when trying to call GeneActivity on this.
I don't understand the meaning of this error message. Does this mean the feature names (gene names of RNA data) can not match the gene names in atac data? Gene names only exist in the annotation of atac data object by the way.
I have tried to ensure the cell and gene naming convention is consistent between the RNA data and ATAC data, but I'm still getting the same error.
Here is my gff annotation:
And this is my session info:
Thanks in advance!
Best, Yao