I need help please

[beginner984]

Hello

I need help please if possible. I have tried a lot of ways

I have 3' multiome 10X from which I have analyzed scRNA-seq part (cancer and control)

> cancer
An object of class Seurat 
36601 features across 18338 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)
> control
An object of class Seurat 
36601 features across 8135 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)
> 

> pbmc<- merge(control, y = cancer)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 20)

> pbmc
An object of class Seurat 
36601 features across 12812 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)
>

As I have cancer and control (although from the same lab and same kit), I integrated the data from SCtransform part like

ifnb.list <- SplitObject(pbmc, split.by = "Condition")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
                                         anchor.features = features)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")

immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)

> immune.combined.sct <- FindNeighbors(immune.combined.sct, dims = 1:30)
Computing nearest neighbor graph
Computing SNN
> immune.combined.sct <- FindClusters(immune.combined.sct, resolution = 0.7)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)

> immune.combined.sct
An object of class Seurat 
51025 features across 12812 samples within 3 assays 
Active assay: integrated (3000 features, 3000 variable features)
 2 other assays present: RNA, SCT
 2 dimensional reductions calculated: pca, umap
>

My question please is: why my heatmap is this unclear with bad blocks? I don't know where I am doing wrong

> DefaultAssay(immune.combined.sct)
[1] "integrated"
> 

Please could inspect my approach to see how I can get a reasonable data because seemingly somewhere I am doing wrong but I can not detect where

[leiliang98]

Hi,

Just looking at your QC plot, the percent.mt is a little high to me. It should be about 5% as shown in the Guided Clustering Tutorial of PBMCs: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

This makes me wonder if there is something wrong with the dataset you are using, specifically too many dead cells.

Good luck!

[cswoboda]

Suggestions:

1. If you leave in the same cutoff for percent.mt, you should regress it out in the PCA by using vars.to.regress = "percent.mt" in the SCTransform command. Might be throwing off the PCA and giving you low quality clusters.
2. This is multiome data, but have you processed it as such? It looks like you have quite a lot of features in your dataset in the RNA Assay. Did you preprocess them like shown here, making sure your antibody assay and gene expression assay are seperate?
   https://satijalab.org/seurat/articles/multimodal_vignette.html
3. Have you tried switching back to the SCT or RNA assay for your heatmap? I normally don't do any visualizations with the integrated assay, as it's a subset of genes (3000 total). You should instead generate markers conserved markers across the two conditions using FindConservedMarkers() if you want a cleaner heatmap with genes shared from both datasets.

[beginner984]

Suggestions:

1. If you leave in the same cutoff for percent.mt, you should regress it out in the PCA by using vars.to.regress = "percent.mt" in the SCTransform command. Might be throwing off the PCA and giving you low quality clusters.
2. This is multiome data, but have you processed it as such? It looks like you have quite a lot of features in your dataset in the RNA Assay. Did you preprocess them like shown here, making sure your antibody assay and gene expression assay are seperate?
   https://satijalab.org/seurat/articles/multimodal_vignette.html
3. Have you tried switching back to the SCT or RNA assay for your heatmap? I normally don't do any visualizations with the integrated assay, as it's a subset of genes (3000 total). You should instead generate markers conserved markers across the two conditions using FindConservedMarkers() if you want a cleaner heatmap with genes shared from both datasets.

Thank you so much

I don't know how to add vars.to.regress = "percent.mt" to this command ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)

[cswoboda]

I don't know how to add vars.to.regress = "percent.mt" to this command ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)

This should work by applying:
ifnb.list <- lapply(X=infb.list, FUN = SCTransform, vars.to.regress = "percent.mt")
lapply should pass the vars.to.regress argument to the SCTransform function, as it passes any argument's it doesn't recognize to the function you're applying to your list object.

optionally, you can write a for loop:

for (i in length(ifnb.list)) {
ifnb.list[[i]] <- SCTransform(ifnb.list[[i]], vars.to.regress = "percent.mt")
}

either or should work in your case

[beginner984]

Thank you so much

About this comment of you making sure your antibody assay and gene expression assay are separate

I am not following the point of this comment :(

I have multiome but I am only analysing the scRNA part, so still do I need to go to this part?

Thanks again

[cswoboda]

Yes, this is likely where your data is getting confused. Can I ask what type of multiome data you have?

In the seurat vignette, it's assuming ADT and RNA data.:

`# Load in the RNA UMI matrix

Note that this dataset also contains ~5% of mouse cells, which we can use as negative

controls for the protein measurements. For this reason, the gene expression matrix has

HUMAN_ or MOUSE_ appended to the beginning of each gene.

cbmc.rna <- as.sparse(read.csv(file = "../data/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz", sep = ",",
header = TRUE, row.names = 1))

To make life a bit easier going forward, we're going to discard all but the top 100 most

highly expressed mouse genes, and remove the 'HUMAN_' from the CITE-seq prefix

cbmc.rna <- CollapseSpeciesExpressionMatrix(cbmc.rna)

Load in the ADT UMI matrix

cbmc.adt <- as.sparse(read.csv(file = "../data/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz", sep = ",",
header = TRUE, row.names = 1))

Note that since measurements were made in the same cells, the two matrices have identical

column names

all.equal(colnames(cbmc.rna), colnames(cbmc.adt))`

As you can see, they load an expression matrix for the rna data "cbmc.rna" and one for the adt data "cbmc.adt"

They then create the seurat object by just using the RNA counts:

creates a Seurat object based on the scRNA-seq data

cbmc <- CreateSeuratObject(counts = cbmc.rna)

[beginner984]

Thank you so much my data is multiome of scATAC-seq + scRNA-seq for which I read filtered_feature_bc_matrix.h5 and use The RNA part of that like

> inputdata.10x <- Read10X_h5("filtered_feature_bc_matrix.h5")
Genome matrix has multiple modalities, returning a list of matrices for this genome
Warning message:
In sparseMatrix(i = indices[] + 1, p = indptr[], x = as.numeric(x = counts[]),  :
  'giveCsparse' has been deprecated; setting 'repr = "T"' for you
> # extract RNA and ATAC data
> rna_counts <- inputdata.10x$`Gene Expression`
> atac_counts <- inputdata.10x$Peaks
> 
> # Create Seurat object
> pbmc <- CreateSeuratObject(counts = rna_counts)
> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
> 
> 

So still I need this process ?

[cswoboda]

So still I need this process ?

Looks like you did this fine. My recommendation is this:

1. Run your pipeline again, normalizing and regressing out percent.mt.
2. Subset your nFeature_RNA a little higher (greater than or equal to 500 for exmaple) maybe, which should bring cells into clusters with higher signals.
3. Post Integration, run FindMarkers using the SCT assay, not the integrated assay.
4. Create your heatmap using the top 20-30 markers for each cluster. You can sort by min.pct1 if you want genes abundandantly expressed in each celltype.

I'd say 3-4 are vital here, as if you used the integrated assay for your heatmap you're likely only testing the 3000 genes you used for clustering.

Based on your UMAP, it also looks like some of your clusters could be merged. 0, 1, 2, 3 for example seem to have markedly similar expression profiles based on your heatmap.

[beginner984]

Thanks a lot

Sorry, I have 4 cancers and 4 controls (from the same lab and protocol though)

This is how I perform

I merge cancer and control Seurat objects before doing any further step, next I use the merged data for integration process like

> cancer
An object of class Seurat 
36601 features across 18338 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)
> control
An object of class Seurat 
36601 features across 8135 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)
> 

> pbmc<- merge(control, y = cancer)

You think this is fine or I should firstly run SCTransform on each Seurat object before integration?

Thank you so much again

[cswoboda]

This is going to be different depending on your data. When you merge control and cancer, and run the standard seurat clustering, do they seperate by batch/replicate, or do celltypes remain? I tend to treat each replicate as it's own dataset and run SCTransform on each individually, and integrate all of them in this manner. In truth, I think this is the most responsible way to normalize data prior to integration, as it won't lead to the blinding or heightening of any technical differences intra-replicate prior to combined analysis and clustering. Different datasets have different levels of noise and technical affects, I work with nuclear data which I believe can have a lot higher technical artefacts between replicates.

[beginner984]

Thank you so much

This is the plot before integration

And this is after integration in which they are overlaid well except for some parts

[cswoboda]

You're interested in the visualization of just the control and condition objects prior to integration (how does just control look merged when visualized by replicate, how does experimental condition look merged when visualized by replicate)

[beginner984]

Thank you so much. Looking at my first plot, as cancer and control groups seem separated, so should I pass them in SCTransform individually before integration, right?

Actually I did not follow your last comment well :(

[cswoboda]

Sorry, I know this is a little confusing.
Let's say you have three replicates for control and three for cancer.
If you make the objects:
ctrl_1 <- CreateSeuratObject(ctrl_counts_1)
ctrl_2 <- CreateSeuratObject(ctrl_counts_2)
ctrl_3 <- CreateSeuratObject(ctrl_counts_3)

cancer_1 <- CreateSeuratObject(cancer_counts_1)
cancer_2 <- CreateSeuratObject(cancer_counts_2)
cancer_3 <- CreateSeuratObject(cancer_counts_3)

and run the normalization, scaling.

Then you merge both conditions seperately
ctrl_merge <- merge(ctrl_1, c(ctrl_2,_ctrl_3)
cancer_merge <- merge(cancer_1, c(cancer_2,_cancer_3)

RunPCA, FindClusters, FindNeighbors, RunUMAP for these merged objects.

If the resulting umap for the merged objects are clustered by celltype, and not by replicate, you're likely good to integrate just ctrl_merge and cancer_merge. If not, you should run SCTransform on each individual object (ctrl_1, ctrl_2, ctrl_3, cancer_1, cancer_2, cancer_3), and integrate objects seperately.

Frankly, I tend to always just integrated my objects seperately by treating each replicate as it's own dataset, and believe this is more common. But if you want to merge prior, this is the correct way to do that prior to integration.

Does this clear things up?