Question

beadarray error conversion

0

Entering edit mode

emiliomastriani • 0

@emiliomastriani-11983

Last seen 6.5 years ago

Dear guys,

after a long job I got my limmaRes as result of limmaDE analysis, as you can see from the next:

str(limmaRes)
Formal class 'limmaResults' [package "beadarray"] with 11 slots
  ..@ DesignMatrix     : num [1:16, 1:2] 0 0 0 0 0 0 1 1 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:16] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "cancer" "normal"
  .. ..- attr(*, "assign")= int [1:2] 1 1
  .. ..- attr(*, "contrasts")=List of 1
  .. .. ..$ as.factor(SampleGroup): chr "contr.treatment"
  ..@ ArrayWeights     : Named num [1:16] 1.82 2.9 1.97 2.15 2.01 ...
  .. ..- attr(*, "names")= chr [1:16] "1" "2" "3" "4" ...
  ..@ ContrastMatrix   : num [1:2, 1] 1 -1
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ Levels   : chr [1:2] "cancer" "normal"
  .. .. ..$ Contrasts: chr "cancer-normal"
  ..@ annotation       : chr "Humanv3"

etc etc etc ...

Because most of tools to conduct the next steps require an expression set (GeneAnswers, GOsummaries, etc), I wish to convert the limmaResults to ExpressionSet, and it could be possible (as you can read):

showMethods("coerce")

from="eSet", to="ExpressionSet"
from="eSet", to="MultiSet"
.....

from="limmaResults", to="eSet"
from="limmaResults", to="GRanges"
from="limmaResults", to="Versioned"

.........

But something is not working:

> eset <- as(limmaRes, "eSet")
> dim(eset)
 Features Contrasts 
    34362         1 
> as(eset, "ExpressionSet")
Error in updateOldESet(from, "ExpressionSet") : 
  no slot of name "pData" for this object of class "AnnotatedDataFrame"

        </td>
    </tr>
    <tr>
        <td>
        <p>Someone can help me?</p>

        <p>Thanks</p>

        <pre>

> sessionInfo() R version 3.3.2 (2016-10-31) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 14.04.4 LTS

locale: [1] LCCTYPE=enUS.UTF-8 LCNUMERIC=C LCTIME=enUS.UTF-8 LCCOLLATE=enUS.UTF-8
[5] LCMONETARY=enUS.UTF-8 LCMESSAGES=enUS.UTF-8 LCPAPER=enUS.UTF-8 LCNAME=C
[9] LCADDRESS=C LCTELEPHONE=C LCMEASUREMENT=enUS.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets methods base

other attached packages: [1] beadarray2.24.0 ggplot22.2.0 illuminaHumanv3.db1.26.0 org.Hs.eg.db3.4.0 AnnotationDbi1.36.0
[6] IRanges2.8.1 S4Vectors0.12.1 limma3.30.6 GEOquery2.40.0 Biobase2.34.0
[11] BiocGenerics0.20.0 BiocInstaller1.24.0

loaded via a namespace (and not attached): [1] Rcpp0.12.8 plyr1.8.4 GenomeInfoDb1.10.1 XVector0.14.0 bitops1.0-6 tools3.3.2
[7] zlibbioc1.20.0 convert1.50.0 base642.0 digest0.6.10 RSQLite1.1-1 memoise1.0.0
[13] tibble1.2 gtable0.2.0 DBI0.5-1 httr1.2.1 stringr1.1.0 grid3.3.2
[19] R62.2.0 marray1.52.0 XML3.98-1.5 BeadDataPackR1.26.0 magrittr1.5 reshape21.4.2
[25] scales0.4.1 GenomicRanges1.26.1 assertthat0.1 colorspace1.3-1 labeling0.3 stringi1.1.2
[31] openssl0.9.5 RCurl1.95-4.8 lazyeval0.2.0 munsell0.4.3 illuminaio_0.16.0

        <p>&nbsp;</p>
        </td>
    </tr>
    <tr>
        <td>&nbsp;</td>
    </tr>
</tbody>

beadarray limma ExpressionSet • 914 views

ADD COMMENT • link updated 7.3 years ago by Gordon Smyth 50k • written 7.4 years ago by emiliomastriani • 0

score 0 · Answer 1 · 2016-12-13

I am not familiar with these classes made by the beadarray package, but it seems nonsensical to me to convert the results of a limma linear model analysis back to an eSet, which is supposed to hold expression data.

It also seems to me that your analysis of this data is flawed in that you are treating all the cancer cell lines as equivalent whereas they are in reality very different to each other.

Have you considered just using the limma package directly? You wouldn't have any conversion issues like this. It's very easy to read and analyze this GEO dataset in just a few minutes. For example:

ftpsite <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE16nnn/GSE16570/matrix"
filename <- "GSE16570_series_matrix.txt.gz"
download.file(paste(ftpsite,filename,sep="/"), destfile=filename)
x <- read.delim(filename, comment.char="!", row.names="ID_REF")

CellLine <- factor(c("NOSE","NOSE","NOSE","NOSE","NOSE","NOSE","ES2","KK","OVAS",
              "OVISE","OVMANA","OVSAYO","OVTOKO","RMG1","SMOV2","TOV21G"))
CellLine <- relevel(CellLine, ref="NOSE")
plotMDS(log2(x), label=CellLine)
design <- model.matrix(~CellLine)
fit <- lmFit(log2(x), design)
fit <- eBayes(fit, robust=TRUE, trend=TRUE)

Then you can easily get DE genes in any cell line. For example, to compare ES2 to normal:

> topTable(fit, coef="CellLineES2")
             logFC AveExpr     t  P.Value adj.P.Val    B
ILMN_1653814  7.43    6.32  66.7 2.04e-13  9.97e-09 17.8
ILMN_2138801  4.33    5.76  54.9 1.16e-12  2.83e-08 17.1
ILMN_2316236  6.11    5.75  50.6 2.41e-12  3.92e-08 16.8
ILMN_1695969  4.14    5.83  47.8 4.01e-12  4.89e-08 16.5
ILMN_1778668 -8.24   11.27 -41.6 1.39e-11  1.33e-07 15.8
ILMN_1738646  4.05    6.02  40.9 1.64e-11  1.33e-07 15.7
ILMN_1713125  5.66    6.03  40.2 1.91e-11  1.33e-07 15.7
ILMN_2361768  5.55    5.81  39.6 2.17e-11  1.33e-07 15.6
ILMN_1798700  5.05    6.00  38.9 2.55e-11  1.38e-07 15.5
ILMN_1696317  4.30    6.12  37.8 3.31e-11  1.62e-07 15.3

Or compare KK to normal:

> topTable(fit, coef="CellLineKK")
             logFC AveExpr     t  P.Value adj.P.Val    B
ILMN_2287543  4.63    5.96  43.2 1.00e-11  1.55e-07 16.7
ILMN_1786720  7.03    6.91  42.5 1.16e-11  1.55e-07 16.6
ILMN_2191428 -8.04   13.66 -42.2 1.24e-11  1.55e-07 16.5
ILMN_1660317  4.92    6.36  41.7 1.38e-11  1.55e-07 16.5
ILMN_1686573  7.66    8.89  40.5 1.80e-11  1.55e-07 16.3
ILMN_1813295  3.56    5.79  39.4 2.26e-11  1.55e-07 16.1
ILMN_1743445  4.19    6.11  39.0 2.48e-11  1.55e-07 16.0
ILMN_1666222  4.49    6.39  38.9 2.53e-11  1.55e-07 16.0
ILMN_1779448  7.29    8.59  37.6 3.45e-11  1.87e-07 15.8
ILMN_1753729  4.67    6.46  36.9 4.08e-11  1.94e-07 15.7

And so on.