Hi there!

I am new at transcriptomic analysis and at posting questions (I usually manage by reading other peoples questions), so let me know if something is unclear (or way too long). I am trying to use edgeR to analyse my transcriptomic data, but I'm running into problems both with regards to design matrix (I can get it to work, but I don't understand it) and as a consequence, I also get problems with my contrasts.

My transcripts come from a bacteria-phage model, and what I am looking at now is the bacterial genes (i.e. how does the phage change the expression of host genes). I have also done the experiment in two different treatments, "full" and "low" (it reflects what kind of medium the bacteria were grown in), and I've sampled each set 6 times. In other words, I have: 3 x bacteria+phage in Full medium for times 0-5 (PXF, where x is a time in minutes - determined by phage progression) 3 x bacteria+water in Full medium for times 0-5 (CXF, x is the same as above, water is just as a control) 3 x bacteria+phage in Low medium for times 0-5 (PYL, where y is a time in minutes - determined by phage progression) 3 x bacteria+water in Low medium for times 0-5 (CYL, y is the same as above) (The values for x and y are not the same, but they are the same proportions of the phages progression, i.e. supposed to be "translatable".)

It totals to 72 samples. It is also important to know that the samples from Full medium comes from 3 sets of starter cultures, i.e. I grew bacteria in a big bottle over night and then I split it in two and treated one part with phage and one with water = pairs. This is the same for the samples in Low medium, they are basically 3 pairs. And to complicate it a bit further, after some crude visualisation of mapping against phage/bacteria genomes - I discovered that three successive (i.e. from time 3, time 4 and time 5 of the same culture) control samples had too much phage genes and therefore had been contaminated, and that one phage sample did not contain any phage genes when that is biologically impossible. Therefore I need to remove these four samples. Is it possible to keep the "corresponding" samples, or do I have to discard them as well? In the analysis I describe here, I've only removed these four samples (which I am very uncertain whether it's correct or not).

I have several hypotheses about the experiment, but first of all I want to compare phage infected bacteria at each time point against the control, for each kind of medium. Secondly, I want to compare the difference between the phage vs control between the two kinds of medium.

To achieve this I have loaded my data and produced a data.frame with sample information:

#reading the cleaned, mapped and "counted" data into R
count <- read_tsv("counts.tsv", col_types = cols(.default = col_character(), count = col_double()))
#parsing the count table to produce metadata for the samples in count
samples <- count %>% distinct(sample) %>%
  mutate(
    treatment = sub('^(.).*', '\\1', sample),
    time = sub('.*t([0-9]+)$', '\\1', sample),
    replicate = sub('.*([0-9])+[LF]t[0-9]+', '\\1', sample),
    medium = sub('.*[0-9]+([LF])t.*', '\\1', sample),
    time = case_when(
      medium == 'L' & time == '0' ~ 1,
      medium == 'L' & time == '1' ~ 16.25,
      medium == 'L' & time == '2' ~ 32.5,
      medium == 'L' & time == '3' ~ 48.75,
      medium == 'L' & time == '4' ~ 65,
      medium == 'L' & time == '5' ~ 113.75,
      medium == 'F' & time == '0' ~ 1,
      medium == 'F' & time == '1' ~ 12.5,
      medium == 'F' & time == '2' ~ 25,
      medium == 'F' & time == '3' ~ 37.5,
      medium == 'F' & time == '4' ~ 50,
      medium == 'F' & time == '5' ~ 87.5,
      TRUE                        ~ -1
     )
  )
#merging three coloumns from data set "samples" - in other words reducing the samples information to be a guide with sample and characteristics
united <- unite(samples, "tr_ti_med", c("treatment", "time", "medium"))
#sorting
united <- united[order(united$sample),]
head(united)
sample  tr_ti_med   replicate
C1Ft0   C_1_F       1       
C1Ft1   C_12.5_F    1       
C1Ft2   C_25_F      1       
C1Ft3   C_37.5_F    1       
C1Ft4   C_50_F      1       
C1Ft5   C_87.5_F    1

And so far, I think I can handle it. Also, creating a DGEList seems fine, I've done the following:

#create DGEList
count_wide <- count %>%
  spread(sample, count, fill = 0)
count_matrix <- as.matrix(count_wide %>% select(-geneid))
rownames(count_matrix) = count_wide$geneid
dgel_count <- DGEList(count_matrix, group=united$tr_ti_med)
#subset for only bacteria genes (phage genes start with "B", no bacterial genes does)
dgel_bact <- dgel_count[grep("B.*",rownames(dgel_count$counts), invert=TRUE), ]
#subsetting bacterial dataset to remove C3Lt3-5 (contaminated), and also P2Ft3 (biologically invalid)
dgel_bact <- dgel_bact[, grep("C3Lt[3-5]", rownames(dgel_bact$samples), invert=TRUE)]
dgel_bact <- dgel_bact[, grep("P2Ft3", rownames(dgel_bact$samples), invert=TRUE)]
#filtering out low counts, then checking the output, applying and calculate nomarlistaion factors
keep_b <- rowSums(cpm(dgel_bact) > 9) >= 3
dgel_bact <- dgel_bact[keep_b, , keep.lib.size=FALSE]
dgel_bact <- calcNormFactors(dgel_bact)

However, when I construct the design matrix I get a bit lost. I've been reading the edgeR's users guide and since I've got paired samples I've been looking at section 3.4.1 and decided I want to create the matrix with an additive formula, and to make it easier for me I'd like to exclude the intercept, I've therefore done the following:

treat <- as.factor(united$tr_ti_med) #contains information on phage/control, time and medium type
subject <- as.factor(united$replicate) #which subject the samples originate from

design <- model.matrix(
  ~0+subject+treat
)
#design needs to contain all samples of interest, therefore subsetting for bacteria,
#I've removed C3Lt3 (row 34), C3Lt4 (row 35), C3Lt5 (row 36) and P2Ft3 (row 52)
design_bact <- design[-c(34:36,52),]
design_bact

My first issue with the design matrix is that I'm lying. In my first attempts, I had split "subject/replicates" into 1F, 2F, 3F, 1L, 2L and 3L to indicate first/second/third pair in Full and in Low medium. However, I ran into "matrix not of full rank" and after trying to understand I read (e.g. here) that I can not have columns that add up to each other because of linear dependencies. Therefore I removed the F/L, and only use "1", "2", "3" - even though the "1" replicates in Full medium have nothing to do with the samples in Low medium. Is this an okay way of getting around the linear dependencies?

My second issue with the design matrix is what columns I get, e.g.:

> colnames(design_bact)
 [1] "subject1"        "subject2"        "subject3"        "treatC_1_L"      "treatC_113.75_L" "treatC_12.5_F"  
 [7] "treatC_16.25_L"  "treatC_25_F"     "treatC_32.5_L"   "treatC_37.5_F"   "treatC_48.75_L"  "treatC_50_F"    
[13] "treatC_65_L"     "treatC_87.5_F"   "treatP_1_F"      "treatP_1_L"      "treatP_113.75_L" "treatP_12.5_F"  
[19] "treatP_16.25_L"  "treatP_25_F"     "treatP_32.5_L"   "treatP_37.5_F"   "treatP_48.75_L"  "treatP_50_F"    
[25] "treatP_65_L"     "treatP_87.5_F"

Since I specified that I did not want an intercept, I get the three subject-columns as I want. But for the treatment, I'm missing "treatC1F". Here, I'm properly lost. For some reason, I think the design assumes that the first "factor" (not sure if that's the right word, here it's "treatC1F") is some kind of baseline, which in my case is not true. I do not have one treatment that could be baseline for all. I have tried forcing a matrix with all my levels, but then I get issues with "design matrix not of full rank". After a while I tried to accept that this is the correct way of doing the matrix and moved on to contrasts, which leads to my issues with contrasts. But I would like to know if I have set up my matrix in a proper way? And could someone try to explain why one of my levels is missing? And what that means? I can understand if that is too big to manage, but maybe there is some not too advanced literature that you could point me to (I'm not very good with statistics, I'm sad to say).

If I assume that the design matrix is fine, I then want to do some actual analysis and I've tried the following:

#getting dispersion
dgel_bact <- estimateDisp(dgel_bact, design_bact, robust=TRUE)
#the data should be fitted with quasi-likelihood
fit_bact <- glmQLFit(dgel_bact, design_bact, robust=TRUE)
#a contrast that works
lowPvsC_t0 <- makeContrasts(treatP_1_L-treatC_1_L, levels=design_bact)
res <- glmQLFTest(fit_bact, contrast=lowPvsC_t0)
topTags(res)
is.de <- decideTestsDGE(res)
summaryis.de)
#the contrast that does not work
fullPvsC_t0 <- makeContrasts(treatP_1_F-treatC_1_F, levels=design_bact)

So, regarding my problems with contrasts, I can compare phage affected bacteria vs control in the Low medium at whatever timepoint I want. But I cannot compare anything to the control bacteria at time 0 (which is "1" minute, therefore the "1" above) in Full medium. If I've done the design matrix in a proper way, how do I manage to "pull" out the information from samples that belong to "treatC1F"?

This post is potentially too long to keep anyones interest, but I hope not, and I would appreciate any help. I've put my questions in bold in an attempt to simplify the reading, but please let me know if I should shorten this.

Best wishes, Emelie.

EDIT: After finding someone I could discuss my experimental setup with, we realised that my "subjects" (the pairs/replicates) are 100% connected to the medium (L/F) that I used (it wouldn't be biologically logical to do it any other way), which means (if I've understood it correctly) that they are linear dependent and that I cannot include both because the model will be unable to decide which causes the variation. (I might get the words wrong, I apologise if so.)

Emelie, I don't have time to check your experimental design or to respond to all your questions, but I think your questions are probably all FAQ and you will find info on this forum if you do a Google search. I will say that if you simply reverse the factors in the linear model:

design <- model.matrix(~0+treat+subject)

then all your contrasts will work.

There is also no problem with the standard design matrix

design <- model.matrix(~ subject + treat)

In that case, you simply use topTags(fit, coef="treatP_1_F") to get the P_1_F vs C_1_F contrast, because all the treat coefficients are relative to C_1_F.