Hi,

I've been playing around with different normalization strategies for scRNA-seq data. Contrary to the header, I think, I have actually two questions, the first one being: Is the computeSumFactors() philosophy really applicable to current Drop-seq data sets?

The sample that initiated the question was generated with Drop-seq, i.e., it covers around 3,000 cells, but with fairly low coverage:

> pData(sceset)$total_features %>% summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
     17     415    1144    1304    1988    6617

The average expression values for each gene across all those cells are quite low, of course, so in order to avoid getting size factors of zero, I use a very small subset of genes (around 370) where the average count across all cells is greater than 1. (I do not filter out cells with low gene counts, which may also be worthy of a discussion). I then use these size factors to normalize the entire data set.

While I was somewhat satisfied with the tSNE results of the normalized counts, I eventually noticed that the data set contained a couple of cell pairs that had exactly the same counts for all their genes. This is clearly an unwanted artifact and I am going to exclude these duplicate cells in the future, but I noticed that computeSumFactors() had assigned quite different size factors to these cells.

For example:

> sizeFactors( sceset[, colnames(duplicated_cells))[c(1:2)] ] )

 I1_1     I2_1
2.646779 5.142105

where I1 and I2 are different cells with the exact same counts for about 5,000 genes. The size factors were quite different, regardless of whether I used the cluster parameter of the function or not.

I just have trouble grasping why that would be - again, I will definitely exclude duplicate cells from future analyses, but I was just trying to get a better understanding of why the sum factors would be so different.

I appreciate any insights - also about whether the approach is actually applicable to Drop-seq data, I have the feeling that a lot of the filtering strategies that are shown in tutorials (such as the one for scater and scran) are not really realistic for Drop-seq data with its fairly low coverage of genes, but relatively high numbers of cells.

Thanks a lot!

Friederike

Even if you have two cells with exactly identical counts, they should get slightly different size factors because they get placed into different pools. A (very small) proportion of equations in the linear system will differ between the two cell copies, which results in slight differences to the estimates. However, these differences should indeed be small; I can't think of a reason for why you should see a >2-fold change between copies. I can only suggest:

• Remove cells with low read counts. These cells don't contribute much information, while their size factors still require estimation. As a result, the solution of the linear system is less stable, such that the few equations that differ between the cell copies may have a large impact on the final estimates.
• The filtering is more aggressive than I would have done, but I doubt this makes much of a difference.
• Make sure you haven't set sizes lower than 20; below that, the pooled counts are not stable.

I've used computeSumFactors for a number of 10X data sets, and it seems to work okay. I usually plot the size factors against the library sizes as a diagnostic, most cells should be scattered around the diagonal (depending on the heterogeneity of the data set).

P.S. If none of that works, you can send me a MWE offline - this behaviour is interesting and should be resolved.

Edit: I'll summarize the offline discussion I had with Friederike after looking at the data. Basically, the problem was that, for the majority of cells, the count matrix was simply too sparse. Most cells express fewer than 700 genes, out of around 6000 genes that are expressed at an average count greater than 0.1. More than 90% of the counts are zero, such that the pooled size factors are also likely to be zero.

The solution was to QC more aggressively, e.g., remove cells with fewer than 500 expressed features. This removes the cells that were previously getting negative size factors, which allows us to avoid using positive=TRUE (see my comments below). In addition, the precision of the size factor estimates can be improved by increasing the number of window sizes, e.g., sizes=2:10*10. After applying these two changes, the cell copies have very close (+/- 5%) estimates, which is what we'd expect.

Of course, this means that we lose those cells with very few expressed features. However, I would argue that all normalization methods will struggle with these cells. TMM and DESeq will return gibberish because there are far too many zeroes. Library size normalization will be entirely driven by a handful of highly expressed genes - composition biases aside, the estimates will be highly imprecise.