jorainer
commented
6 years ago Dear @apastore , thanks for reporting. In the current Bioconductor development version I've moved (and rewritten) the mapping code from Pbase to ensembldb and this also fixed the error you describe.
ensembldb has now 6 functions, proteinToTranscript, proteinToGenome, transcriptToProtein, transcriptToGenome, genomeToTranscript and genomeToProtein to map between protein, transcript and genome coordinates. In order to use these functions you would need to install a current developmental version of R (R-3.5 devel) and install Bioc3.7-devel (using biocLite as you normally would do).
You can do the mapping as described below:
library(ensembldb)
library(AnnotationHub)
edb <- query(AnnotationHub(), c("Ensembl 90 EnsDb", "Homo sapiens"))[[1]]
hgnc_symbol = "MLLT3"
## Get proteins along with protein domains.
prts <- proteins(edb, filter = ~ symbol == hgnc_symbol,
columns = c("protein_id", "tx_id", listColumns(edb, "protein_domain")))
## You've got a DataFrame, each row representing one protein domain.
prts
DataFrame with 6 rows and 8 columns
protein_id tx_id protein_domain_id protein_domain_source
<character> <character> <character> <character>
1 ENSP00000369678 ENST00000380321 NA NA
2 ENSP00000369695 ENST00000380338 PS51037 pfscan
3 ENSP00000369695 ENST00000380338 PF03366 pfam
4 ENSP00000485996 ENST00000630269 PS51037 pfscan
5 ENSP00000485996 ENST00000630269 PF03366 pfam
6 ENSP00000486193 ENST00000629733 NA NA
interpro_accession prot_dom_start prot_dom_end symbol
<character> <integer> <integer> <character>
1 NA NA NA MLLT3
2 IPR005033 8 112 MLLT3
3 IPR005033 29 108 MLLT3
4 IPR005033 5 109 MLLT3
5 IPR005033 26 105 MLLT3
6 NA NA NA MLLT3
## Remove proteins without protein domains
prts <- prts[!is.na(prts$protein_domain_id), ]
## Define an IRanges representing the protein domains (their start and end coords
## within the protein sequence) and use the protein ID as their names
dmns <- IRanges(start = prts$prot_dom_start, end = prts$prot_dom_end,
names = prts$protein_id)
## Now map the domains to the genome
dmns_gnm <- proteinToGenome(dmns, edb)
## The result is a list of GRanges, each GRanges representing the genomic coordinates of
## the protein domain
length(dmns_gnm)
[1] 4
## The genomic coordinates for the first protein domain
dmns_gnm[[1]]
GRanges object with 3 ranges and 7 metadata columns:
seqnames ranges strand | protein_id tx_id
<Rle> <IRanges> <Rle> | <character> <character>
[1] 9 [20620654, 20620825] - | ENSP00000369695 ENST00000380338
[2] 9 [20456704, 20456786] - | ENSP00000369695 ENST00000380338
[3] 9 [20448207, 20448266] - | ENSP00000369695 ENST00000380338
exon_id exon_rank cds_ok protein_start protein_end
<character> <integer> <logical> <integer> <integer>
[1] ENSE00003665608 2 TRUE 8 112
[2] ENSE00003628570 3 TRUE 8 112
[3] ENSE00003670233 4 TRUE 8 112
-------
seqinfo: 1 sequence from GRCh38 genomeHave also a look at the new coordinate_mapping vignette of ensembldb for more details.
My sessionInfo:
> sessionInfo()
R Under development (unstable) (2017-12-14 r73916)
Platform: x86_64-apple-darwin17.3.0/x86_64 (64-bit)
Running under: macOS High Sierra 10.13.3
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.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
[8] methods base
other attached packages:
[1] AnnotationHub_2.11.2 ensembldb_2.3.7 AnnotationFilter_1.3.0
[4] GenomicFeatures_1.31.1 AnnotationDbi_1.41.4 Biobase_2.39.1
[7] GenomicRanges_1.31.3 GenomeInfoDb_1.15.1 IRanges_2.13.4
[10] S4Vectors_0.17.16 BiocGenerics_0.25.1
loaded via a namespace (and not attached):
[1] SummarizedExperiment_1.9.5 progress_1.1.2
[3] lattice_0.20-35 htmltools_0.3.6
[5] rtracklayer_1.39.4 yaml_2.1.16
[7] interactiveDisplayBase_1.17.0 blob_1.1.0
[9] XML_3.98-1.9 rlang_0.1.6
[11] pillar_1.0.1 DBI_0.7
[13] BiocParallel_1.13.1 bit64_0.9-7
[15] matrixStats_0.52.2 GenomeInfoDbData_1.1.0
[17] stringr_1.2.0 zlibbioc_1.25.0
[19] ProtGenerics_1.11.0 Biostrings_2.47.1
[21] memoise_1.1.0 biomaRt_2.35.5
[23] BiocInstaller_1.29.3 httpuv_1.3.5
[25] curl_3.1 Rcpp_0.12.14
[27] xtable_1.8-2 DelayedArray_0.5.15
[29] XVector_0.19.1 mime_0.5
[31] bit_1.1-12 Rsamtools_1.31.1
[33] RMySQL_0.10.13 digest_0.6.13
[35] stringi_1.1.6 shiny_1.0.5
[37] grid_3.5.0 tools_3.5.0
[39] bitops_1.0-6 magrittr_1.5
[41] RCurl_1.95-4.10 lazyeval_0.2.1
[43] tibble_1.4.1 RSQLite_2.0
[45] pkgconfig_2.0.1 Matrix_1.2-12
[47] prettyunits_1.0.2 assertthat_0.2.0
[49] httr_1.3.1 R6_2.2.2
[51] GenomicAlignments_1.15.2 compiler_3.5.0 
apastore
Hi, thanks a lot for this package is extremely useful and save me a lot of code.
I am experiencing some issue when I try to map protein domain to GRCh38. many protein domain get mapped correctly, but for some protein it return a empty GrangeList even if multiple protein domain are presented. One example is human MLLT3.
Do you have any idea?
Thanks a lot!