[BioC] How to extract expressed genes from Affymetrix data

michelle_low michelle_low at zoho.com
Mon Jul 23 07:56:01 CEST 2012


Hi all,
 
I have done the differential expression analysis below but I wanna know the total number of expressed genes/transcripts in each cell (wild type and the mir223 knockout cell). Can somebody help? Thanks in advance..
 
 
Regards,
Michelle
 
 
 
 R version 2.13.2 (2011-09-30)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[R.app GUI 1.42 (5910) i386-apple-darwin9.8.0]

[Workspace restored from /Users/mlow/.RData]
[History restored from /Users/mlow/.Rapp.history]

> library(affy)
Loading required package: Biobase

Welcome to Bioconductor

  Vignettes contain introductory material. To view, type
  'browseVignettes()'. To cite Bioconductor, see
  'citation("Biobase")' and for packages 'citation("pkgname")'.

>  library(limma)
>  pd=read.AnnotatedDataFrame("phenodata.txt",header=TRUE,sep="",row.names=1)
> a=ReadAffy(filenames=rownames(pData(pd)),phenoData=pd,verbose=TRUE)
1 reading wildrep1.CEL ...instantiating an AffyBatch (intensity a 1004004x6 matrix)...done.
Reading in : wildrep1.CEL
Reading in : wildrep2.CEL
Reading in : wildrep3.CEL
Reading in : miR223rep1.CEL
Reading in : miR223rep2.CEL
Reading in : miR223rep3.CEL
>  x=rma(a)
Background correcting
Normalizing
Calculating Expression
> c=paste(pd$treatment,pd$n,sep="")
> f=factor(c)
> design=model.matrix(~0+f)
>  colnames(design)=levels(f)
> fit=lmFit(x,design)
> library(mouse4302.db)
Loading required package: AnnotationDbi
Loading required package: org.Mm.eg.db
Loading required package: DBI


> fit$genes$Symbol <- getSYMBOL(fit$genes$ID,"mouse4302.db")
Error: could not find function "getSYMBOL"
> library(annotate)
> fit$genes$Symbol <- getSYMBOL(fit$genes$ID,"mouse4302.db")
> contrast.matrix=makeContrasts(E="present-absent")
Error in .Internal(inherits(x, what, which)) : 'x' is missing
> contrast.matrix=makeContrasts(E="present-absent",levels=design)
> fit2=contrasts.fit(fit,contrast.matrix)
> fit2=eBayes(fit2)
> results1 <-topTable (fit2, coef=1, p.value=0.3,number=nrow(fit2))
> dim(results1)
[1] 889   8
> results1 <-topTable (fit2, coef=1, p.value=0.5,number=nrow(fit2))
> dim(results1)
[1] 2997    8
> results1 <-topTable (fit2, coef=2, p.value=0.5,number=nrow(fit2))
Error in toptable(fit = fit[c("coefficients", "stdev.unscaled")], coef = coef,  : 
  subscript out of bounds
> write.table(results1, file="WT-KO.txt")
> results1 <-topTable (fit2, coef=2, p.value=0.4,number=nrow(fit2))
Error in toptable(fit = fit[c("coefficients", "stdev.unscaled")], coef = coef,  : 
  subscript out of bounds
> results1 <-topTable (fit2, coef=1, p.value=0.4,number=nrow(fit2))
> dim(results1)
[1] 1668    8
> write.table(results1, file="WT-KO.txt")



More information about the Bioconductor mailing list