[BioC] pileupAsGRanges hangs

Martin Morgan mtmorgan at fhcrc.org
Mon Dec 10 17:27:46 CET 2012


On 12/10/2012 05:03 PM, Mark Dunning wrote:
> Hi,
>
> I'm trying to generate some pileups from targetted-resequence data,
> but finding that the pileupAsGRanges command just hangs on some bam
> file and region combinations. The following works fine;

There was a bug like this fixed in Rsamtools 1.9.17, so a solution might be to 
upgrade to the current release

   source("http://bioconductor.org/biocLite.R")
   biocLite("BiocUpgrade")

You might be able to narrow this to the Rsamtools::applyPileups code vs. code 
from biovizBase (where I'm guessing pileupAsGRanges is coming from?) by invoking

   pileupFiles <- PileupFiles(bams)
   pileupParams <- PileupParam(which = gr, ...)
   applyPileups(pileupFiles, function(...) NULL, param=pileupParams)

If it is applyPileups, and it is not fixed in the current release, then it would 
be great if you could send me a small bam file (e.g., using filterBam to select 
the problem region).

Martin


>
>
> 	gr <- ampRng[j,]
>> gr
> GRanges with 1 range and 1 elementMetadata col:
>          seqnames               ranges strand |   Symbol
>             <Rle>            <IRanges>  <Rle> | <factor>
>    amp11     chr8 [13167119, 13167317]      * |    Lamp1
>    ---
>    seqlengths:
>      chr1 chr10 chr11 chr18  chr2  chr3  chr4  chr6  chr8
>        NA    NA    NA    NA    NA    NA    NA    NA    NA
>
>
>         bam <- "CleanedBams/2424-N1-2.bam"
>   	pu <- pileupAsGRanges(bam, region=gr)
>
>> head(pu)
> GRanges with 6 ranges and 7 elementMetadata cols:
>        seqnames               ranges strand |         A         C         G
>           <Rle>            <IRanges>  <Rle> | <integer> <integer> <integer>
>    [1]     chr8 [13167119, 13167119]      + |         0       176         0
>    [2]     chr8 [13167120, 13167120]      + |         0       206         0
>    [3]     chr8 [13167121, 13167121]      + |         0         1         0
>    [4]     chr8 [13167122, 13167122]      + |         0         0         0
>    [5]     chr8 [13167123, 13167123]      + |         0       251         0
>    [6]     chr8 [13167124, 13167124]      + |         0         3         0
>                T         N     depth                       bam
>        <integer> <integer> <numeric>               <character>
>    [1]         1         0       177 CleanedBams/2424-N1-2.bam
>    [2]         0         0       206 CleanedBams/2424-N1-2.bam
>    [3]       248         0       249 CleanedBams/2424-N1-2.bam
>    [4]       251         0       251 CleanedBams/2424-N1-2.bam
>    [5]         2         0       253 CleanedBams/2424-N1-2.bam
>    [6]       251         0       254 CleanedBams/2424-N1-2.bam
>    ---
>    seqlengths:
>     chr8
>       NA
>
> But if I try and get the pileup of the same region, but a different
> bam file, R just hangs forever and never produces any output. Nor does
> it give me an error message so I'm not sure what it's actually doing.
>
>> bam <- "CleanedBams/2424-N1-1.bam"
>
>> pu <- pileupAsGRanges(bam, region=gr)
>
>
> Any idea how I can debug the problem.? I am able to run samtools
> pileup at the command line for the offending bam file without any
> problem.
>
>
> samtools mpileup -r chr8:13167119-13167317 CleanedBams/2424-N1-1.bam | head
>
> [mpileup] 1 samples in 1 input files
> <mpileup> Set max per-file depth to 8000
> chr8	13167119	N	96	^>C^]C^>C^>C^>C^2C^]C^>C^2T^>C^]C^>C^]C^]C^>C^>C^>T^]C^>C^>C^>C^>C^>C^>C^]C^>C^>C^>C^>C^]C^]C^]C^>C^2C^>C^>C^]C^>C^]C^>C^]C^]C^>C^>C^>C^>T^>C^]T^]C^>C^]C^>C^]C^]C^]C^]C^]C^>C^]T^>T^>C^]C^]C^]C^]C^>T^>C^>C^>C^>C^>C^]C^>T^>C^]C^]C^]T^]C^]C^]C^]C^]C^>C^]C^>C^>C^>C^]C^2T^>C^>C^>C^]C^>C^]C^]T	)$.,),*.$,)$-,.)$)-)--,.,.,,-,,,,,*.,.,,..$)-$-$,-$.,,,,,+$$$,,.,$$,$...*-++$$-..-,+),-*$-$..$-$
> chr8	13167120	N	96	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC	-0++0+0-6+-0-+--;-.2-+*++**-+*+*++-*+9)+++--06.1+--9+++++);;0++++10+0.9+6-));0+-+2*)2-02120-.0-6
> chr8	13167121	N	96	TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT	666686;;;66;666;;;;=;6666;;66;6;666;6;36;6A6;=;63;;;666363=;;66;6;86;6=6;633=;666=66A68;;;;;=86;
> chr8	13167122	N	96	TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTTTTTTTTTTTTTTTTTTTTTT	777777747377337788787(777477773777778837748777773777377383877377378772877*337777787373782747-777
> chr8	13167123	N	96	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC	88881888:888888718888888=888888888888888:88888888888888888888888888888888888=8=88888=88=88888888
> chr8	13167124	N	96	TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT	<=@=:<<;<=><==<::><<;><<=<<7><=<<<@<>;=<<<><<><==<;<=<<<7=<=<><;<<<=<:<=<;==<<<;=>==<<<=<<:<;><=
> chr8	13167125	N	96	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCACCCCCCCCCCCCCCCCCCCCCCCCCCC	9<:<78:::;:0;:::86:73::6;5:38:::::;;:8;:8:8:7:9;;:35;:79:;:::::8;:6;93:;:8;;7:8:;8:;3::97:8383<;
> chr8	13167126	N	96	AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAAAAAAAAAAAAAA	2:::-:)9:=::==44:6:::::6::::::=::-::::=:::84:8)==:::=::::===9:2::::=+92*86==:84=:5:=4;8:*468884=
> chr8	13167127	N	96	TTTTTTTTTTTTTT-1NTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT;9<9=.95<>>89<..>9>>>;:><9<<;9899<8;;<8><99.8<59999>9<9:99:9<9;99:99:<<:;;98><;<9>9959>;.<:==:>9
> chr8	13167128	N	96	GGGGGGGGGGGGG*GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG	<<<<<<8-<<<8<<6<4-;<;<<<;<<<:<<<<6<<6;<6<<6<;<<<<<<<:6<<47<<<<<<<<6<8;<<<-<<<<8<<::<6<<6<<<<86;<
>
>
> Regards,
>
> Mark
>
>
>> sessionInfo()
> R version 2.15.1 (2012-06-22)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>   [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>   [7] LC_PAPER=C                 LC_NAME=C
>   [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] biovizBase_1.4.2     Rsamtools_1.8.6      Biostrings_2.24.1
> [4] GenomicRanges_1.8.13 IRanges_1.14.4       BiocGenerics_0.2.0
>
> loaded via a namespace (and not attached):
>   [1] AnnotationDbi_1.18.3  Biobase_2.16.0        biomaRt_2.12.0
>   [4] bitops_1.0-4.1        BSgenome_1.24.0       cluster_1.14.2
>   [7] colorspace_1.1-1      DBI_0.2-5             dichromat_1.2-4
> [10] GenomicFeatures_1.8.3 grid_2.15.1           Hmisc_3.9-3
> [13] labeling_0.1          lattice_0.20-10       munsell_0.4
> [16] plyr_1.7.1            RColorBrewer_1.0-5    RCurl_1.91-1
> [19] RSQLite_0.11.2        rtracklayer_1.16.3    scales_0.2.2
> [22] stats4_2.15.1         stringr_0.6.1         tools_2.15.1
> [25] XML_3.9-4             zlibbioc_1.2.0
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>


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



More information about the Bioconductor mailing list