Using Speckle to assess differential abundance of cells from four different treatment groups within a cluster

[hannaminns]

Hello!

I'm hoping to use your tool to assess differential abundance of cells from different treatment groups within a cluster.

For our project, we have four different treatment groups (tx1, tx2, tx3, tx4), each tagged with a hashtag antibody so that they can be run together on the 10X.

We ran four batches, each with (1) replicate of each treatment group like so:
Batch 1 - tx1, tx2, tx3, tx4
Batch 2 - tx1, tx2, tx3, tx4
Batch 3 - tx1, tx2, tx3, tx4
Batch 4 - tx1, tx2, tx3, tx4
Total = 16 samples, 4 replicates for each tx group per 4 batches

I have also identified my cell types through Seurat clustering, so I have identities for "txgroup", "orig.ident" (batch ID), and "celltype" in my Seurat object.

i.e.
clusters = seurat_obj$celltype
sample = seurat_obj$orig.ident
group = seurat_obj$txgroup

I want to assess statistically if, for example, "in cluster 1, are there more cells from tx1 compared to the other three tx groups"? And look at that sort of differential abundance testing based on treatment group for all clusters in my dataset.

Is it feasible to use Speckle for this given my experimental design and question? If so, do you have any guidance on how to run the function? I tried to run propeller.anova(), but must not be calling/understanding the parameters correctly because it is not working.

Thank you!!

Best,
Hanna

[bphipson]

Hi Hanna You have a great experimental design! Would you like to take Batch into account in the analysis as well? I'm wondering exactly what comparisons you would like to do in this experimental design set-up. propeller.anova will test for ANY differences between the four treatment groups. If you want to specifically test tx1 vs (tx2,tx3,tx4) then you will need to set up a contrast to do that and then call propeller.ttest. I have an example of how to do that in the speckle vignette. First, you need to get the transformed proportions: props <- getTransformedProps(seurat_obj$celltype, seurat_obj$orig.ident, transform="logit") This will output a list object containing the cell type counts, proportions and transformed proportions. Based on the columns of these matrices, you will need to identify the batch and treatment information for each of your samples (i.e. the columns of the matrix). The rows are the cell types. Assuming the batch and treatment variables have been created that match the order of the columns of the proportions matrix, you need to then set up the design matrix as follows: design <- model.matrix(~0 + treatment + batch) The first four columns will correspond to the treatment. The remaining columns will correspond to batch. Assuming the columns of the design matrix are labeled Tx1, Tx2, Tx3 Tx4, you then need to set up your contrast matrix (i.e. the comparison you want to do) mycontr <- makeContrasts(Tx1 - (Tx2 + Tx3 + Tx4)/3, levels=design) This contrast will compare Tx1 to the average of the other three treatments. You can set up Tx2, Tx3 and Tx4 vs the rest in the same way. In the current implementation you can only test one contrast at a time. Next you test the contrast using propeller.ttest: propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE) Hope this helps. Belinda

[hannaminns]

Hi Belinda,

Thanks so much for your response. I've reviewed it and the vignette in more detail and got an output, but I'm not sure if I did it correctly.

This data has been batch corrected using Harmony, but I assume I should still take in account batch when looking at differential abundance. Does that make sense?

For props <- getTransformedProps() I ended up with 16 sample columns because for each batch, there are cells pertaining to each of the four treatment groups:
i.e. the first four columns of props are:
batch1_tx1, batch1_tx2, batch1_tx3, batch1_tx4
And then on for batches 2, 3, and 4

treatment <- rep(c("Tx1","Tx2","Tx3", "Tx4"),each=4)
batch <- rep(c(1,2,3,4),4)
design <- model.matrix(~ 0 + treatment + batch)

Then first ran:
propeller.anova(prop.list=props, design=design, coef = c(1,2,3,4), robust=TRUE, trend=FALSE, sort=TRUE)

And also did the ttest:
mycontr <- makeContrasts(treatmentTx1 - (treatmentTx2+ treatmentTx3 + treatmentTx4)/3, levels=design)
propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE)

However, for both the anova and also the ttest (tested a few different contrasts), I'm getting no p values < 0.05, which should not be the case for all the clusters. Looking at the raw counts coming from each batch and each treatment group, there should definitely be some differential abundance within certain clusters that is statistically significant.

I'm wondering if I set this up wrong. Is there anything I need to change?

Best,
Hanna

[apal6]

Hi,

I have three treatment groups (Untreated, treated and mixed) and untreated and mixed has 4 samples and treated has 3 sample. I want to test for the sample specific proportion of celltypes in each cluster. Not sure what am I doing incorrectly but here is my script and the error:

props <- getTransformedProps(clusters = Idents(seurat),
sample = seurat$orig.ident)

#clusters <- Idents(seurat) # Clusters information
#sample <- seurat$orig.ident # Sample type information

props <- getTransformedProps(clusters = clusters, sample = sample)

treatment <- c(rep("Untreated", 4), rep("treatment", 3), rep("mixed_treatment", 4))
batch <- rep(c(1, 2,3), times = c(4, 3, 4))
design <- model.matrix(~ 0 + treatment + batch)
design

propeller.anova(prop.list=props, design=design, robust=TRUE, trend=FALSE, sort=TRUE)

Coefficients not estimable: batch
Error: (converted from warning) Partial NA coefficients for 10 probe(s)

What should I change?

Thanks in advance,
Aastha