[BioC] pileup function in ShortRead package

[Patrick Aboyoun]

Nora,
In cases like yours, it is useful to provide the error message along  
with the code. In particular you reference objects aln_Pat and  
pos_Patient that are not created in your code snippet. After I made  
some modifications to your code, I was able to run it without issue:

> suppressMessages(library(ShortRead))
> sp <- SolexaPath(system.file("extdata", package="ShortRead"))
> filt <- compose(chromosomeFilter("chr[1-5].fa"), strandFilter("+"))
> aln_Pat <- readAligned(sp, "s_2_export.txt", filter=filt)
> aln_Pat_subset<-aln_Pat[1:10]
> pos_Patient_subset<-position(aln_Pat_subset)
> str_Patient_subset<-strand(aln_Pat_subset)
> subset_pileuptest<-pileup(pos_Patient_subset, 36,  
> max(pos_Patient_subset,na.rm=TRUE)+37,dir=str_Patient_subset)
> head(subset_pileuptest)
[1] 0 0 0 0 0 0
> sessionInfo()
R version 2.9.2 (2009-08-24)
x86_64-unknown-linux-gnu

locale:
LC_CTYPE=en_US;LC_NUMERIC=C;LC_TIME=en_US;LC_COLLATE=en_US;LC_MONETARY=C;LC_MESSAGES=en_US;LC_PAPER=en_US;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US;LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ShortRead_1.2.1    lattice_0.17-26    BSgenome_1.12.4     
Biostrings_2.12.10
[5] IRanges_1.2.3

loaded via a namespace (and not attached):
[1] Biobase_2.4.1 grid_2.9.2    hwriter_1.1   tools_2.9.2




Quoting Nora Rieber <n.rieber at dkfz-heidelberg.de>:

> Hi there,
>
> I'm using the pileup function in the ShortRead package and it seems I
> call it with the wrong arguments, or that there is some kind of bug in
> the function. All my reads are piled up in forward direction only, even
> if they were actually mapped to the backward strand.
>
> Here's a small (nonfunctional since I assume I cannot append data to
> this e-mail) example of how I call the pileup function:
>
> #read in the MAQ alignment file:
> aln_s_7_MCIP4<-readAligned("/home","sequenceALN.txt", type="MAQMapview")
> aln_Pat_subset<-aln_Pat[1:10]
> pos_Patient_subset<-position(aln_Pat_subset)
> str_Patient_subset<-strand(aln_Pat_subset)
> #use only a small subset (all positions are on the same chromosome 1)
>
> subset_pileuptest<-pileup(pos_Patient_subset, 36,
> max(pos_Patient,na.rm=TRUE)+37,dir=str_Patient_subset)
>
> Is there anything wrong with my dir argument? I checked several places
> all over the genome for the complete alignment and my reads all seem to
> be considered on the + strand, even when they are not. Can anybody help?
>
> This is my session info:
>
>> sessionInfo()
> R version 2.9.2 (2009-08-24)
> x86_64-unknown-linux-gnu
>
> locale:
> LC_CTYPE=de_DE.UTF-8;LC_NUMERIC=C;LC_TIME=de_DE.UTF-8;LC_COLLATE=de_DE.UTF-8;LC_MONETARY=C;LC_MESSAGES=de_DE.UTF-8;LC_PAPER=de_DE.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=de_DE.UTF-8;LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] ShortRead_1.2.1   lattice_0.17-25   BSgenome_1.12.3   Biostrings_2.12.9
> [5] IRanges_1.2.3
>
> loaded via a namespace (and not attached):
> [1] Biobase_2.4.1 grid_2.9.2    hwriter_1.1
>
> Thanks,
> Nora
>
>