[BioC] Error mapping RNA-seq data to genome

[Martin Morgan]

Hi Paul --

On 08/17/2010 08:09 AM, Paul Geeleher wrote:
> I'm attempting to read some aligned RNA-seq paired end reads in
> Bioconductor. The reads were aligned with maq (.fastq to .bqf to .map)
> and converted to SAM then BAM format by samtools as I couldn't get the
> paired end reads into biocondcutor in MAQ format.
> 
> I'm basically following this tutorial:
> http://bioconductor.org/help/course-materials/2010/CSAMA10/Lab-8-RNAseqUseCase.pdf
> 
> The following step causes me an error:
> 
>> seqlengths(alnRanges) <- seqlengths(hgTxDb)[names(seqlengths(alnRanges))]
> Error in validObject(.Object) :
>   invalid class "GRanges" object: slot 'ranges' contains values
> outside of sequence bounds
> 
> I think it means that there is something aligned to outside the range
> of my reference genome but I don't understand why this would be as I
> used HG18 for everything.

probably one or more of your reads are 'hanging over' one of the
chromosome bounds. You might

  tapply(alnRanges, seqnames(alnRanges), f
         unction(x) c(min(start(x)), max(end(x))))

to get a sense of whether you're on the right track, and maybe something
like

  refLen <- seqlengths(hgTxDb)[as.character(seqnames(alnRanges)]
  alnRanges[refLen < end(alnRanges)

to get the actual ranges? You could trim or discard these

Martin

> Any suggestions as to how I can possibly figure out what went wrong
> here would be greatly appreciated. I can provide more information if
> necessary. Summaries of my GenomicRanges and GenomicFeatures objects
> are below.
> 
> 
>> alnRanges
> GRanges with 11634719 ranges and 1 elementMetadata value
>            seqnames         ranges strand   | pData.alignData.from...notNA...
>               <Rle>      <IRanges>  <Rle>   |                       <integer>
>        [1]     chr1   [3345, 3381]      +   |                             163
>        [2]     chr1   [3391, 3427]      -   |                             147
>        [3]     chr1   [4223, 4259]      +   |                             163
>        [4]     chr1   [4247, 4283]      -   |                             147
>        [5]     chr1   [4275, 4311]      +   |                             163
>        [6]     chr1   [4304, 4340]      +   |                             163
>        [7]     chr1   [4320, 4356]      +   |                             163
>        [8]     chr1   [4343, 4379]      -   |                             147
>        [9]     chr1   [4343, 4379]      -   |                             147
>        ...      ...            ...    ... ...                             ...
> [11634711]     chrM [16531, 16567]      -   |                             147
> [11634712]     chrM [16532, 16568]      +   |                             161
> [11634713]     chrM [16532, 16568]      -   |                             147
> [11634714]     chrM [16532, 16568]      -   |                             147
> [11634715]     chrM [16532, 16568]      -   |                             147
> [11634716]     chrM [16534, 16570]      -   |                             147
> [11634717]     chrM [16534, 16570]      -   |                             147
> [11634718]     chrM [16534, 16570]      -   |                             147
> [11634719]     chrM [16541, 16577]      -   |                             147
> 
> seqlengths
>   chr1  chr2  chr3  chr4  chr5  chr6 ... chr20 chr21 chr22  chrX  chrY  chrM
>     NA    NA    NA    NA    NA    NA ...    NA    NA    NA    NA    NA    NA
> 
>> hgTxDb <- makeTranscriptDbFromUCSC(genome="hg18", tablename="ensGene")
>> hgTxDb
> TranscriptDb object:
> | Db type: TranscriptDb
> | Data source: UCSC
> | Genome: hg18
> | UCSC Table: ensGene
> | Type of Gene ID: Ensembl gene ID
> | Full dataset: yes
> | transcript_nrow: 63280
> | exon_nrow: 276075
> | cds_nrow: 225373
> | Db created by: GenomicFeatures package from Bioconductor
> | Creation time: 2010-08-17 15:28:19 +0100 (Tue, 17 Aug 2010)
> | GenomicFeatures version at creation time: 1.0.6
> | RSQLite version at creation time: 0.9-1
>>
> 
> 
> 
> 
> 


-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793