I am trying to analyze some data with the tidybulk package. This tutorial is great

I have a paired design as described in section 3.4.1 in the edgeR guide. I have 4 independent subjects and each has a pre-treatment (control) and post-treatment sample. I want to see what genes were differential expressed in each subject after treatment.

My base R code is as such:

library(edgeR)
data=read.csv("~/Rmd_docs/transcriptome_counts.csv")
data=mutate(data, new_sample_name=paste(subjectID, ifelse(Visit=='Screening',"Pre","Post"),sep="-"))

samplenames <- c("115-20407-Post", "115-20407-Pre",
                 "115-20413-Post", "115-20413-Pre",
                 "126-20408-Post", "126-20408-Pre", 
                 "142-20104-Post", "142-20104-Pre")
subjectids <- c("115-20407", "115-20407",
                "115-20413", "115-20413",
                "126-20408", "126-20408",
                "142-20104", "142-20104")
treatments <- c("T","C",
                "T", "C",
                "T", "C",
                "T","C")
targets <- data.frame (samplenames, subjectids, treatments)

Subject <- factor(targets$subjectids)
Treat <- factor(targets$treatments, levels=c("C","T"))
design <- model.matrix(~ Subject + Treat)

data_mat=dcast(data, gene~ new_sample_name, value.var = "raw_count")

y <- DGEList(counts=data_mat[,2:9], genes=data_mat[,1], group = Treat)

#filter and normalize
keep <- filterByExpr(y, design)
y <- y[keep,,keep.lib.sizes=FALSE]


y_normed <- calcNormFactors(y)
rownames(design) <- colnames(y_normed)
y_normed <- estimateDisp(y_normed, design, robust=TRUE)

fit <- glmQLFit(y_normed, design)
qlf <- glmQLFTest(fit)
topTags(qlf)
# include your problematic code here with any corresponding output 
# please also include the results of running the following in an R session 

sessionInfo( )

For my tidybulk code, I have a file with sampleID, subjectID, gene, raw count, and treated (logical column)

library(tidybulk)
tb_counts=tidybulk(counts_indication_df, .sample = sampleID, .transcript = gene, .abundance = raw_count)
counts_scaled_df <- tb_counts %>% 
    keep_abundant(factor_of_interest = Treated) %>%
    scale_abundance()

counts_de <- counts_scaled_df %>%
  test_differential_abundance( ~ subjectID + Treated) # is this correct for paired-design

But when I compare p-values from tidybulk edgeR with my base R edgeR results, they are not any where close to being the same:

I must be misunderstanding how to correctly specify the design matrix for paired experiments. I would appreciate any help.

Hello Amit,

few things I would check

• do you expect the distribution of p.value edgeR to be flat? (top-left quadrant)?
• when you compare edgeR with tidybulk wrapper, is the gene ordering the same?

A feedback on the code

• scale_abundance() is not needed/used for test_differential_abundance()
• with tidybulk, you don't need to create multiple variables, something like this

counts_de <- 
    tb_counts %>% 
    keep_abundant(factor_of_interest = Treated) %>%
    scale_abundance() %>% # Not needed
    test_differential_abundance( ~ ...)

does the job

[EDIT: PLEASE IGNORE BELOW COMMENT: see reply below]

A comment

if you don't use contrasts, usually the formula has the factor of interest as the first covariate

~ Treated + subjectID

[EDIT: PLEASE IGNORE ABOVE COMMENT: see reply below]