On Wed, Dec 7, 2011 at 7:22 AM, Tatjana <tatjana.milic86@gmail.com> wrote:

> I am working on RNA Seq analysis. I mapped the reads using bwa. As a
> reference,
> I used genome from Ensembl. As I wanted to get a count table, where I will
> have
> a list of genes, and the number of reads which were uniquely mapped to it,
> I
> tried to use countOverlaps in R.
>
> First question I have is regarding difference between 2 functions:
> readAligned
> and readGappedAlignments.
>
> I inputed both: bam files (using function readAligned) and indexed bam
> files
> (using function readGappedAlignments).  When I tried do use countOverlaps,
> I
> encountered some problems. Errors are:
>
> Error in countOverlaps(fcatusDb, tmp[countOverlaps(tmp, fcatusDb) == 1]) :
>  error in evaluating the argument 'subject' in selecting a method for
> function
> 'countOverlaps': Error in tmp[countOverlaps(tmp, fcatusDb) == 1] :
>  error in evaluating the argument 'i' in selecting a method for function
> '[':
> Error in function (classes, fdef, mtable)  :
>  unable to find an inherited method for function "countOverlaps", for
> signature
> "GRanges", "TranscriptDb"
>
>
This means that countOverlaps() does not support passing it a GRanges and
TranscriptDb. In fact, countOverlaps does not support TranscriptDb at all,
because there is no obvious translation from TranscriptDb to a set of
ranges on the genome. If you're trying to overlap with transcripts,
considering only the exonic regions, you will want to call exonsBy() on the
db object, to get a GRangesList. Also, if your reads have intron gaps in
their alignments, you should call grglist() on bwa_output to get the
appropriate GRangesList. Both of those GRangesLists are then passed to
countOverlaps. If you're counting overlaps, you might also want to have a
look at summarizeOverlaps().

Michael


>
> The whole code I used is:
>
> bwa_output=readGappedAlignments("out_trial.prefix.bam", format="BAM")
> fcatusDb=makeTranscriptDbFromBiomart(biomart="ensembl", dataset =
> "fcatus_gene_ensembl")
>
> tmp=as(bwa_output, "Granges")
> #also tried with "GrangesList", the same errors occur
> rc=countOverlaps(fcatusDb, tmp[countOverlaps(tmp, fcatusDb)==1])
>
>

> If someone has experience with this,or know how to solve this, please
> reply to
> this post.
>
> Thank you.
> Tatjana
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

