I tried to use pdInfoBuilder to make a package with the platform information for a budding yeast tiling array from Affymetrix. I wanted to use this package together with Bioconductor's "oligo" package.
I downloaded these files from Affymetrix:
4008948_yeast_3X_A.CEL
Sc03b_MF_v04.bpmap
Sc03b_MR_v04.bpmap
I made a package with the following code:
library(pdInfoBuilder)
baseDir <- '/Users/EricF/projects/ISW2_FUN30/Yoshida_et_al_Pasero_BrdU_Mol_Cell_2014/s_cerevisiae_libraryfile/BPMAP'
bpmap <- list.files(path = baseDir, pattern = '^\\S+_MR_\\S+bpmap$', full.names = TRUE)
stopifnot(length(bpmap) == 1)
cel <- list.files(path = baseDir, pattern = 'CEL$', full.names = TRUE)
stopifnot(length(cel) == 1)
seed <- new("AffyTilingPDInfoPkgSeed",
bpmapFile = bpmap, celFile = cel,
author = "Eric Foss",
email = "efoss@fredhutch.org",
biocViews = "AnnotationData",
genomebuild = "sacCer1",
organism = "Yeast", species = "Saccharomyces cerevisiae")
makePdInfoPackage(seed, destDir = ".")
install.packages('pd.sc03b.mr.v03', repos = NULL, type = 'source')
(I have done this both with the "_MF_" and with the "_MR_" .bpmap files, and I had identical problems with each.) Everything seemed to go well, although I did get 11 warnings, each of which was the following:
1: In result_fetch(res@ptr, n = n) :
SQL statements must be issued with dbExecute() or dbSendStatement() instead of dbGetQuery() or dbSendQuery().
When I loaded the package with "library(pd.sc03b.mr.v03)", everything went smoothly. However, when I look at the lengths of the chromosomes, something is wrong: 5 out of the 16 chromosomes had maximum coordinates that were far greater than the length of the chromosome (e.g. almost twice the length of the entire genome for a single chromosome; the other 11 chromosomes had maximum coordinates that were reasonable). The code below shows an example of how I determined this.
fls <- list.files(pattern = 'CEL$') # 15 .CEL files from PMID 24856221
affyRaw <- read.celfiles(filenames = fls,
pkgname = 'pd.sc03b.mr.v03')
max(pmPosition(affyRaw)) # 23,944,040 -almost twice the length of the cerevisiae genome
My session info is pasted below:
> sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.4
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] pd.sc03b.mr.v03_0.0.1 DBI_1.1.1 oligo_1.56.0 Biobase_2.52.0 oligoClasses_1.54.0 RSQLite_2.2.7
[7] Biostrings_2.60.1 GenomeInfoDb_1.28.1 XVector_0.32.0 IRanges_2.26.0 S4Vectors_0.30.0 BiocGenerics_0.38.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.6 BiocManager_1.30.16 compiler_4.1.0 MatrixGenerics_1.4.0 bitops_1.0-7
[6] iterators_1.0.13 tools_4.1.0 zlibbioc_1.38.0 bit_4.0.4 preprocessCore_1.54.0
[11] memoise_2.0.0 lattice_0.20-44 ff_4.0.4 pkgconfig_2.0.3 rlang_0.4.11
[16] Matrix_1.3-4 foreach_1.5.1 DelayedArray_0.18.0 rstudioapi_0.13 fastmap_1.1.0
[21] GenomeInfoDbData_1.2.6 affxparser_1.64.0 vctrs_0.3.8 bit64_4.0.5 grid_4.1.0
[26] blob_1.2.1 splines_4.1.0 codetools_0.2-18 matrixStats_0.59.0 GenomicRanges_1.44.0
[31] SummarizedExperiment_1.22.0 RCurl_1.98-1.3 cachem_1.0.5 crayon_1.4.1 affyio_1.62.0
>
I would very much appreciate any suggestions for fixing this problem.
Thank you.
Eric
Now that I think about it, you can just get around the issue yourself. This is definitely a bug and should be fixed, but the days of people using
oligoon the regular, and in particular doing Chip-ChIP are probably behind us now, so maybe it's not that pressing. Anyway, you can do the following.Hi James,
Thanks so much! This worked perfectly. And the Arabidopsis thing explains why it was specifically chromosomes 1-5 whose lengths were so far off, since Arabidopsis has 5 chromosomes. I think that this also suggests that the accessor function called "pmChr" in the oligo package treats, for example, "At:TIGRv5;chr1" as equivalent to "Sc:Oct_2003;chr1", which seems like a bad thing. That's probably part of the bug you referred to.
I have one follow-up question: Your package ended up being named "pd.sc03b.mf.v04" whereas mine was called "pd.sc03b.mr.v03". I started with a .bpmap file with a name identical to that of the one you used: Sc03b_MF_v04.bpmap
This .bpmap file has "MF" in its name. As I understand it, "MF" records strands that hybridize to the Forward strand (so reverse strand sequences), whereas "MR" does the opposite. Also, I understand that it doesn't matter whether one chooses the _MF_ or the _MR_ .bpmap file if the technique you are using is not strand-aware (as mine is not). Also, the .bpmap file has "v04" in its name. It seems pretty reasonable to me that you ended up with a package with "mf" and "v04" in its name; so why did I end up with a package with "mr" and "v03" in its name?
Thanks again for the help. There is no way I could have figured this out on my own.
Best wishes,
Eric
Yes, that's exactly the problem. Whenever one tries to write a package that parses a type of file (in this case bpmap files from Affy tiling arrays), usually a canonical example file is used to infer what to expect in the file, and that then defines how you will parse and interpret the contents. The original tiling arrays were for human and mouse, so I would assume that Benilton used one or both of those as examples, and assumed all other tiling arrays would follow the same convention. But Affy has a history (IMO) of breaking that sort of convention and consequently breaking our tools. This is just one example of that. Maybe they decided they have the real estate on the array for three different species, so they put them all on there and now they can make one array that they can sell as three different tiling arrays. Which is smart!
But anyway.
The name of the package comes from the CEL file you used. To see the code for an S4 method like
makePdInfoPackage, you useselectMethod(orshowMethods).If you do the same using your tiling array I am confident you will get 'pd.sc03b.mf.v03'.
Hi James,
Yes, everything checks out exactly as you expected:
Thanks again. This is really helpful.
Best wishes,
Eric
I have run into one more problem here with chromosome lengths: The information for this Affy chip is supposed to refer to the sacCer1 build of the saccharomyces cerevisiae genome. I created a TilingFeatureSet object called "affyRaw":
The maximum coordinates for almost all of my chromosomes in affyRaw are just a little shorter than the lengths of the chromosomes in the sacCer1 genome assembly, which is reasonable. However, for two chromosomes, the maximum coordinates are longer than the length of the chromosome in the genome assembly - in one case longer by 16 base pairs, and in the other by 193 base pairs:
I don't particularly care about this, except if it is an indication that I am somehow misunderstanding which genome build the Affy chip is supposed to be referring to. "getPlatformDesign" indicates that it is sacCer1:
It's possible that this is going back far enough that the system of yeast genome builds wasn't completely formalized, so that there was some ambiguity about what genome "sacCer1" actually refers to. If that's the case, then everything is fine, and I will just accept my results as they are, but I want to make sure that this discrepancy doesn't reflect a mistake that I am making. For example, if somehow these coordinates referred to the most current build of the yeast genome (sacCer3), then all of the chromosome lengths would be longer than the maximum coordinates that I find for each chromosome, so everything would make sense. (My sessionInfo is pasted below.)
Thanks.
Eric
Note that you are the one who provided the genome build information when you built the package. However, it appears this array was developed in 2003, which is circa sacCer1, so that makes sense.
I wouldn't be too worried that some of the probes are off the end of two chromosomes - it's not uncommon for genomes to be updated between builds, and the chromosomal extent when Affy designed this array may have been different from what is reported currently.
Good point - I had forgotten that I was the one who provided that "sacCer1" information, so my "confirmation" that I was working with the correct genome assembly was meaningless.
I think that you are right and that I shouldn't be concerned that the lengths of 2 out of 16 chromosomes are a bit off. Thanks very much for the help.
Eric