Integrating two Seurat Objects

Question

[roberts]

Hello- I am trying two integrate two seurat objects of pbmc but am having some issues

 Sample1 <- readRDS("Desktop/final_data/Sample1.RDS")
Sample1$sample <- "Sample 1"

Sample2 <- readRDS("Desktop/final_data/Sample2.RDS")
Sample2$sample <- "Sample 2"

patient.merge <- merge(Sample1, y=c(Sample2))
DefaultAssay(patient.merge)<- "RNA"
patient.list <- SplitObject(patient.merge, split.by="sample")
patient.list <- patient.list[c( "Sample 1", "Sample 2")]
for (i in 1:length(patient.list)) {
  patient.list[[i]] <- SCTransform(patient.list[[i]], verbose = FALSE)

}
patient.features <- SelectIntegrationFeatures(object.list = patient.list, nfeatures = 3000)
patient.list <- PrepSCTIntegration(object.list = patient.list, anchor.features = patient.features, 
                                   verbose = FALSE)
patient.anchors <- FindIntegrationAnchors(object.list = patient.list, normalization.method = "SCT", 
                                          anchor.features = patient.features, verbose = FALSE)

patient.integrated <- IntegrateData(anchorset = patient.anchors, normalization.method = "SCT", 
                                    verbose = FALSE)

patient.integrated <- ScaleData(patient.integrated, verbose = FALSE)
unwanted_genes2 <- grep(pattern = "^HLA*|^IGHV*|^IGHJ*|^IGHD*|^IGKV*|^IGLV*|^TRBV*|^TRBD*|^TRBJ*|^TRDV*|^TRDD*|^TRDJ*|^TRAV*|^TRAJ*|^TRGV*|^TRGJ*", x = patient.integrated@assays$[e-mail], value = T)
remove_genes2 <- patient.integrated @assays$[e-mail] %in% unwanted_genes2
patient.integrated @assays$[e-mail] = patient.integrated@assays$[e-mail][!remove_genes2]

patient.integrated <- RunPCA(patient.integrated, features = VariableFeatures(patient.integrated), npcs = 30, verbose=F)
patient.integrated <- JackStraw(patient.integrated, num.replicate = 100)
patient.integrated <- ScoreJackStraw(patient.integrated, dims = 1:20)
ElbowPlot(patient.integrated)

patient.integrated <- FindNeighbors(patient.integrated, dims = 1:15)
patient.integrated <- FindClusters(patient.integrated, resolution = 0.5)
patient.integrated <-RunUMAP(patient.integrated, dims=1:15)
DimPlot(patient.integrated, reduction= "umap")
DimPlot(patient.integrated, reduction= "umap", group.by = "sample")

patient.integrated.markers <- FindAllMarkers(patient.integrated, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
patient.integrated.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
FeaturePlot(patient.integrated, features = c("CD3G","CD8A", "CD4","CTLA4" ,"FOXP3", "CCR7","CD27", "KLRG1",
                                             "CCR5", "PRF1", "GZMA", "GZMB", "IFNG", "CD19", "PDCD1", "CXCR5",
                                             "KLRC1", "KLRC3","FCGR2A","CSF1R", "FLT3", "CLEC4C", "HBA", "HBB",
                                             "HBD", "IL6", "IL17", "EOMES"))
top10 <- patient.integrated.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(patient.integrated, features = top10$gene, group.by = "sample") + NoLegend()

Specifically the issues I have are that when I run integrate dataI get the warning -- adding a command log without an assay associated with it and when I run feature plot I get

    Warning: Found the following features in more than one assay, excluding the default. We will not include these in the final dataframe: CD4, FLT3, IL6
Warning message:
In FetchData(object = object, vars = c(dims, "ident", features),  :
  The following requested variables were not found: CD4, FLT3, IL6, IL17
    enter code here

which does not make sense because CD4 and FLT3 are found in the original seurat objects?

Any help would be greatly appreciated! Thanks!

Answer 1

[igor]

Each object can have multiple assays. Make sure you specify the correct assay. Most commands have as assay parameter. It should probably be RNA instead of integrated. You can check what the default one is with DefaultAssay().

[roberts]

Thanks, I changed my code to

 DefaultAssay(patient.integrated)<- "RNA"
patient.integrated.markers <- FindAllMarkers(patient.integrated, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
patient.integrated.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
FeaturePlot(patient.integrated, features = c("CD3G","CD8A", "CD4","CTLA4" ,"FOXP3", "CCR7","CD27", "KLRG1",
                                             "CCR5", "PRF1", "GZMA", "GZMB", "IFNG", "CD19", "PDCD1", "CXCR5",
                                             "KLRC1", "KLRC3","FCGR2A","CSF1R", "FLT3", "CLEC4C", "HBA", "HBB",
                                             "HBD", "IL6", "IL17", "EOMES"))
top10 <- patient.integrated.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(patient.integrated, features = top10$gene, group.by = "sample") + NoLegend()

but now I get an error when I try to produce a heatmap saying--

Error in DoHeatmap(patient.integrated, features = top10$gene, group.by = "sample") : 
  No requested features found in the scale.data slot for the RNA assay

.

[igor]

The scale.data slot only has the variable genes by default. It seems none of your genes were part of that list.