[BioC] R: Time course experiment....
Manca Marco (PATH)
m.manca at maastrichtuniversity.nl
Wed Feb 9 17:19:36 CET 2011
Hi Viritha,
the first thing popping to my mind is that you could use the lfc condition in topTable (from limma's Reference Manual -> lfc specifies cutoff value for log2-fold-change. Only genes with larger fold changes are listed.)
Thus your code would look like
design <- model.matrix(~factor(rep(1:2, each=3)))
fit <- lmFit(eset, design)
fit2 <- eBayes(fit)
# If you want the first 10000 results... but usually I put "number=THE_NUMBER_OF_GENES_ON_MY_CHIP"
top.list<-topTable(fit2,number=10000,coef=2,lfc=1,adjust="BH")
your.results <- top.list[top.list$adj.P.Val < 0.05, ]
and you would have, at once, only genes whose logFC is higher than 1 or lower than -1 and with and adjusted p value of less than .05
In your code the mistake is this (look the line with the ash code for comment):
> your.results1<-your.results[your.results$logFC < -1,]
#you are selecting only genes with a logFC lower than -1 and saving them as your.results1
> your.results2<-your.results1[your.results1$logFC > 1,]
#now from the previously saved list of genes with logFC < -1 you are asking to sort out those that have a logFC > 1 ...since the previous step was executed correctly, this results in a 0 list (no gene with logFC<-1 can also have logFC>1 of course)
> write.table(your.results2,"Panc1_final_list.txt",sep='\t')
#well... you are saving an empty object...
I hope this helps...
You may find other functions of some interest for you, as decideTest, or options for topTable if you read the reference manual -> http://www.bioconductor.org/packages/release/bioc/manuals/limma/man/limma.pdf
All the best, Marco
--
Marco Manca, MD
University of Maastricht
Faculty of Health, Medicine and Life Sciences (FHML)
Cardiovascular Research Institute (CARIM)
Mailing address: PO Box 616, 6200 MD Maastricht (The Netherlands)
Visiting address: Experimental Vascular Pathology group, Dept of Pathology - Room5.08, Maastricht University Medical Center, P. Debyelaan 25, 6229 HX Maastricht
E-mail: m.manca at maastrichtuniversity.nl
Office telephone: +31(0)433874633
Personal mobile: +31(0)626441205
Twitter: @markomanka
*********************************************************************************************************************
This email and any files transmitted with it are confidential and solely for the use of the intended recipient.
It may contain material protected by privacy or attorney-client privilege. If you are not the intended recipient or the person responsible for
delivering to the intended recipient, be advised that you have received this email in error and that any use is STRICTLY PROHIBITED.
If you have received this email in error please notify us by telephone on +31626441205 Dr Marco MANCA
*********************************************************************************************************************
________________________________________
Da: viritha kaza [viritha.k at gmail.com]
Inviato: mercoledì 9 febbraio 2011 16.24
A: Manca Marco (PATH)
Cc: james.reid at ifom-ieo-campus.it; bioconductor at stat.math.ethz.ch
Oggetto: Re: [BioC] Time course experiment....
Hi Marco,
Thanks...
That happens.
I wanted to include one more condition on the top.table list of logFC to consider only those which are greater than or equal to1 and less than or equal to -1.
I just tried to expand the code
> your.results1<-your.results[your.results$logFC < -1,]
> your.results2<-your.results1[your.results1$logFC > 1,]
> write.table(your.results2,"Panc1_final_list.txt",sep='\t')
but just got a blank txt file with headings.I also tried using & and | operators, but got errors.
I hope you could help on this.
Thanks,
Viritha
On Tue, Feb 8, 2011 at 3:32 PM, Manca Marco (PATH) <m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl>> wrote:
Sorry, I had overlooked that part of the question... I thought you just wanted help in setting the threshold so I just added that line to your code.
As James has already mentioned to have the full list of genes in your results you have to set something like:
design <- model.matrix(~factor(rep(1:2, each=3)))
fit <- lmFit(eset, design)
fit2 <- eBayes(fit)
# If you want the first 10000 results... but usually I put "number=THE_NUMBER_OF_GENES_ON_MY_CHIP"
top.list<-topTable(fit2,number=10000,coef=2,adjust="BH")
your.results <- top.list[top.list$adj.P.Val < 0.05, ]
...I can only say to my defense "you only see what you look for" =P I should have given a more complete answer since the beginning.
All the best, Marco
--
Marco Manca, MD
University of Maastricht
Faculty of Health, Medicine and Life Sciences (FHML)
Cardiovascular Research Institute (CARIM)
Mailing address: PO Box 616, 6200 MD Maastricht (The Netherlands)
Visiting address: Experimental Vascular Pathology group, Dept of Pathology - Room5.08, Maastricht University Medical Center, P. Debyelaan 25, 6229 HX Maastricht
E-mail: m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl>
Office telephone: +31(0)433874633
Personal mobile: +31(0)626441205
Twitter: @markomanka
*********************************************************************************************************************
This email and any files transmitted with it are confidential and solely for the use of the intended recipient.
It may contain material protected by privacy or attorney-client privilege. If you are not the intended recipient or the person responsible for
delivering to the intended recipient, be advised that you have received this email in error and that any use is STRICTLY PROHIBITED.
If you have received this email in error please notify us by telephone on +31626441205 Dr Marco MANCA
*********************************************************************************************************************
________________________________________
Da: viritha kaza [viritha.k at gmail.com<mailto:viritha.k at gmail.com>]
Inviato: martedì 8 febbraio 2011 19.01
A: Manca Marco (PATH)
Cc: bioconductor at stat.math.ethz.ch<mailto:bioconductor at stat.math.ethz.ch>
Oggetto: Re: [BioC] Time course experiment....
Hi Marco,
Thanks for your reply.
But topTable by default only gives top 10.So how do I know the number in topTable fuction which satisfies the criteria of 0.05 that I need to give.
waiting for your reply...
Thanks,
Viritha
On Tue, Feb 8, 2011 at 12:35 PM, Manca Marco (PATH) <m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl><mailto:m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl>>> wrote:
Hi Viritha,
just do
code:
design <- model.matrix(~factor(rep(1:2, each=3)))
fit <- lmFit(eset, design)
fit2 <- eBayes(fit)
top.list<-topTable(fit2,coef=2,adjust="BH")
your.results <- top.list[top.list$adj.P.Val < 0.05, ]
...it should do the trick...
All the best, Marco
--
Marco Manca, MD
University of Maastricht
Faculty of Health, Medicine and Life Sciences (FHML)
Cardiovascular Research Institute (CARIM)
Mailing address: PO Box 616, 6200 MD Maastricht (The Netherlands)
Visiting address: Experimental Vascular Pathology group, Dept of Pathology - Room5.08, Maastricht University Medical Center, P. Debyelaan 25, 6229 HX Maastricht
E-mail: m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl><mailto:m.manca at maastrichtuniversity.nl<mailto:m.manca at maastrichtuniversity.nl>>
Office telephone: +31(0)433874633
Personal mobile: +31(0)626441205
Twitter: @markomanka
*********************************************************************************************************************
This email and any files transmitted with it are confidential and solely for the use of the intended recipient.
It may contain material protected by privacy or attorney-client privilege. If you are not the intended recipient or the person responsible for
delivering to the intended recipient, be advised that you have received this email in error and that any use is STRICTLY PROHIBITED.
If you have received this email in error please notify us by telephone on +31626441205 Dr Marco MANCA
*********************************************************************************************************************
________________________________________
Da: bioconductor-bounces at r-project.org<mailto:bioconductor-bounces at r-project.org><mailto:bioconductor-bounces at r-project.org<mailto:bioconductor-bounces at r-project.org>> [bioconductor-bounces at r-project.org<mailto:bioconductor-bounces at r-project.org><mailto:bioconductor-bounces at r-project.org<mailto:bioconductor-bounces at r-project.org>>] per conto di viritha [viritha.k at gmail.com<mailto:viritha.k at gmail.com><mailto:viritha.k at gmail.com<mailto:viritha.k at gmail.com>>]
Inviato: martedì 8 febbraio 2011 18.14
A: bioconductor at stat.math.ethz.ch<mailto:bioconductor at stat.math.ethz.ch><mailto:bioconductor at stat.math.ethz.ch<mailto:bioconductor at stat.math.ethz.ch>>
Oggetto: Re: [BioC] Time course experiment....
Naomi Altman <naomi at ...> writes:
>
> I did not get the original posting, but doesn't Sohail just need
> "TopTable" for this?
>
> --Naomi
Hi Naomi,
I had a similiar issue and used TopTable but I am getting only first 10 but how
should I write the code to get selected probesets based on an adjusted p-value
less than or equal to 0.05.
code:
design <- model.matrix(~factor(rep(1:2, each=3)))
fit <- lmFit(eset, design)
fit2 <- eBayes(fit)
topTable(fit2,coef=2,adjust="BH")
I have 3 control and 3 treated samples.The eset has the probenames and the
expression matrix.
So my question is how do I get the data of all the probesets which pass this
criteria?
Thanks,
Viritha
>
> At 08:51 AM 1/2/2006, Gordon Smyth wrote:
> >Dear Sohail,
> >
> >Well, there are lots of ways to generate such a table. Perhaps the simplest
is
> >
> > fitsel <- fit2[sel.dif, ]
> > as.data.frame( fitsel )
> >
> >Best wishes
> >Gordon
> >
> > >Date: Tue, 20 Dec 2005 14:03:43 -0500
> > >From: "Khan, Sohail" <khan at ...>
> > >Subject: [BioC] Time course experiment....
> > >To: <bioconductor at ...>
> > >
> > >Dear List,
> > >
> > >I have performed a time course analysis using limma, as described in
> > >"Bioinformatics and Computational Biology Solutions ........".
> > >How can I get a list of differentially expressed genes? I've tried
> > >the code below:
> > >sel.dif<-p.adjust(fit2$F.p.vlaue,method="fdr") <0.05
> > >This produces a logical vector, right?. I would like a table of
> > >differentially expressed genes with p vales etc. Sorry, if I missed
> > >this in the limma user's guide. Thanks for any suggestions.
> > >
> > >
> > >Sohail Khan
> > >Scientific Programmer
> > >COLD SPRING HARBOR LABORATORY
> > >Genome Research Center
> > >500 Sunnyside Boulevard
> > >Woodbury, NY 11797
> > >(516)422-4076
> >
> >_______________________________________________
> >Bioconductor mailing list
> >Bioconductor at ...
> >https://stat.ethz.ch/mailman/listinfo/bioconductor
>
> Naomi S. Altman 814-865-3791 (voice)
> Associate Professor
> Dept. of Statistics 814-863-7114 (fax)
> Penn State University 814-865-1348 (Statistics)
> University Park, PA 16802-2111
>
_______________________________________________
Bioconductor mailing list
Bioconductor at r-project.org<mailto:Bioconductor at r-project.org><mailto:Bioconductor at r-project.org<mailto:Bioconductor at r-project.org>>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list