write function or forloop for contrasts and MA plots

raynamharris commented 5 years ago

I would like to construct contrasts for edgeR using a function or for loop. I don't yet know how to do this, yet....

The problem: Currently, I am constructing contrasts to make MA plots by hand in 01_limma_all.Rmd based on the levels in colData$group. For example, here are 6 contrasts comparing hatch and lay stages for a given sex and tissue.

my.contrasts <- makeContrasts(
             FG_HL = female.gonad.hatch - female.gonad.lay,
             FH_HL = female.hypothalamus.hatch - female.hypothalamus.lay,
             FP_HL = male.pituitary.hatch - male.pituitary.lay,
             MP_HL = male.pituitary.hatch - male.pituitary.lay,
             MH_HL = male.hypothalamus.hatch - male.hypothalamus.lay,
             MG_HL = male.gonad.hatch - male.gonad.lay,          
levels=parentaldesign)

Hand-typing the contrasts is a problem because:

There are 74 levels, so it would take forever to make all the sex and tissue and treatment contrasts necessary to explore all the two-way comparisons.
More importantly, this is typo prone! Notice I wrote FP_HL = male.pituitary.hatch - male.pituitary.lay... So my third figure in 01_limma_all.md is miss-labelled female when it should be male.

raynamharris commented 5 years ago

I haven't sorted out how to use a loop or function to create my.contrasts.

But, I did write a function and loop to create plots for all the contrasts! For the above example it would look like this:

mycontrasts <- c("FG_HL", "FH_HL", "FP_HL", "MP_HL", "MH_HL", "MG_HL")

printplotcontrasts <- function(whichcontrast){
  cont <- whichcontrast
  summary(decideTestsDGE(
    glmTreat(fit, contrast=my.contrasts[,cont], lfc = 1), 
    adjust.method="fdr", p.value=0.01))
  print(kable(topTags(glmTreat(fit, contrast=my.contrasts[,cont]), n=5), digits=2, lfc = 1))
  print(plotMD(glmTreat(fit, contrast=my.contrasts[,cont], lfc=1), main=whichcontrast, frame.plot=F))
}

for(i in mycontrasts){
  printplotcontrasts(i)
}

This is implemented in

raynamharris commented 5 years ago

Correction: print and kable don't seem to agree, and I forgot to print the summary.

printplotcontrasts <- function(whichcontrast){
  cont <- whichcontrast
  print(summary(decideTestsDGE(
    glmTreat(fit, contrast=my.contrasts[,cont], lfc = 1), 
    adjust.method="fdr", p.value=0.01)))
  print(topTags(glmTreat(fit, contrast=my.contrasts[,cont]), n=5), digits=2, lfc = 1)
  print(plotMD(glmTreat(fit, contrast=my.contrasts[,cont], lfc=1), main=whichcontrast, frame.plot=F))
}

raynamharris commented 4 years ago

Currently have a R script that goes through all the contrast that are related to the aims of the project

https://github.com/macmanes-lab/DoveParentsRNAseq/blob/master/analysis/02_DESeq2.R

macmanes-lab / DoveParentsRNAseq

write function or forloop for contrasts and MA plots #3