Open cb4github opened 7 years ago
Hello Carl, Glad to hear you are finding Goldmine to be useful.
To answer your questions:
Q1. Given a query range that overlaps a promoter of a subject gene for which there is no entry in the resulting goldmine gene table (such as row 1 in the example dmrQuery below), then is it true that there is no such entry in the resulting goldmine gene table since the given query range does not overlap with the subject gene body?
Yes this is correct that there will not be a gene table entry for a region from the query regions if it does not overlap with any gene. All query regions do appear in the "context" table, but will only appear in the genes or features tables if the query region does indeed overlap with one of those subjects.
Q2. Given a query range that results in a value of promoter_per > 100 (such as row 2 in the example dmrQuery below), what is the meaning of context$promoter_per in this case?
This is a good question, as promoter_per
by definition should have a max of 100. I ran your code example and looked at the data. I the correct value should be "100". I noticed that the row 2 region actually overlaps promoters of 2 genes (ZCCHC17 at 100% and SNRNP40 at 3.62%), and I suspect that somehow the values have been summed. This is something I'll need to look into further.
Q3. What is the meaning of gene$Promoter?
This is the percentage of the query region (given by chr/start/end of that row) that is overlapped by the promoter of the isoform given in "isoform.id" of that row.
I should point out that the genes
table only shows a row for a query region if that query region directly overlaps with the gene body. Thus, overlaps restricted to promoters and gene 3' regions do not appear as rows in genes
. This is why there is now a row for chr1:31296092-31296616 and ZCCHC17.
I believe I have found the source of this bug, and it was due to strand not being ignored in the calcPercentOverlap
function. When a query region overlapped with two promoters that themselves also overlapped and those two promoters were on different strands the bug would be induced.
Note that Goldmine was written to ignore strand information for simplicity of implementation, and strand should have been ignored here as well. (Stranded analyses are however supported by pre-splitting the query, genes, and features and running goldmine()
once for each strand.) In this instance, stranded information was being sent from the upstream getGeneModels()
function. The fix (commit at https://github.com/jeffbhasin/goldmine/commit/da4370b2ef59dbd853681241f34d786fa50d66c0) was to ignore this information, and then compute the percent overlaps. Overlapping promoters were already dealt with correctly by the code on the presumption that the function is fed an unstranded subject.gr
argument.
I created a branch called issue-3
. Please try it out and if this fixes the issue and induces no other issues I can pull it to master
.
library(devtools)
install_github("jeffbhasin/goldmine",ref="issue-3")
Thanks very much for finding this and please let me know if you have any additional questions!
Hi Jeff,
Thanks for your reply, and I now have a couple of new questions.
Q4. Re SNRNP40, I'm not sure how to get the result 3.62% for overlap of the promoter. I keep getting ~0.9%.
library(intervals)
dmr2 <- Intervals(
matrix(
dmrToRefseqOverlap$genes[,c(query.start, query.end)],
ncol = 2, byrow = FALSE
),
closed = c( TRUE, TRUE ),
type = "Z"
)
promoter2 <- Intervals(
matrix(
dmrToRefseqOverlap$genes[, isoform.end] + c(-200, 2000),
ncol = 2, byrow = FALSE
),
closed = c( TRUE, TRUE ),
type = "Z"
)
> print(dmr2)
Object of class Intervals
1 interval over Z:
[31296092, 31296616]
> print(promoter2)
Object of class Intervals
1 interval over Z:
[31296597, 31298797]
> print(dmrIntersectPromoter2 <- interval_intersection(dmr2, promoter2))
Object of class Intervals
1 interval over Z:
[31296597, 31296616]
> print(sizeDmrIntersectPromoter2 <- size(dmrIntersectPromoter2))
[1] 20
> print(sizePromoter2 <- size(promoter2))
[1] 2201
> print(100*sizeDmrIntersectPromoter2/sizePromoter2)
[1] 0.9086779
>
Q5. In my original example code above, I meant to use GENOME <- "hg19"
, so that when I rerun my example code with that modification, I get a new error as in the below...
...
> dmrToRefseqOverlap <-
+ goldmine(
+ query = dmrQuery,
+ genes = refseqGenes,
+ promoter = c(
+ NUMBER_BP_UPSTREAM_OF_PROMOTER,
+ .... [TRUNCATED]
Computing gene models
Generating context annotation - genes
Generating genes report
Error in `$<-.data.frame`(`*tmp*`, "type", value = "E") :
replacement has 1 row, data has 0
> traceback()
9: stop(sprintf(ngettext(N, "replacement has %d row, data has %d",
"replacement has %d rows, data has %d"), N, nrows), domain = NA)
8: `$<-.data.frame`(`*tmp*`, "type", value = "E")
7: `$<-`(`*tmp*`, "type", value = "E")
6: reportGenes2(query.gr, genes.gr)
5: goldmine(query = dmrQuery, genes = refseqGenes, promoter = c(NUMBER_BP_UPSTREAM_OF_PROMOTER,
NUMBER_BP_DOWNSTREAM_OF_PROMOTER), genome = GENOME, cachedir = CACHE_DIR) at debugGoldmine.R#22
4: eval(ei, envir)
3: eval(ei, envir)
2: withVisible(eval(ei, envir))
1: source("~/research/epigenetics/20161123/lebron/debugGoldmine.R",
echo = TRUE)
Please advise, thanks.
Best, Carl B.
Hello Carl,
I investigated your question Q5 first, and found that this was due to the set of query ranges being a case where there were no intron/exon overlaps. Some code was running that only applies in case there are intron/exon overlaps, so I corrected the error by conditioning the processing over the intron/exon overlap table on there being actual overlap. Before I had really only run Goldmine on datasets with at least one exon/intron overlap where this case didn't arise.
You can see the commits to fix this issue here: https://github.com/jeffbhasin/goldmine/commit/138a66ffd202b66628fd02cf432303143126e2c8 and https://github.com/jeffbhasin/goldmine/commit/c2bae07aefb3fb5fcfb8a3d50847e50e596c088e and you can test it by pulling the latest commit from the "issue-3" branch as before.
Please let me know if any further issues. I will investigate your Q4 next.
Jeff
Hello Carl, I ran your code for Q4 and also did my own arithmetic. Goldmine computes length 19 bp for the overlap because it includes the transcription start bp as part of the upstream. Asides from that, I get the same values as you do. I found that the 3.62 is the percent of the DMR that is overlapped by promoter. The percentages are the fraction of overlap / width of query. I believe this is consistent with my response to Q3. At design time, I thought it was more logical to make these with respect to the DMR, so you see what % of the DMR is in each gene model portion, rather than the other way around.
The below code appended to the end of your own gives the 3.62 value, accounting for the off-by-in in our promoter conventions.
100*(sizeDmrIntersectPromoter2-1)/size(dmr2)
Jeff
Updated: I edited my comment because I think this is consistent with my answer to Q3
Hello Jeff,
I tried your updated version...
source("http://bioconductor.org/biocLite.R")
biocLite(c("GenomicRanges", "IRanges", "devtools"))
# Then, install Goldmine from GitHub. Be sure to accept installation of any additional pre-requisite packages from CRAN.
library(devtools)
install_github("jeffbhasin/goldmine",ref = "issue-3")
... and it seems that another, perhaps old, problem as (re)arisen.
With the following code (tl;dr - please excuse my admittedly lengthy attempt at a simple, reproducible example)...
CACHE_DIR <- "gbcache"
GENOME <- "hg19"
REFSEQ <- "refseq"
refseqGenesDt <- getGenes(REFSEQ, genome = GENOME, cachedir = CACHE_DIR)
dmr1 <- GRanges("chr1", IRanges(start = 754738, end = 755255))
dmr2 <- GRanges("chr1", IRanges(start = 783146, end = 783802))
testQuery <- dmr2
testSubject <- refseqGenesDt[gene.id == "LINC01128"][2:3]
testQuery
testSubject
testOverlap <-
goldmine(
query = testQuery,
genes = testSubject,
promoter = c(
NUMBER_BP_UPSTREAM_FROM_TSS_IN_PROMOTER,
NUMBER_BP_DOWNSTREAM_FROM_TSS_IN_PROMOTER
),
genome = GENOME,
cachedir = CACHE_DIR,
)
testOverlap
testQuery2 <- c(dmr1, dmr2)
testSubject2 <- refseqGenesDt[gene.id == "FAM87B"]
testSubject2 <- rbind(testSubject2, testSubject)
testQuery2
testSubject2
testOverlap2 <-
goldmine(
query = testQuery2,
genes = testSubject2,
promoter = c(
NUMBER_BP_UPSTREAM_FROM_TSS_IN_PROMOTER,
NUMBER_BP_DOWNSTREAM_FROM_TSS_IN_PROMOTER
),
genome = GENOME,
cachedir = CACHE_DIR,
)
testOverlap2
...I got the following output...
> CACHE_DIR <- "gbcache"
> GENOME <- "hg19"
> REFSEQ <- "refseq"
> refseqGenesDt <- getGenes(REFSEQ, genome = GENOME, cachedir = CACHE_DIR)
> dmr1 <- GRanges("chr1", IRanges(start = 754738, end = 755255))
> dmr2 <- GRanges("chr1", IRanges(start = 783146, end = 783802))
> testQuery <- dmr2
> testSubject <- refseqGenesDt[gene.id == "LINC01128"][2:3]
> testQuery
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 [783146, 783802] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> testSubject
chr start end strand name gene.id isoform.id cdsStart cdsEnd exonCount exonStarts
1: chr1 762971 794826 + LINC01128 LINC01128 NR_047519 794826 794826 6 762970,764382,783033,787306,788050,788770,
2: chr1 762971 794826 + LINC01128 LINC01128 NR_047521 794826 794826 5 762970,764382,787306,788050,788770,
exonEnds
1: 763155,764484,783186,787490,788146,794826,
2: 763155,764484,787490,788146,794826,
> testOverlap <-
+ goldmine(
+ query = testQuery,
+ genes = testSubject,
+ promoter = c(
+ NUMBER_BP_UPSTREAM_FROM_TSS_IN_PROMOTER .... [TRUNCATED]
Computing gene models
Generating context annotation - genes
Generating genes report
Generating exon/intron overlap diagrams
> testOverlap
$context
chr start end width strand qrow promoter_per end3_per exon_per intron_per intergenic_per utr5_per utr3_per call call_genes overlapped_genes
1: chr1 783146 783802 657 * 1 0 0 6.24 100 0 0 0 exon LINC01128 LINC01128
nearest_genes distance_to_nearest_gene url
1: LINC01128 0 http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
$genes
qrow srow query.chr query.start query.end gene.symbol gene.id isoform.id isoform.chr isoform.start isoform.end isoform.strand overlap.bp
1: 1 1 chr1 783146 783802 LINC01128 LINC01128 NR_047519 chr1 762971 794826 + 657
2: 1 2 chr1 783146 783802 LINC01128 LINC01128 NR_047521 chr1 762971 794826 + 657
query.overlap.per isoform.overlap.per noncoding Promoter
1: 100 2.06 TRUE 0
2: 100 2.06 TRUE 0
ExonIntron 3' End
1: E1 (0), I1 (0), E2 (0), I2 (0), E3 (6.24), I3 (93.91), E4 (0), I4 (0), E5 (0), I5 (0), E6 (0) 0
2: E1 (0), I1 (0), E2 (0), I2 (100), E3 (0), I3 (0), E4 (0), I4 (0), E5 (0) 0
url
1: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
2: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
$features
list()
> testQuery2 <- c(dmr1, dmr2)
> testSubject2 <- refseqGenesDt[gene.id == "FAM87B"]
> testSubject2 <- rbind(testSubject2, testSubject)
> testQuery2
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 [754738, 755255] *
[2] chr1 [783146, 783802] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> testSubject2
chr start end strand name gene.id isoform.id cdsStart cdsEnd exonCount exonStarts
1: chr1 752751 755214 + FAM87B FAM87B NR_103536 755214 755214 2 752750,754102,
2: chr1 762971 794826 + LINC01128 LINC01128 NR_047519 794826 794826 6 762970,764382,783033,787306,788050,788770,
3: chr1 762971 794826 + LINC01128 LINC01128 NR_047521 794826 794826 5 762970,764382,787306,788050,788770,
exonEnds
1: 753582,755214,
2: 763155,764484,783186,787490,788146,794826,
3: 763155,764484,787490,788146,794826,
> testOverlap2 <-
+ goldmine(
+ query = testQuery2,
+ genes = testSubject2,
+ promoter = c(
+ NUMBER_BP_UPSTREAM_FROM_TSS_IN_PROMO .... [TRUNCATED]
Computing gene models
Generating context annotation - genes
Generating genes report
Generating exon/intron overlap diagrams
> testOverlap2
$context
chr start end width strand qrow promoter_per end3_per exon_per intron_per intergenic_per utr5_per utr3_per call call_genes
1: chr1 754738 755255 518 * 1 0 100 92.08 0 7.92 0 0 3' end FAM87B
2: chr1 783146 783802 657 * 2 0 0 6.24 100 0.00 0 0 exon LINC01128
overlapped_genes nearest_genes distance_to_nearest_gene url
1: FAM87B FAM87B 0 http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A754738-755255
2: LINC01128 LINC01128 0 http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
$genes
qrow srow query.chr query.start query.end gene.symbol gene.id isoform.id isoform.chr isoform.start isoform.end isoform.strand overlap.bp
1: 1 1 chr1 754738 755255 FAM87B FAM87B NR_103536 chr1 752751 755214 + 477
2: 2 2 chr1 783146 783802 LINC01128 LINC01128 NR_047519 chr1 762971 794826 + 657
3: 2 3 chr1 783146 783802 LINC01128 LINC01128 NR_047521 chr1 762971 794826 + 477
query.overlap.per isoform.overlap.per noncoding Promoter
1: 92.08 19.36 TRUE 0
2: 100.00 2.06 TRUE 0
3: 92.08 2.06 TRUE 0
ExonIntron 3' End
1: E1 (0), I1 (0), E2 (92.08) 100
2: E1 (0), I1 (0), E2 (0), I2 (0), E3 (6.24), I3 (93.91), E4 (0), I4 (0), E5 (0), I5 (0), E6 (0) 0
3: E1 (0), I1 (0), E2 (0), I2 (100), E3 (0), I3 (0), E4 (0), I4 (0), E5 (0) 0
url
1: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A754738-755255
2: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
3: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr1%3A783146-783802
$features
list()
> testOverlap$genes[, c(1:14)]
qrow srow query.chr query.start query.end gene.symbol gene.id isoform.id isoform.chr isoform.start isoform.end isoform.strand overlap.bp
1: 1 1 chr1 783146 783802 LINC01128 LINC01128 NR_047519 chr1 762971 794826 + 657
2: 1 2 chr1 783146 783802 LINC01128 LINC01128 NR_047521 chr1 762971 794826 + 657
query.overlap.per
1: 100
2: 100
> testOverlap2$genes[, c(1:14)]
qrow srow query.chr query.start query.end gene.symbol gene.id isoform.id isoform.chr isoform.start isoform.end isoform.strand overlap.bp
1: 1 1 chr1 754738 755255 FAM87B FAM87B NR_103536 chr1 752751 755214 + 477
2: 2 2 chr1 783146 783802 LINC01128 LINC01128 NR_047519 chr1 762971 794826 + 657
3: 2 3 chr1 783146 783802 LINC01128 LINC01128 NR_047521 chr1 762971 794826 + 477
query.overlap.per
1: 92.08
2: 100.00
3: 92.08
Warning messages:
1: In data.table(qrow = qrow, srow = srow, query.chr = query.chr, query.start = query.start, :
Item 13 is of size 2 but maximum size is 3 (recycled leaving remainder of 1 items)
2: In data.table(qrow = qrow, srow = srow, query.chr = query.chr, query.start = query.start, :
Item 14 is of size 2 but maximum size is 3 (recycled leaving remainder of 1 items)
>
...wherein it appears that the second of the two results from goldmine
appears to have some corrupted values for testOverlap2$genes[, 13:14]
as compared to that of testOverlap$genes[, 13:14]
.
In particular in forming testOverlap2
, it appears that the internal results for the two columns 13 and 14 (both of length 2) have been recycled in order to fit the result columns (both of length 3). Looking at the code further in the updated version of the method goldmine:::calcPercentOverlap
, it appears that the erstwhile code #subject.gr <- reduce(subject.gr)
has been reinstated - perhaps intended to fix a prior error.
Hope this helps.
Best, Carl
Dear Team,
First, thanks for all your efforts - in particular for encapsulating in goldmine and abbreviating what has been many a previous effort in mapping DM results to a RefSeq genome.
Please excuse if the following are not so much issues but questions about the code, which presumably then is working as designed.
Q1. Given a query range that overlaps a promoter of a subject gene for which there is no entry in the resulting goldmine
gene
table (such as row 1 in the exampledmrQuery
below), then is it true that there is no such entry in the resulting goldminegene
table since the given query range does not overlap with the subject gene body?Q2. Given a query range that results in a value of
promoter_per
> 100 (such as row 2 in the exampledmrQuery
below), what is the meaning ofcontext$promoter_per
in this case?Q3. What is the meaning of
gene$Promoter
?Best, Carl B.
Here's example code.
Here's the output.