const-ae / const-ae/glmGamPoi

Too many DEGs?

Open
#35 1 comment 0 reactions 0 assignees View on GitHub
Dominant language
R
Stars
126
Forks
18
PR merge metrics
No merged PRs in 30d

Description

Hi,

I am a fresh user of glmGamPoi and thanks for providing this awesome tool! I tried to identify DEGs in scRNA and successfully run the following code:
```
load('/Volumes/kfsd/GSE200903_new_combined_cluster_FibroSamples_2022-01-14.RData')
combine_cluster <- as.SingleCellExperiment(combined_cluster)

set.seed(1)
colData(combine_cluster)
# Take highly expressed genes and proper cells:
combine_cluster_subset <- combine_cluster[rowSums(counts(combine_cluster)) > 100,
which( ! is.na(combine_cluster$CellType) & combine_cluster$CellType
%in% c("fibroblasts"))]

# add a stim2 based on stim
colData(combine_cluster_subset)$stim2 <- colsplit(colData(combine_cluster_subset)$stim,
split = "_", names = c('cell', 'trt','time'))$trt
# Convert counts to dense matrix
counts(combine_cluster_subset) <- as.matrix(counts(combine_cluster_subset))
# Remove empty levels because glm_gp() will complain otherwise
colData(combine_cluster_subset)$stim2 <- as.factor(colData(combine_cluster_subset)$stim2)
colData(combine_cluster_subset)$stim2 <- relevel(colData(combine_cluster_subset)$stim2, ref = "normal")

fit <- glm_gp(combine_cluster_subset, design = ~ stim2 ,
reference_level = "normal")
summary(fit)
colnames(fit$Beta)

# The contrast argument specifies what we want to compare
# We test the expression difference of stimulated and control T-cells
#
# There is no sample label in the colData, so we create it on the fly
# from `stim` and `ind` columns in colData(fit$data).
de_res <- test_de(fit, contrast = `stim2PDAC`,
pseudobulk_by = paste0(stim2, "-", rownames(colData(combine_cluster_subset))))

# The large `lfc` values come from groups were nearly all counts are 0
# Setting them to Inf makes the plots look nicer
de_res$lfc <- ifelse(abs(de_res$lfc) > 20, sign(de_res$lfc) * Inf, de_res$lfc)

# Most different genes
de_res<-de_res[order(de_res$pval), ]
```
However, I identified 11752 DEGs out of 14537 genes and I found many genes has df1=1.
```
head(de_res)
name pval adj_pval f_statistic df1 df2 lfc
35 Rpl7 0 0 3244.666 1 20005.31 -0.5195890
44 Pi15 0 0 2180.603 1 20005.31 -1.6835958
124 Il1r1 0 0 3883.918 1 20005.31 2.2249821
132 Fhl2 0 0 5068.944 1 20005.31 4.8176423
158 Gls 0 0 1635.089 1 20005.31 0.9511179
164 1700019D03Rik 0 0 3380.809 1 20005.31 -2.7193859
```

I wondered how could I address this problem? Thanks in advance!

Best,
Kun

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.