Help with metatranscriptomics: makeContrasts design and checking for housekeeping gene expression levels?

Question

[pearl2070]

I have 2 sets of questions about a metatranscriptomics analysis I'm working on.

(1) My experimental design was to compare OG environmental samples (samples kept 30 mins at original temperature) to HT samples (environmental samples kept high temperature for 30 mins), where both OG and HT samples were in glass bottles. I had a CN/control sample group that was taken at the original temperature, immediately prior to bottling, to be able to tell how much of an effect being bottled had on the samples. RNA-seq was done on all the samples and I have metatranscriptomic data. Unfortunately, it looks like all 3 sample groups (OG, HT, CN) have differential expression, where I was hoping that OG and CN would have similar expression. Is there a way to tease apart what the effect of the temperature treatment was, with something like the following setup?

#limma's makeContrasts
glm.contrasts<-makeContrasts(
  TreatmentEffect= HT- (OG-CN),
  levels=design.mat)

(2) Is there a way for me to look at the expression level of housekeeping genes in the metatranscriptome samples so that I can see if they look consistent between samples? I expect these metatranscriptomes to be dominated by cyanobacteria, which are my organism of interest in the samples, so I tried searching up housekeeping genes in cyanobacteria and found 3 that are used in PCR normalization. Are those appropriate to use?

Here is how I tried to assess if their expression was consistent between some of the samples:

list<-(grepl(c("RNPA|SECA|PPC"), rownames(dds[["counts"]])))

ddsSmall <- dds[list, ]

fit.glm2 <- glmFit(ddsSmall, design.mat)

glm.contrasts2<-makeContrasts(
  OGvHT= OG-HT,
  levels=design.mat)

GS.allvSB.glm2 <- glmLRT(fit.glm2, contrast=glm.contrasts[,"OGvHT"])

plotMD(OGvHT)

The plot reveals that some of these "housekeeping genes" are actually upregulated. But I think the filtering went wrong? The rownames of the dds counts tables are from a Trinity assembled metatranscriptome with annotations from Trinnotate and have IDs like TRINITY_DN4765_c2_g1 and TRINITY_DN20080_c0_g1^SP:SECA_ROSS1^sigP. The filtered dds1 also include some rownames which don't have any of the 3 genes I've filtered for. Is there a better annotation tool for Trinity assembled metatranscriptomes and/or a better way to filter for specific genes in the data?

[Ram]

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. If your code has long lines with a single command, break those lines into multiple lines with proper escape sequences so they're easier to read and still run when copy-pasted. I've done it for you this time.