Hi everyone,
I am new to microarray data and I am having a really hard time finding a good work flow so far. My experiment is with human transcriptome data. I have 4 time points and I am trying to find DEG between them. This is just a test run with two times points. I have no idea what to do and so far my script hasn't worked so well. This is what I have so far. Thanks in advance. My volcano plot also looks really weird.
source("http://bioconductor.org/biocLite.R") biocLite("oligo") biocLite("limma") biocLite("pd.hugene.2.1.st") library(genefilter) library(limma) library(oligo) library(pd.hugene.2.0.st) #directory for the data mydir <- "C:\\Users\\hakim\\Desktop\\Bioinformatics_Thesis_Concordia_Microarry_Data\\Control" #setting seed for reproducibility set.seed(1) #listing the files from directory using special CEL file read function celList <- list.celfiles(mydir, full.names=TRUE) #reading data from cellist and setting annotation package to approiate one for this microarray rawData <- read.celfiles(celList, pkgname='pd.hugene.2.0.st') #normalizing the data using RMA algorithm normData <- rma(rawData) #checking boxplot of raw data par(mar=c(10,4.5,2,1)) boxplot(rawData,las=3) #checking boxplot of normalized data boxplot(normData,las=3) #the annotation package biocLite("hugene20sttranscriptcluster.db") library(hugene20sttranscriptcluster.db) #retreaving feature data featureData(normData) <- getNetAffx(normData, "transcript") design = cbind(control = 1, controlvstreatment = c((rep.int(0,13)),rep.int(1,13))) fit<- lmFit(normData, design) normfit <-eBayes(fit) omg <- topTable(normfit, coef="controlvstreatment", adjust="BH") write.table(omg,file="finally.txt",sep= "/" ) results <- decideTests(fit) vennDiagram(results) volcanoplot(normfit)
What do you mean by "your script hasn't worked so well"? Has it crashed, or thrown an error? If you don't detect any DE genes, that doesn't mean that the script has failed; it could just mean that you've got no DE. You need to give more information here, as well as something more useful regarding the "weirdness" of the volcano plot.
The script works, I should have worded that differently. I am the one who doesn't understand so well the output of my script.
My first question is if it is possible to use the limma package with dependent data (time series). If so then this script may be alright; I still do not understand my output so much.
This is the summary of my results.
control controlvstreatment
-1 0 474
0 0 52791
1 53617 352
This is my volcano plots and venn diagram. I am not really sure what the volcano plot is telling me and if this is normal.
http://imgur.com/a/CStnP
Thanks for the help!