Wrong avg_logFC?

[JerryZhang-1222]

Dear developer,

I've been using Libra for several projects and I noticed that the avg_logFC does not equal to log(group1.exp/grou2.exp).

It's also weird that in this screenshot, the first row, foldchange equals to 5.4788/1 while the second row equals to 1.9859155/0.3636364 = 5.461267, which is smaller. However, the logFoldchange of the second row is bigger than that in the first row.

I wonder how's that happened and how did foldchange been calculated? Thank you very much!

Best,
Feiyang.

[mph270]

Hello,

I am also getting odd results with avg_logFC. For some rows, the directionality of avg_logFC matches that of the change in expression (WT.exp vs KO.exp) but for other rows it is the opposite directionality (for example below, rows 1 and 3). N.B. This is a separate issue from the assignment of the reference group.

Thank you for your help!

Best,
Molly

[mph270]

Here is another example after making genotype and animal ID (orig.ident) factor levels both alphabetical and filtering out lowly expressed genes from the results. The avg_logFC value for the third row seems very off.

I used the following code:

DE_edgeR_LRT = run_de(mapped_seurat, cell_type_col = "new_ident", label_col = "genotype", replicate_col = "orig.ident")

Thanks again for your help and for creating this package!

-Molly

[AlanTeoYueYang]

Hi, with the latest pull, you set factor for reference in your log fold change calculations else it would use label1 and label2 alphabetically. The log fold-change calculations are made with Seurat here.

[IanevskiAleksandr]

I guess for all methods (single-cell AND pseudobulk), Libra calculates the avg_logFC values using Seurat::FoldChange function:

rtmp_stats = Seurat::FoldChange(sc, label1_barcodes, label2_barcodes, base=exp(1))

Also explicit factor level setting will ensure consistent fold change directions in your analysis, regardless of the internal calculation method.

[nigiord]

I guess for all methods (single-cell AND pseudobulk), Libra calculates the avg_logFC values using Seurat::FoldChange function:

rtmp_stats = Seurat::FoldChange(sc, label1_barcodes, label2_barcodes, base=exp(1))

Also explicit factor level setting will ensure consistent fold change directions in your analysis, regardless of the internal calculation method.

Not really, the result from Seurat::FoldChange() is actually droped at line 292 here (for pseudobulk methods) :

Libra/R/run_de.R

Lines 284 to 293 in 388d41c

	} else {
	DE %<>%
	# calculate adjusted p values
	group_by(cell_type) %>%
	mutate(p_val_adj = p.adjust(p_val, method = 'BH')) %>%
	# make sure gene is a character not a factor
	mutate(gene = as.character(gene)) %>%
	dplyr::select(cell_type, gene, p_val, p_val_adj, de_family, de_method, de_type, avg_logFC) %>%
	dplyr::left_join(out_stats %>% dplyr::select(-avg_logFC), by = c('gene', 'cell_type'))
	}

The reported avg_logFC is actually from pseudobulk_de :

Libra/R/run_de.R

Lines 269 to 271 in 388d41c

	'avg_logFC' = 'log2FoldChange', ## DESEeq2
	'avg_logFC' = 'logFC', ## limma/edgeR
	'avg_logFC' = 'avg_log2FC' # Seurat V4

They do not match at all in my case... (wrong value, wrong directionnality) I’m trying to understand why.

EDIT : apparently the design created here:

Libra/R/pseudobulk_de.R

Line 84 in 388d41c

design = model.matrix(~ group, data = targets)

wrongfully uses the last level as the test group, hence why the logFC are all reversed. Not sure why they also differ in value (edgeR gives log2fc while Seurat::FoldChange used by Libra is on base exp(1) but it’s still different after making the correction).

Did anyone find a way to make sense of what Libra returns? Is it safer to use directly Seurat/Signac objects?

No judgment but is Libra still maintained? I see there's a lot of open issues (some of them I already encountered) and I wonder if it's worth suggesting pull requests or if I should look for an alternative.

[Random-x-1]

1. Change base from exp(1) to 2, and you can get consistent results to Seurat::FindMarkers with correct logfc and P value in singlecell and wilcox settings, also should change colname to make it consistent or you will get an error.

2. For edgeR, or other pseudobulk based methods, the design = model.matrix(~ group, data = targets) returned wrong design matrix and log2fc reversed, manually define reference level by relevel(targets$group, ref = levels(meta$label)[2]) can solve this problem.