I have been stuck a few days on this problem, so I thought I'd better ask. This is not just any edgeR question as the input data is different than usual (non-RNA). So if anyone knows of a better tool to run this analysis please tell me, otherwise just assume the data is RNA counts data.
Data: unnormalized counts (derived from depth i.e. mapped reads against a metagenomic bin). A metagenomic bin corresponds to an organism (bacterium species) so obviously its abundance is correlated with the abundance of other species within the same sample, but it may not be relevant to the question here.
Data structure: (every subject belongs to either C or T group and is sampled twice at t=0 and t=1)
df <- data.frame(timepoint = rep_len(0:1, length.out=16), subjectID = c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8), group = c(rep("C", 8), rep("T", 8)))
Aim: find differentially abundant metagenomic bins between the two cohorts (C vs T) i.e. which genomes go up/down between t=0 and t=1 in a different way between the two cohorts?
My code so far:
#clean dataframe by selecting rows with only a certain amount of NA values input_df <- df[rowSums( is.na(df ) <=80, ] #replace missing values with zeros input_df[is.na(input_df)] <- 0
Prepare the count data (convert countdata first column (genomes) to rownames):
count_data <- data.frame(input_df[,-1], row.names=input_df[,1], check.names = FALSE)
count_data_2 <- voom(count_data)
Prepare design matrix: (I am creating 4 groups: C_t0, C_t1, T_t0, T_t1)
Treat <- factor(paste(metadata$cohort,metadata$date,sep=".")) design <- model.matrix(~0+Treat) colnames(design) <- levels(Treat)
Block on subject: (gives me warnings)
corfit <- duplicateCorrelation(count_data_2,design,block=metadata$subject) corfit$consensus
Fit and make contrasts:
fit <- lmFit(count_data_2,design,block=metadata$subject,correlation=corfit$consensus) cm <- makeContrasts( T1vT0 = Treatment.1-Treatment.0, C1vC0 = Control.2-Control.1, T0vC0 = Treatment.1-Control.1, T1vC1 = (Treatment.1-Treatment.0)-(Control.2-Control.1), levels = design) fit2 <- contrasts.fit(fit, cm) fit2 <- eBayes(fit2) top20 <- topTable(fit2, n = 20, coef="T1vC1", adjust.method = "BH")
Q: Is this setup/my code correct? Can it be improved? Do I get in this way what I am after?
I also tried edgeR (instead of limma/voom) and I was getting seemingly more real p-values (as in: closer to p=0.005, whereas now with limma/voom I am getting suspiciously tiny p-values).
Thank you so much for your help!
Kind Regards Dany