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" 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

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]]