An alternative approach would be to summarize the Rsamtools::applyPileup()
or gmapR::bam_tally output.

Michael


On Wed, Jun 19, 2013 at 12:45 PM, Hervé Pagès <hpages@fhcrc.org> wrote:

> Hi Jacques,
>
> I just added a tool to the Rsamtools package, that makes this easy.
> It's the sequenceLayer() function. Here is how to use it:
>
> (1) Load the BAM file into a GAlignments object:
>
>       library(Rsamtools)
>       bamfile <- system.file("extdata", "ex1.bam", package="Rsamtools")
>       param <- ScanBamParam(what="seq")
>       gal <- readGAlignmentsFromBam(**bamfile, param=param)
>       qseq <- mcols(gal)$seq  # the query sequences
>
> (2) Use sequenceLayer() to "lay" the query sequences on the reference
>     space. This will remove the parts from the query sequences that
>     correspond to insertions and soft clipping, and it will fill them
>     with - where deletions and/or skipped regions occurred:
>
>       qseq_on_ref <- sequenceLayer(qseq, cigar(gal),
>                                    from="query", to="reference")
>
> (3) Compute 1 consensus matrix per chromosome:
>
>       qseq_on_ref_by_chrom <- splitAsList(qseq_on_ref, seqnames(gal))
>       qseq_pos_by_chrom <- splitAsList(start(gal), seqnames(gal))
>
>       cm_by_chrom <- lapply(names(qseq_pos_by_**chrom),
>           function(seqname)
>               consensusMatrix(qseq_on_ref_**by_chrom[[seqname]],
>                               as.prob=TRUE,
>                               shift=qseq_pos_by_chrom[[**seqname]]-1,
>                               width=seqlengths(gal)[[**seqname]]))
>       names(cm_by_chrom) <- names(qseq_pos_by_chrom)
>
>       ## 'cm_by_chrom' is a list of consensus matrices. Each matrix
>       ## has 17 rows (1 per letter in the DNA alphabet) and 1 column
>       ## per chromosome position.
>
> (4) Compute the consensus string from each consensus matrix. We'll
>     put "+" ins the strings wherever there is no coverage for that
>     position, and "N" where there is coverage but no consensus.
>
>       cs_by_chrom <- lapply(cm_by_chrom,
>           function(cm) {
>               ## Because consensusString() doesn't like consensus
>               ## matrices with columns that contain only zeroes (and
>               ## you will have columns like for chromosome positions
>               ## that don't receive any coverage), we need to "fix"
>               ## 'cm' first.
>               idx <- colSums(cm) == 0
>               cm["+", idx] <- 1
>               DNAString(consensusString(cm, ambiguityMap="N"))
>           })
>
> (5) Write 'cs_by_chrom' to a FASTA file:
>
>       writeXStringSet(DNAStringSet(**cs_by_chrom), "consensus.fa")
>
> Please let us know how this goes with your 454 data.
>
> Note that you will need the 1.13.18 version of Rsamtools for this,
> which is the latest devel version of the package. This means you need
> to install BioC devel. See our website for how to do this.
> Rsamtools 1.13.18 will be available in the next 24 hours or so thru
> biocLite().
>
> Cheers,
> H.
>
>
>
> On 06/19/2013 02:54 AM, Maintainer wrote:
>
>>
>> Hello,
>> I have an indexed bam files from 454 sequencing on short reference
>> sequences. I would like to create a consensus sequence in fasta format from
>> this bam file using R Bioconductor. Is there a way to do that simply ?
>> Thanks in advance
>> Jacques
>>
>>   -- output of sessionInfo():
>>
>> R version 2.15.1 (2012-06-22)
>> Platform: i386-pc-mingw32/i386 (32-bit)
>>
>> locale:
>> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
>> [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
>> [5] LC_TIME=French_France.1252
>>
>> attached base packages:
>> [1] grDevices datasets  splines   graphics  stats     tcltk     utils
>> [8] methods   base
>>
>> other attached packages:
>> [1] Rsamtools_1.10.2     Biostrings_2.26.3    GenomicRanges_1.10.7
>> [4] IRanges_1.16.6       BiocGenerics_0.4.0   TinnR_1.0-5
>> [7] R2HTML_2.2.1         Hmisc_3.10-1.1       survival_2.36-14
>>
>> loaded via a namespace (and not attached):
>> [1] bitops_1.0-5    cluster_1.14.4  grid_2.15.1     lattice_0.20-15
>> [5] parallel_2.15.1 stats4_2.15.1   tools_2.15.1    zlibbioc_1.4.0
>>
>> --
>> Sent via the guest posting facility at bioconductor.org.
>>
>> ______________________________**______________________________**
>> ____________
>> devteam-bioc mailing list
>> To unsubscribe from this mailing list send a blank email to
>> devteam-bioc-leave@lists.**fhcrc.org <devteam-bioc-leave@lists.fhcrc.org>
>> You can also unsubscribe or change your personal options at
>> https://lists.fhcrc.org/**mailman/listinfo/devteam-bioc<https://lists.fhcrc.org/mailman/listinfo/devteam-bioc>
>>
>>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages@fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>
>
> ______________________________**_________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
> Search the archives: http://news.gmane.org/gmane.**
> science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>

	[[alternative HTML version deleted]]

