
Hi list,

sorry for a simple question, but I am a newbie a bit lost reading all the information on how to handle NGS data using R tools. I have a set of BAM files from Junior sequencer, one BAM per amplicon per sample (Roche's software does not output one single BAM file per sample including all amplicons). I also have a bed file with the features sequenced. At the moment I am only testing with two amplicons for one sample. My final aim is to get the coverage per amplicon per sample

I am using R version 2.14.2 (2012-02-29)

I read in the bed file and tried to create a GRanges object:

> bed.frame=read.table("bed.frame",sep="\t",as.is=TRUE, header=TRUE) 
> bed.frame
     chromosome start end
1 APOBex26_M13      1 478
2 APOBex29_M13      1 448
> gr<-GRanges(seqnames=bed.frame[,1],ranges=IRanges(bed.frame[,2],end=bed.frame[,3]))
> gr                                                                                                                                                                                
GRanges with 2 ranges and 0 elementMetadata values:                                                                                                                                 
           seqnames    ranges strand                                                                                                                                                
              <Rle> <IRanges>  <Rle>                                                                                                                                                
  [1] APOBex26_M13   [1, 478]      *                                                                                                                                                
  [2] APOBex29_M13   [1, 448]      *                                                                                                                                                
  ---                                                                                                                                                                               
  seqlengths:                                                                                                                                                                       
   APOBex26_M13 APOBex29_M13                                                                                                                                                       
              NA            NA                             


I also read in the list of BAM files vand convert to BamFileList:

> fls = list.files(pattern="*bam",full=TRUE)
> fls
 [1] "./Sample_Mult_1_MID-151_vs_APOBex26_M13.bam"     
 [2] "./Sample_Mult_1_MID-151_vs_APOBex29_M13.bam"
bfs <- BamFileList(fls)

I then try to summarizeOverlaps but gives me the following error.

>olap <-summarizeOverlaps(gr, bfs)       
Error in .Call2("cigar_to_width", cigar, PACKAGE = "GenomicRanges") : 
  in 'cigar' element 1: CIGAR operation 'P' (at char 7) is not supported yet, sorry!                                                                                        

What is the meaning of this error and how can I overcome it?
Is the gr object not created properly? Is Metadata necessary? Why is not seqlengths filled automatically using the IRanges?

Thanks for your help,

Dave
 		 	   		  
	[[alternative HTML version deleted]]

