Transcriptomic Differential analysis of mouse mammary gland RNA-Seq

In ViSEAGO: ViSEAGO: a Bioconductor package for clustering biological functions using Gene Ontology and semantic similarity

###################
# knitr document options
knitr::opts_chunk$set(eval=TRUE,fig.pos = 'H',
fig.width=9,message=FALSE,comment=NA,warning=FALSE)

###################
# library
library(edgeR)
library(limma)
library(graphics)

###################
# Get Raw Count data
###################

  ###################
  # temp file
  temp<-paste(
    tempfile(),
    "gz",
    sep="."
  )

  ###################
  # load the file
  download.file(
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60450/suppl/GSE60450_Lactation-GenewiseCounts.txt.gz",
    destfile =temp,
    quiet=TRUE
  )

  ###################
  # unzip
  gunzip(temp)

  ###################
  # load raw data
  rawdata<-data.table::fread(sub("\\.gz","",temp))

  ###################
  # genecount
  geneCount <-as.matrix(rawdata[,-c(1,2),with=FALSE])

  ###################
  # add rownames
  rownames(geneCount)<-rawdata$EntrezGeneID

  ###################
  # design and factors
  targets<-data.table::data.table(
    File=c(
      "GSM1480297_MCL1-DG_BC2CTUACXX_ACTTGA_L002_R1.txt",
      "GSM1480298_MCL1-DH_BC2CTUACXX_CAGATC_L002_R1.txt",
      "GSM1480299_MCL1-DI_BC2CTUACXX_ACAGTG_L002_R1.txt",
      "GSM1480300_MCL1-DJ_BC2CTUACXX_CGATGT_L002_R1.txt",
      "GSM1480301_MCL1-DK_BC2CTUACXX_TTAGGC_L002_R1.txt",
      "GSM1480302_MCL1-DL_BC2CTUACXX_ATCACG_L002_R1.txt",
      "GSM1480291_MCL1-LA_BC2CTUACXX_GATCAG_L001_R1.txt",
      "GSM1480292_MCL1-LB_BC2CTUACXX_TGACCA_L001_R1.txt",
      "GSM1480293_MCL1-LC_BC2CTUACXX_GCCAAT_L001_R1.txt",
      "GSM1480294_MCL1-LD_BC2CTUACXX_GGCTAC_L001_R1.txt",
      "GSM1480295_MCL1-LE_BC2CTUACXX_TAGCTT_L001_R1.txt",
      "GSM1480296_MCL1-LF_BC2CTUACXX_CTTGTA_L001_R1.txt"
    ),
    Sample=c(
      "GSM1480297",
      "GSM1480298",
      "GSM1480299",
      "GSM1480300",
      "GSM1480301",
      "GSM1480302",
      "GSM1480291",
      "GSM1480292",
      "GSM1480293",
      "GSM1480294",
      "GSM1480295",
      "GSM1480296"
    ),
    CellType=rep(c("Basal","Luminal"),each=6),
    Status=rep(rep(c("virgin","pregnant","lactate"),2),each=2)
  )

###################
# design and factors
group <- factor(
  paste0(
    targets$CellType, ".", targets$Status)
)

###################
# create DGElist
y <-DGEList(geneCount, group=group)

###################
# minima count filtering
keep <- rowSums(cpm(y) > 0.5) >= 2

###################
# filter y
y <- y[keep, , keep.lib.sizes=FALSE]

###################
# Normalization
y <-calcNormFactors(y)

###################
# graphic parameters
colors<-rep(c("blue", "darkgreen", "red"),2)
points<-c(0,1,2,15,16,17)

###################
# MDSplot
plotMDS(
  y,
  col=colors[group],
  pch=points[group]
)

###################
# add MDSplot legend
legend(
  "topleft",
  legend=levels(group),
  pch=points,
  col=colors,
  ncol=2
)

###################
# model
design <-stats::model.matrix(~ 0 + group)

###################
# renames
colnames(design) <- levels(group)

###################
# Estimate Common, Trended and Tagwise Negative Binomial dispersions by weighted likelihood empirical Bayes robustified against outliers
y <- estimateDisp(y, design, robust=TRUE)
y$common.dispersion

###################
# Plot Biological Coefficient of Variation
plotBCV(y)

###################
# Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood Tests
fit <-glmQLFit(y, design, robust=TRUE)

###################
# Plot the quasi-likelihood dispersion
plotQLDisp(fit)

###################
# Construct Matrix of Custom Contrasts
con<-makeContrasts(
  L.PvsL = Luminal.pregnant - Luminal.lactate,
  L.VvsL = Luminal.virgin - Luminal.lactate,
  L.VvsP = Luminal.virgin - Luminal.pregnant,
  levels=design
)

###################
# Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood Tests
qlf2<-glmQLFTest(fit, contrast=con)

###################
# Construct Matrix of Custom Contrasts
con.L.pregnantvslactate <- makeContrasts(
  L.PvsL = Luminal.pregnant - Luminal.lactate,
  levels=design
)

con.L.virginvslactate <-makeContrasts(
  L.VvsL = Luminal.virgin - Luminal.lactate,
  levels=design
)

con.L.virginvspregnant <-makeContrasts(
  L.VvsP = Luminal.virgin - Luminal.pregnant,
  levels=design
)

###################
# Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood Tests
res<-lapply(ls(pattern="con."),function(x){

  ###################
  # Genewise Negative Binomial Generalized Linear Models with Quasi-likelihood Tests
  qlf<-glmQLFTest(
    fit,
    contrast=get(x)
  )

  ###################
  # Table of the Top Differentially Expressed Genes/Tags
   topTags(
     qlf,
     n=nrow(y$counts),
     adjust.method="BH",
     sort.by="none"
    )$table
})
names(res)<-gsub("^con\\.","",ls(pattern="con."))

###################
# exporting results of the differential analysis
write(
  rownames(res$L.pregnantvslactate),
  file="background_L.txt"
)

write(
  rownames(res$L.pregnantvslactate)[res$L.pregnantvslactate$FDR<0.05],
  file="pregnantvslactateDE.txt"
)

write(
  rownames(res$L.virginvslactate)[res$L.virginvslactate$FDR<0.05],
  file="virginvslactateDE.txt"
)

write(
  rownames(res$L.virginvspregnant)[res$L.virginvspregnant$FDR<0.05],
  file="virginvspregnantDE.txt"
)

ViSEAGO documentation built on Nov. 8, 2020, 6:51 p.m.