Hi James, Am I still making any silly mistake over here. I am comparing normal (21) with lesionous patient samples (34). This is what I get in R console. Is there a problem with the script? Any advice would greatly help for further analysis.Thanks, Roopa > getwd() [1] "C:/Documents and Settings/rsubbaiaih/My Documents" > setwd(dir="/CRSP 406-11/DEMOS/GSE14905-a") > ls() character(0) > #-----------------------------------------------# > library(affy) > eset = justRMA() Loading required package: AnnotationDbi > eset ExpressionSet (storageMode: lockedEnvironment) assayData: 54675 features, 55 samples element names: exprs, se.exprs protocolData sampleNames: GSM372286.CEL GSM372287.CEL ... GSM372367.CEL (55 total) varLabels: ScanDate varMetadata: labelDescription phenoData sampleNames: GSM372286.CEL GSM372287.CEL ... GSM372367.CEL (55 total) varLabels: sample varMetadata: labelDescription featureData: none experimentData: use 'experimentData(object)' Annotation: hgu133plus2 > f <- factor(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, + 2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2), labels=c("Healthy", "affected")) > design <- model.matrix(~ 0 + f) > design fHealthy faffected 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 6 1 0 7 1 0 8 1 0 9 1 0 10 1 0 11 1 0 12 1 0 13 1 0 14 1 0 15 1 0 16 1 0 17 1 0 18 1 0 19 1 0 20 1 0 21 1 0 22 0 1 23 0 1 24 0 1 25 0 1 26 0 1 27 0 1 28 0 1 29 0 1 30 0 1 31 0 1 32 0 1 33 0 1 34 0 1 35 0 1 36 0 1 37 0 1 38 0 1 39 0 1 40 0 1 41 0 1 42 0 1 43 0 1 44 0 1 45 0 1 46 0 1 47 0 1 48 0 1 49 0 1 50 0 1 51 0 1 52 0 1 53 0 1 54 0 1 55 0 1 attr(,"assign") [1] 1 1 attr(,"contrasts") attr(,"contrasts")$f [1] "contr.treatment" > colnames(design) <-c("Healthy","affected") > library(limma) > fit <- lmFit(eset, design) > fit An object of class "MArrayLM" $coefficients Healthy affected 1007_s_at 10.081309 9.548286 1053_at 6.807501 7.482849 117_at 5.969921 6.147594 121_at 7.403842 7.733666 1255_g_at 2.804475 2.827041 54670 more rows ... $rank [1] 2 $assign [1] 1 1 $qr $qr Healthy affected 1 -4.5825757 0.000000 2 0.2182179 -5.830952 3 0.2182179 0.000000 4 0.2182179 0.000000 5 0.2182179 0.000000 50 more rows ... $qraux [1] 1.218218 1.000000 $pivot [1] 1 2 $tol [1] 1e-07 $rank [1] 2 $df.residual [1] 53 53 53 53 53 54670 more elements ... $sigma 1007_s_at 1053_at 117_at 121_at 1255_g_at 0.2866489 0.4801618 0.3499880 0.2332555 0.1053397 54670 more elements ... $cov.coefficients Healthy affected Healthy 0.04761905 0.00000000 affected 0.00000000 0.02941176 $stdev.unscaled Healthy affected 1007_s_at 0.2182179 0.1714986 1053_at 0.2182179 0.1714986 117_at 0.2182179 0.1714986 121_at 0.2182179 0.1714986 1255_g_at 0.2182179 0.1714986 54670 more rows ... $pivot [1] 1 2 $genes [1] "1007_s_at" "1053_at" "117_at" "121_at" [5] "1255_g_at" 54670 more rows ... $Amean 1007_s_at 1053_at 117_at 121_at 1255_g_at 9.751804 7.224988 6.079756 7.607733 2.818425 54670 more elements ... $method [1] "ls" $design Healthy affected 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 50 more rows ... > contrast.matrix <-makeContrasts(affected-Healthy,levels = design) > fit2 <- contrasts.fit(fit, contrast.matrix) > fit2 <- eBayes(fit2) > results <- decideTests(fit2, adjust="fdr",lfc=1) > summary(results) affected - Healthy -1 1187 0 52504 1 984 > results TestResults matrix 1007_s_at 1053_at 117_at 121_at 1255_g_at 0 0 0 0 0 54670 more rows ... > On Thu, Jan 31, 2013 at 4:48 PM, James W. MacDonald <jmacdon@uw.edu> wrote: > If you just use the expression values from the original authors, I get > just under 9K probesets for this comparison at an FDR of 0.05 and no fold > change criterion. It drops to just under 900 with a 2-fold difference added > in. > > So yeah, seems like a lot to me as well. > > Best, > > Jim > > > On 1/31/2013 4:02 PM, Steve Lianoglou wrote: > >> ... what Jim said. >> >> But also, this 20k differentially expressed (likely probe sets, not >> genes) is raising a red flag for me, no? Am I alone here? >> >> That's .. what's the word I'm looking for ... "a lot". >> >> -steve >> >> On Thu, Jan 31, 2013 at 3:56 PM, James W. MacDonald<jmacdon@uw.edu> >> wrote: >> >>> Hi Roopa, >>> >>> >>> On 1/31/2013 3:45 PM, Roopa Subbaiaih wrote: >>> >>>> Hi Steve, >>>> >>>> This was the script I used- >>>> getwd() >>>> setwd(dir="/CRSP 406-11/DEMOS/GSE14905-a") >>>> ls() >>>> #-----------------------------**------------------# >>>> library(affy) >>>> eset = justRMA() >>>> f<- factor(c(1,1,1,1,1,1,1,1,1,1,**1,1,1,1,1,1,1,1,1,1,1, >>>> >>>> 2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,**2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,** >>>> 2,2,2,2), >>>> labels=c("Healthy", "unaffected")) >>>> design<- model.matrix(~ 0 + f) >>>> design >>>> colnames(design)<-c("Healthy",**"unaffected") >>>> design >>>> library(limma) >>>> fit<- lmFit(eset, design) >>>> library(hgu133plus2.db) >>>> fit$genes$Symbol<- getSYMBOL(fit$genes$ID,"**hgu133plus2.db") >>>> contrast.matrix<-**makeContrasts(affected-**Healthy,levels = design) >>>> fit2<- contrasts.fit(fit, contrast.matrix) >>>> fit2<- eBayes(fit2) >>>> topTable(fit2,coef=1,p=0.05, adjust="fdr") >>>> results<- decideTests(fit2, adjust="fdr", p=0.05) >>>> summary(results) >>>> write.table(results,file="**myresults.txt") >>>> write.fit(). >>>> >>>> I had identified ~54,000 genes of which ~ 20K were differentially >>>> expressed. >>>> >>>> But when I use these genes for pathway analysis the software asks for >>>> fold >>>> change values but not p value so it is easier to analyze the data. >>>> >>>> What I did was - I used the differentially expressed gene table for >>>> further >>>> analysis. That is I converted logFC values to FC(test/control) assuming >>>> that >>>> >>>> FC= FCmean(test)-FCmean(blank) and LogFC is log2 of FC values. >>>> >>>> Once I got test/control values I converted them to fold changes using >>>> "IF" >>>> function in excel sheet to eliminate genes with fold changes between -2 >>>> to >>>> +2. >>>> >>>> Once I did this the number of significant genes drastically reduced to ~ >>>> 2,000. >>>> >>>> Is this the right method? >>>> >>> >>> No. Note that the range of fold changes after 'unlogging' will be 0-INF, >>> and >>> the down-regulated genes will be in the range 0-1 whereas the upregulated >>> genes will be in the range 1-INF. (e.g. two fold up will be 2, whereas 2 >>> fold down will be 1/2 or 0.5). >>> >>> The easiest way to filter is to keep the logFC and filter on -1 and 1. Or >>> you can use the lfc argument to decideTests(). >>> >>> Best, >>> >>> Jim >>> >>> >>> >>> Please advice, thanks, Roopa >>>> >>>> On Thu, Jan 31, 2013 at 3:23 PM, Steve Lianoglou< >>>> mailinglist.honeypot@gmail.com**> wrote: >>>> >>>> Hi, >>>>> >>>>> On Thu, Jan 31, 2013 at 2:54 PM, Roopa Subbaiaih<rss115@case.edu> >>>>> wrote: >>>>> >>>>>> Hi All, >>>>>> >>>>>> Thanks for the reply, I could pull out the the whole information for >>>>>> differentially expressed genes. The criteria used was adjust="fdr", >>>>>> >>>>> p=0.05. >>>>> >>>>>> I came up with ~ 20,000 genes to be differentially expressed. >>>>>> >>>>> Hmm ... surely 20k cannot be correct? >>>>> >>>>> Since I wanted to analyze these genes for deregulated pathways I had to >>>>>> come up with fold change values for further analysis. >>>>>> >>>>>> I assume that for each gene FC= FCmean(test)-FCmean(blank). LogFC is >>>>>> log2 >>>>>> of FC values. >>>>>> >>>>>> When I convert the FC values (test/blank) to foldchanges using IF >>>>>> >>>>> function >>>>> >>>>>> I get lesser number of genes to be deregulated. The criteria was =>2 >>>>>> foldchanges and =<-2 fold changes. >>>>>> >>>>> I'm missing previous context to this email, so -- not sure what the >>>>> "IF function" is, but if you're using limma, the log2fold changes are >>>>> reported for you in the logFC column that is returned from >>>>> `topTable(...)` >>>>> >>>>> -steve >>>>> >>>>> -- >>>>> Steve Lianoglou >>>>> Graduate Student: Computational Systems Biology >>>>> | Memorial Sloan-Kettering Cancer Center >>>>> | Weill Medical College of Cornell University >>>>> Contact Info: http://cbio.mskcc.org/~lianos/**contact<http: cbi="" o.mskcc.org="" ~lianos="" contact=""> >>>>> >>>>> >>>> -- >>> James W. MacDonald, M.S. >>> Biostatistician >>> University of Washington >>> Environmental and Occupational Health Sciences >>> 4225 Roosevelt Way NE, # 100 >>> Seattle WA 98105-6099 >>> >>> >> >> > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > > -- --------------------------------------- Roopa Shree Subbaiaih Post Doctoral Fellow Department of Dermatology School of Medicine Case Western Reserve University Cleveland, OH-44106 Tel:+1 216 368 0211 [[alternative HTML version deleted]]