Limma/Affy - Include expression Intensity in Final 'decideTests' results
2
0
Entering edit mode
Bade ▴ 310
@bade-5877
Last seen 3.4 years ago
Delaware
Hi All, I am using 'affy' and 'limma' to analyze affymetrix data. Simplified code: mydata <- ReadAffy() esetRMA <- rma(mydata) sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') group <-factor(sampletype) design<- model.matrix(~0+group) fit <- lmFit(esetRMA, design) contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) fit2<- contrasts.fit(fit,contrast.matrix) fit2 <- eBayes(fit2) results.default <- decideTests(fit2) I get the results without any problem. But what I need is expression intensity of probes across all samples to be included in the results. I tried: Intensity <- (2^(exprs(esetRMA))) But cannot fit data to the results fData(results.default) <- Intensity Could somebody please tell me how can include expression intensities in the 'decideTests' output? AK
• 1.3k views
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 18 hours ago
United States
Hi Atul, On 4/28/2013 11:55 PM, Atul Kakrana wrote: > Hi All, > > I am using 'affy' and 'limma' to analyze affymetrix data. > > Simplified code: > > mydata <- ReadAffy() > esetRMA <- rma(mydata) > sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') > group <-factor(sampletype) > design<- model.matrix(~0+group) > fit <- lmFit(esetRMA, design) > contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) > fit2<- contrasts.fit(fit,contrast.matrix) > fit2 <- eBayes(fit2) > results.default <- decideTests(fit2) > > I get the results without any problem. But what I need is expression > intensity of probes across all samples to be included in the results. > I tried: > Intensity <- (2^(exprs(esetRMA))) > > But cannot fit data to the results > fData(results.default) <- Intensity > > Could somebody please tell me how can include expression intensities > in the 'decideTests' output? The output from decideTests() is just a matrix of 0 and 1, so you wouldn't want to append anything to that to begin with. You might take a look at limma2annaffy() in affycoretools, which does essentially what I think you want. Best, Jim > > AK > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099
ADD COMMENT
0
Entering edit mode
Hi James, Thanks for the suggestion. By using limma2annaffy I am being able to get expression intensity in final results but the problem is that a separate file is generated for each comparison For Ex. if we have two time points and one control than it will generate two different files. I was thinking of some way to add expression intensity to the results from 'decideTests' that has all the comparisons at one place as that is important for time series comparison. I think one way to do so is just extract the expression intensities from normalized data and add to final result file manually using MS Excel. Besides that I see an error when 'text=TRUE' option is used in limma2annafyy: >limma2annaffy(esetRMA.Reduced, fit2, design, adjust = "none",contrast.matrix, annotation(esetRMA), pfilt = 0.01, fldfilt= 1, text = TRUE, html =FALSE) You are going to output 3 tables, With this many genes in each: 1918 1950 2424 Do you want to accept or change these values? [ a/c ] a Warning message: In chkPkgs(chip) : The mouse4302.db package does not appear to contain annotation data. The error does not appears with just html output. Thanks Atul On 29-Apr-13 9:49 AM, James W. MacDonald wrote: > Hi Atul, > > > On 4/28/2013 11:55 PM, Atul Kakrana wrote: >> Hi All, >> >> I am using 'affy' and 'limma' to analyze affymetrix data. >> >> Simplified code: >> >> mydata <- ReadAffy() >> esetRMA <- rma(mydata) >> sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') >> group <-factor(sampletype) >> design<- model.matrix(~0+group) >> fit <- lmFit(esetRMA, design) >> contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) >> fit2<- contrasts.fit(fit,contrast.matrix) >> fit2 <- eBayes(fit2) >> results.default <- decideTests(fit2) >> >> I get the results without any problem. But what I need is expression >> intensity of probes across all samples to be included in the results. >> I tried: >> Intensity <- (2^(exprs(esetRMA))) >> >> But cannot fit data to the results >> fData(results.default) <- Intensity >> >> Could somebody please tell me how can include expression intensities >> in the 'decideTests' output? > > The output from decideTests() is just a matrix of 0 and 1, so you > wouldn't want to append anything to that to begin with. > > You might take a look at limma2annaffy() in affycoretools, which does > essentially what I think you want. > > Best, > > Jim > > >> >> AK >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor@r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD REPLY
0
Entering edit mode
Hi Atul, You might then try writeFit(), which is probably more to your liking. Something like write.table(writeFit(fit2, "mouse4302", esetRMA.Reduced), "All_data.txt", sep = "\t", quote=FALSE, row.names = FALSE) Also please note that Excel has an unfortunate habit of changing any gene symbol that looks like a date into a date, and this is irreversible (so e.g., the 1Apr gene will be converted to 4/1/2013, because that makes sense...). You want to import the data into Excel rather than opening, and in the third step of the wizard ensure that you change the format of the symbol column to text rather than general. Best, Jim On 4/30/2013 12:32 PM, Atul Kakrana wrote: > Hi James, > > Thanks for the suggestion. By using limma2annaffy I am being able to > get expression intensity in final results but the problem is that a > separate file is generated for each comparison For Ex. if we have two > time points and one control than it will generate two different files. > > I was thinking of some way to add expression intensity to the results > from 'decideTests' that has all the comparisons at one place as that > is important for time series comparison. I think one way to do so is > just extract the expression intensities from normalized data and add > to final result file manually using MS Excel. > > Besides that I see an error when 'text=TRUE' option is used in > limma2annafyy: > > > limma2annaffy(esetRMA.Reduced, fit2, design, adjust = "none",contrast.matrix, annotation(esetRMA), pfilt = 0.01, fldfilt= 1, text = TRUE, html =FALSE) > > You are going to output 3 tables, > With this many genes in each: > 1918 1950 2424 > Do you want to accept or change these values? [ a/c ] > a > Warning message: > In chkPkgs(chip) : > The mouse4302.db package does not appear to contain annotation data. > > The error does not appears with just html output. > > Thanks > > Atul > > > > > On 29-Apr-13 9:49 AM, James W. MacDonald wrote: >> Hi Atul, >> >> >> On 4/28/2013 11:55 PM, Atul Kakrana wrote: >>> Hi All, >>> >>> I am using 'affy' and 'limma' to analyze affymetrix data. >>> >>> Simplified code: >>> >>> mydata <- ReadAffy() >>> esetRMA <- rma(mydata) >>> sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') >>> group <-factor(sampletype) >>> design<- model.matrix(~0+group) >>> fit <- lmFit(esetRMA, design) >>> contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) >>> fit2<- contrasts.fit(fit,contrast.matrix) >>> fit2 <- eBayes(fit2) >>> results.default <- decideTests(fit2) >>> >>> I get the results without any problem. But what I need is expression >>> intensity of probes across all samples to be included in the >>> results. I tried: >>> Intensity <- (2^(exprs(esetRMA))) >>> >>> But cannot fit data to the results >>> fData(results.default) <- Intensity >>> >>> Could somebody please tell me how can include expression intensities >>> in the 'decideTests' output? >> >> The output from decideTests() is just a matrix of 0 and 1, so you >> wouldn't want to append anything to that to begin with. >> >> You might take a look at limma2annaffy() in affycoretools, which does >> essentially what I think you want. >> >> Best, >> >> Jim >> >> >>> >>> AK >>> >>> _______________________________________________ >>> Bioconductor mailing list >>> Bioconductor at r-project.org >>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>> Search the archives: >>> http://news.gmane.org/gmane.science.biology.informatics.conductor >> > -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099
ADD REPLY
0
Entering edit mode
Howdy, On Tue, Apr 30, 2013 at 9:59 AM, James W. MacDonald <jmacdon at="" uw.edu=""> wrote: > Hi Atul, > > You might then try writeFit(), which is probably more to your liking. > Something like > > write.table(writeFit(fit2, "mouse4302", esetRMA.Reduced), "All_data.txt", > sep = "\t", quote=FALSE, row.names = FALSE) > > Also please note that Excel has an unfortunate habit of changing any gene > symbol that looks like a date into a date, and this is irreversible (so > e.g., the 1Apr gene will be converted to 4/1/2013, because that makes > sense...). Tim had previously pointed my attention to this package: http://cran.r-project.org/web/packages/HGNChelper/index.html which has a function called `findExcelGeneSymbols` with the following description: """ This function identifies gene symbols which may have been mogrified by Excel or other spreadsheet programs. If output is assigned to a variable, it returns a vector of the same length where symbols which could be mapped have been mapped. """ I've never used it, so I'm not sure how well it works or not -- I mean, I wouldn't trust it, but gotta give credit to the author for trying :-) -steve -- Steve Lianoglou Computational Biologist Department of Bioinformatics and Computational Biology Genentech
ADD REPLY
0
Entering edit mode
Hi Jim, You suggestion did exactly what I wanted but number of genes in output is different for both: results.default <- decideTests(fit2) write.fit(fit2,results.default,file='RESULTS2.tsv', sep = '\t') AND write.table(writeFit(fit2,"mouse4302", esetRMA.Reduced), "TestData.tsv", sep = '\t', quote = FALSE, row.names = FALSE) The decideTest gives total of ~16000 genes where as 'writeFit' gives ~20,000 genes. I believe this is because of the default p-val cutoff in decideTests : decideTests(object,method="separate",adjust.method="BH",p.value=0.05,l fc=0) As 'writeFit' currently does not support such functions (like decideTests) I think it would be better to add expression intensities manually to the results. Thanks Atul On 04/30/2013 12:59 PM, James W. MacDonald wrote: > Hi Atul, > > You might then try writeFit(), which is probably more to your liking. > Something like > > write.table(writeFit(fit2, "mouse4302", esetRMA.Reduced), > "All_data.txt", sep = "\t", quote=FALSE, row.names = FALSE) > > Also please note that Excel has an unfortunate habit of changing any > gene symbol that looks like a date into a date, and this is > irreversible (so e.g., the 1Apr gene will be converted to 4/1/2013, > because that makes sense...). > > You want to import the data into Excel rather than opening, and in the > third step of the wizard ensure that you change the format of the > symbol column to text rather than general. > > Best, > > Jim > > On 4/30/2013 12:32 PM, Atul Kakrana wrote: >> Hi James, >> >> Thanks for the suggestion. By using limma2annaffy I am being able to >> get expression intensity in final results but the problem is that a >> separate file is generated for each comparison For Ex. if we have two >> time points and one control than it will generate two different files. >> >> I was thinking of some way to add expression intensity to the results >> from 'decideTests' that has all the comparisons at one place as that >> is important for time series comparison. I think one way to do so is >> just extract the expression intensities from normalized data and add >> to final result file manually using MS Excel. >> >> Besides that I see an error when 'text=TRUE' option is used in >> limma2annafyy: >> >> > limma2annaffy(esetRMA.Reduced, fit2, design, adjust = >> "none",contrast.matrix, annotation(esetRMA), pfilt = 0.01, fldfilt= >> 1, text = TRUE, html =FALSE) >> >> You are going to output 3 tables, >> With this many genes in each: >> 1918 1950 2424 >> Do you want to accept or change these values? [ a/c ] >> a >> Warning message: >> In chkPkgs(chip) : >> The mouse4302.db package does not appear to contain annotation data. >> >> The error does not appears with just html output. >> >> Thanks >> >> Atul >> >> >> On 29-Apr-13 9:49 AM, James W. MacDonald wrote: >>> Hi Atul, >>> >>> >>> On 4/28/2013 11:55 PM, Atul Kakrana wrote: >>>> Hi All, >>>> >>>> I am using 'affy' and 'limma' to analyze affymetrix data. >>>> >>>> Simplified code: >>>> >>>> mydata <- ReadAffy() >>>> esetRMA <- rma(mydata) >>>> sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') >>>> group <-factor(sampletype) >>>> design<- model.matrix(~0+group) >>>> fit <- lmFit(esetRMA, design) >>>> contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) >>>> fit2<- contrasts.fit(fit,contrast.matrix) >>>> fit2 <- eBayes(fit2) >>>> results.default <- decideTests(fit2) >>>> >>>> I get the results without any problem. But what I need is >>>> expression intensity of probes across all samples to be included in >>>> the results. I tried: >>>> Intensity <- (2^(exprs(esetRMA))) >>>> >>>> But cannot fit data to the results >>>> fData(results.default) <- Intensity >>>> >>>> Could somebody please tell me how can include expression >>>> intensities in the 'decideTests' output? >>> >>> The output from decideTests() is just a matrix of 0 and 1, so you >>> wouldn't want to append anything to that to begin with. >>> >>> You might take a look at limma2annaffy() in affycoretools, which >>> does essentially what I think you want. >>> >>> Best, >>> >>> Jim >>> >>> >>>> >>>> AK >>>> >>>> _______________________________________________ >>>> Bioconductor mailing list >>>> Bioconductor at r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>> Search the archives: >>>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>> >> >
ADD REPLY
0
Entering edit mode
Hi Atul, The writeFit function simply outputs ALL of the probesets, along with statistics and annotation. If you just want the probesets that are significant in one or more contrast, you could subset using the results.default object. Or add things manually. Your choice. Best, Jim On 4/30/2013 3:05 PM, Atul Kakrana wrote: > Hi Jim, > > You suggestion did exactly what I wanted but number of genes in output > is different for both: > > results.default<- decideTests(fit2) > write.fit(fit2,results.default,file='RESULTS2.tsv', sep = '\t') > > AND > > write.table(writeFit(fit2,"mouse4302", esetRMA.Reduced), "TestData.tsv", > sep = '\t', quote = FALSE, row.names = FALSE) > > The decideTest gives total of ~16000 genes where as 'writeFit' gives > ~20,000 genes. I believe this is because of the default p-val cutoff in > decideTests : > decideTests(object,method="separate",adjust.method="BH",p.value=0.05 ,lfc=0) > > As 'writeFit' currently does not support such functions (like > decideTests) I think it would be better to add expression intensities > manually to the results. > > Thanks > > Atul > > > On 04/30/2013 12:59 PM, James W. MacDonald wrote: >> Hi Atul, >> >> You might then try writeFit(), which is probably more to your liking. >> Something like >> >> write.table(writeFit(fit2, "mouse4302", esetRMA.Reduced), >> "All_data.txt", sep = "\t", quote=FALSE, row.names = FALSE) >> >> Also please note that Excel has an unfortunate habit of changing any >> gene symbol that looks like a date into a date, and this is >> irreversible (so e.g., the 1Apr gene will be converted to 4/1/2013, >> because that makes sense...). >> >> You want to import the data into Excel rather than opening, and in the >> third step of the wizard ensure that you change the format of the >> symbol column to text rather than general. >> >> Best, >> >> Jim >> >> On 4/30/2013 12:32 PM, Atul Kakrana wrote: >>> Hi James, >>> >>> Thanks for the suggestion. By using limma2annaffy I am being able to >>> get expression intensity in final results but the problem is that a >>> separate file is generated for each comparison For Ex. if we have two >>> time points and one control than it will generate two different files. >>> >>> I was thinking of some way to add expression intensity to the results >>> from 'decideTests' that has all the comparisons at one place as that >>> is important for time series comparison. I think one way to do so is >>> just extract the expression intensities from normalized data and add >>> to final result file manually using MS Excel. >>> >>> Besides that I see an error when 'text=TRUE' option is used in >>> limma2annafyy: >>> >>>> limma2annaffy(esetRMA.Reduced, fit2, design, adjust = >>> "none",contrast.matrix, annotation(esetRMA), pfilt = 0.01, fldfilt= >>> 1, text = TRUE, html =FALSE) >>> >>> You are going to output 3 tables, >>> With this many genes in each: >>> 1918 1950 2424 >>> Do you want to accept or change these values? [ a/c ] >>> a >>> Warning message: >>> In chkPkgs(chip) : >>> The mouse4302.db package does not appear to contain annotation data. >>> >>> The error does not appears with just html output. >>> >>> Thanks >>> >>> Atul >>> >>> >>> On 29-Apr-13 9:49 AM, James W. MacDonald wrote: >>>> Hi Atul, >>>> >>>> >>>> On 4/28/2013 11:55 PM, Atul Kakrana wrote: >>>>> Hi All, >>>>> >>>>> I am using 'affy' and 'limma' to analyze affymetrix data. >>>>> >>>>> Simplified code: >>>>> >>>>> mydata<- ReadAffy() >>>>> esetRMA<- rma(mydata) >>>>> sampletype<- c('1','1','1','2','2','2','3','3','3','4','4','4') >>>>> group<-factor(sampletype) >>>>> design<- model.matrix(~0+group) >>>>> fit<- lmFit(esetRMA, design) >>>>> contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) >>>>> fit2<- contrasts.fit(fit,contrast.matrix) >>>>> fit2<- eBayes(fit2) >>>>> results.default<- decideTests(fit2) >>>>> >>>>> I get the results without any problem. But what I need is >>>>> expression intensity of probes across all samples to be included in >>>>> the results. I tried: >>>>> Intensity<- (2^(exprs(esetRMA))) >>>>> >>>>> But cannot fit data to the results >>>>> fData(results.default)<- Intensity >>>>> >>>>> Could somebody please tell me how can include expression >>>>> intensities in the 'decideTests' output? >>>> The output from decideTests() is just a matrix of 0 and 1, so you >>>> wouldn't want to append anything to that to begin with. >>>> >>>> You might take a look at limma2annaffy() in affycoretools, which >>>> does essentially what I think you want. >>>> >>>> Best, >>>> >>>> Jim >>>> >>>> >>>>> AK >>>>> >>>>> _______________________________________________ >>>>> Bioconductor mailing list >>>>> Bioconductor at r-project.org >>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>>> Search the archives: >>>>> http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099
ADD REPLY
0
Entering edit mode
boczniak767 ▴ 720
@maciej-jonczyk-3945
Last seen 20 days ago
Poland
Hi Atul, > > Hi All, > > I am using 'affy' and 'limma' to analyze affymetrix data. > > Simplified code: > > mydata <- ReadAffy() > esetRMA <- rma(mydata) > sampletype <- c('1','1','1','2','2','2','3','3','3','4','4','4') > group <-factor(sampletype) > design<- model.matrix(~0+group) > fit <- lmFit(esetRMA, design) > contrast.matrix<-makeContrasts(B-A,C-A,D-A, levels=design) > fit2<- contrasts.fit(fit,contrast.matrix) > fit2 <- eBayes(fit2) > results.default <- decideTests(fit2) > > I get the results without any problem. But what I need is expression > intensity of probes across all samples to be included in the results. > I > tried: > Intensity <- (2^(exprs(esetRMA))) > > But cannot fit data to the results > fData(results.default) <- Intensity > > Could somebody please tell me how can include expression intensities > in > the 'decideTests' output? You could try make new object with binded columns with fold-changes. If you used BH p-value correction and "separate" method you could write: result_table=cbind(fit2$coefficients,fit2$p.value,p.adjust(fit2$p.valu e,"BH"),separate) So you'll get table with: gene id, fold-change (coefficients), p-value, adjusted p-value, significance code. Where "separate" is a result of: separate=decideTests(fit2,method="separate",adjust.method="BH",p.value =0.05) HTH, Maciej > > AK -- Dr Maciej Jonczyk, Department of Plant Molecular Ecophysiology Faculty of Biology, University of Warsaw 02-096 Warsaw, Miecznikowa 1 Poland -- This email was Anti Virus checked by Astaro Security Gateway. http://www.astaro.com
ADD COMMENT

Login before adding your answer.

Traffic: 778 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6