[BioC] summarizeOverlaps - circular MT

[Valerie Obenchain]

Hi Stefanie,

On 08/01/12 05:41, Stefanie Tauber wrote:
> Dear list,
>
> when I try to count human reads vs a human transcript database, I get the
> following error:
>
>
> library(GenomicFeatures)
> library(biomaRt)
>
> Here I construct my human transcript database:
>
> ensembl = useDataset("hsapiens_gene_ensembl",mart=useMart("ensembl"))
>
> res = getBM(attributes = "ensembl_transcript_id", filters =
> ("biotype","status"),values=list("protein_coding","known"),mart = ensembl)
>
> humanDb = makeTranscriptDbFromBiomart(biomart = "ensembl", dataset =
> "hsapiens_gene_ensembl",transcript_ids = as.character(res[,1]))
>
> tx = transcriptsBy(humanDb)
>
>
> The reads are an GappedAlignments object:
> reads
>
> GappedAlignments with 11791172 alignments and 0 elementMetadata cols:
>                      seqnames strand       cigar    qwidth     start       end
>                         <Rle>   <Rle>  <character>  <integer>  <integer>  <integer>
>          SRR015293.3        3      *         32M        32 186338939 186338970
>          SRR015293.5       16      *         32M        32  72094409  72094440
>          SRR015293.7        1      *         32M        32 159683461 159683492
>          SRR015293.9        4      *         32M        32 155529684 155529715
>         SRR015293.10        4      *         32M        32 110632783 110632814
>
>
> This is working fine:
> counts = assays(summarizeOverlaps(tx, reads, mode = "Union"))$counts
>
> Method "IntersectionStrict" does not work:
>   counts1 = assays(summarizeOverlaps(tx, liver, mode = "IntersectionStrict"))$counts
> Error in assays(summarizeOverlaps(tx, liver, mode = "IntersectionStrict")) :
>    error in evaluating the argument 'x' in selecting a method for function
> 'assays': Error in queryHits(findOverlaps(query, subject, maxgap = maxgap,
> minoverlap = minoverlap,  :
>    error in evaluating the argument 'x' in selecting a method for function
> 'queryHits': Error in .findOverlaps.circle(circle.length, seqselect(queryRanges,
> qIdxs),  :
>    overlap type "within" is not yet supported for circular sequence MT
>
>
> What can I do? Do I have to omit all genes on the MT?
Thanks for reporting this. The problem was that findOverlaps for GRanges 
does not currently support type="within" for circular chromosomes. The 
mode 'IntersectionStrict' uses findOverlaps in this way whereas 'Union' 
and 'IntersectionNotEmpty' do not. I have added a check for circular 
chromosomes to the 'IntersectionStrict' mode. This will issue a warning 
when circular chromosomes are encountered and they will be omitted from 
the counting.

This has been fixed in GenomicRanges release 1.8.11 and devel 1.9.41. 
These versions will be available through biocLite() tomorrow morning (~9 
am PST) or through svn immediately.

Valerie
>
> Best,
> Stefanie
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor