The tallyVariants function allows binned counts by read position. So I
think all you need to do here is create a bin for the first position [1,1]
and consider only the counts from that. The only down-side is that the
counts are not quality filtered. Unfortunately, there is no cross-tallying
of read position and quality, so there's no way to get this without a
significant change to bam_tally. Also note that tallyVariants only reports
positions with at least one alt read; not sure if this matters for your use
case. You can turn this off if you use the low-level gmapR interface to
bam_tally.

Michael


On Thu, Jul 11, 2013 at 1:24 PM, Valerie Obenchain <vobencha@fhcrc.org>wrote:

> On 07/11/2013 01:11 PM, Valerie Obenchain wrote:
>
>> As an fyi, Martin just checked in a change to Rsamtools (devel) that
>> allows filtering by range. This should speed up the filtering
>> substantially. I'll use the new syntax in the filter example below.
>>
>> I'm not sure of the best way to tallyVariants for all groups of reads
>> for each unique start position. Here's an approach I might take - maybe
>> others will chime in with their ideas. Start by reading in the bam and
>> identify the unique starts.
>>
>>    library(Rsamtools)
>>    fl <- system.file("extdata", "ex1.bam", package="Rsamtools")
>>    gal <- readGAlignments(fl)
>>    starts <- unique(start(gal))
>>
>> mclapply() over the seqlevels containing unique starts.
>> (code not tested)
>>
>
> Nope. This won't work. Trying to be too clever with the mclapply().
>
> You could create two vectors, one of seqlevels and one of starts (unique
> combinations). Pass these to Map() and use the same logic as below. This
> would allow you to create the correct param, filter and get the tally
> results for each seqlevel-start combo.
>
> Valerie
>
>
>>    lst <- split(start(gal), as.factor(seqnames(gal)))
>>    starts <- lapply(lst, unique)
>>    vtparam <- VariantTallyParam(...)
>>    mclapply(starts,
>>             function(start, rname, fl, vtparam) {
>>               ## 'what' imports only the fields used in the filter
>>               ## 'which' is the position(s) to overlap
>>               param <- ScanBamParam(
>>                                 what=c("rname", "pos"),
>>                                 which=GRanges(rname,
>>                                   IRanges(start, width=1)))
>>               filter <- FilterRules(list(atPos = function(x) {
>>                           (x$rname %in% rname) &
>>                           (x$pos %in% start)}))
>>               dest0 <- filterBam(fl, tempfile(),
>>                                  filter=filter, param=param)
>>               tallyVariants(dest0, vtparam)},
>>               rname=names(starts), fl=fl, vtparam=vtparam)
>>
>>
>> Rsamtools in devel is broken today but should be ok tomorrow.
>>
>> Valerie
>>
>> On 07/11/2013 10:27 AM, Taylor, Sean D wrote:
>>
>>> Thanks Valerie, I'll give that a try. And I see that you just
>>> clarified my question about the input file for tallyVariants.
>>>
>>> I had thought about doing this before but hadn't been able to figure
>>> out the proper filter syntax. Thanks!  The only downside to this
>>> approach is that I ultimately want to do this for all unique start
>>> sites that align to my reference. That could entail generating
>>> thousands of tempfiles that all have to be read back in. It could be
>>> done with lapply, but it sounds rather memory and time intensive.
>>> Another approach that I tried was reading the bam file as a
>>> GappedAlignments object and then splitting it by start position into a
>>> GAlignmentsList. That was really easy, but so far as I can tell
>>> tallyVariants will not accept GappedAlignments objects as an input. I
>>> wonder if there is another way to vectorize this approach?
>>>
>>> Thanks,
>>> Sean
>>>
>>>
>>> -----Original Message-----
>>> From: Valerie Obenchain [mailto:vobencha@fhcrc.org]
>>> Sent: Thursday, July 11, 2013 10:02 AM
>>> To: Taylor, Sean D
>>> Cc: bioconductor@r-project.org; Michael Lawrence
>>> Subject: Re: [BioC] Filtering BAM files by start position for
>>> VariantTools
>>>
>>> Hi Sean,
>>>
>>> As you've discovered, the 'which' in the 'param' (for reading bam
>>> files) specifies positions to overlap, not start or end. One approach
>>> to isolating reads with a specific starting position would be to
>>> filter the bam by 'pos'.
>>>
>>>       library(VariantTools)
>>>       fl <- LungCancerLines::**LungCancerBamFiles()$H1993
>>>
>>>       mystart <- 1110426
>>>       filt <- list(setStart=function(x) x$pos %in% mystart)
>>>       dest <- tempfile()
>>>       filterBam(fl, dest, filter=FilterRules(filt))
>>>       scn <- scanBam(dest)
>>>
>>> Confirm all reads start with 'mystart':
>>>
>>>   > table(scn[[1]]$pos)
>>>
>>> 1110426
>>>      2388
>>>
>>> If you want a tally of all nucleotides for all sequences starting with
>>> 'mystart' then no need to supply a 'which':
>>>       param <- VariantTallyParam(gmapR::**TP53Genome(),
>>>                                  readlen=100L,
>>>                                  high_base_quality=23L)
>>>       tally <- tallyVariants(fl, param)
>>>
>>>
>>> Valerie
>>>
>>>
>>> On 07/09/2013 02:06 PM, Taylor, Sean D wrote:
>>>
>>>> I am trying to read a specific set of records from a bam file for use
>>>> in the VariantTools package. I'm trying to construct a which argument
>>>> (a GRanges object) that will pull in a set of records from all reads
>>>> that only start at a specified position. (i.e. all reads that start
>>>> at position 100). So far I have only been able to specify reads that
>>>> overlap position 100, but have not been able to find a way to define
>>>> the start site.
>>>>
>>>> #Example code:
>>>>
>>>>> bams <- LungCancerLines::**LungCancerBamFiles()
>>>>> bam <- bams$H1993
>>>>> which<-GRanges(seqnames=c('**TP53'), IRanges(1110426, width=1), '+')
>>>>> tally.param <- VariantTallyParam(gmapR::**TP53Genome(),
>>>>>
>>>> + readlen = 100L,
>>>> + high_base_quality = 23L,
>>>> + which = which)
>>>>
>>>>> raw.variants <- tallyVariants(bam, tally.param)
>>>>>
>>>>
>>>> This code shows all the variants at position 1110426, but not all the
>>>> variants from the reads that start at position 1110426.
>>>>
>>>> Ultimately, I am trying to do this for all start positions in my data
>>>> set, so I would want something that looks like this pseudocode:
>>>>
>>>>> raw.variants<-lapply (start(bam), function (x){
>>>>>
>>>>     which<-GRanges(seqnames=c('**chrM'), '+', start=x)
>>>>     tally.param<-**VariantTallyParam(gmap, readlen=100L, which=which)
>>>>     tallyVariants(bamfile, tally.param)
>>>> })
>>>>
>>>> Thanks,
>>>> Sean
>>>>
>>>>  sessionInfo()
>>>>>
>>>> R version 3.0.1 (2013-05-16)
>>>> Platform: x86_64-unknown-linux-gnu (64-bit)
>>>>
>>>> locale:
>>>> [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>>>> LC_TIME=en_US.UTF-8
>>>>    [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8
>>>> LC_MESSAGES=en_US.UTF-8
>>>>    [7] LC_PAPER=C                 LC_NAME=C
>>>> LC_ADDRESS=C
>>>> [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8
>>>> LC_IDENTIFICATION=C
>>>>
>>>> attached base packages:
>>>> [1] grid      parallel  stats     graphics  grDevices utils
>>>> datasets  methods   base
>>>>
>>>> other attached packages:
>>>> [1] RSQLite_0.11.4          DBI_0.2-7               BiocInstaller_1.10.2
>>>>    [4] LungCancerLines_0.0.8   GenomicFeatures_1.12.1
>>>> AnnotationDbi_1.22.5
>>>>    [7] Biobase_2.20.0          gmapR_1.2.0
>>>> latticeExtra_0.6-24
>>>> [10] lattice_0.20-15         RColorBrewer_1.0-5      genoPlotR_0.8
>>>> [13] ade4_1.5-2              VariantTools_1.2.2
>>>> VariantAnnotation_1.6.6
>>>> [16] Rsamtools_1.12.3        Biostrings_2.28.0
>>>> GenomicRanges_1.12.4
>>>> [19] IRanges_1.18.1          BiocGenerics_0.6.0
>>>>
>>>> loaded via a namespace (and not attached):
>>>> [1] biomaRt_2.16.0                          bitops_1.0-5
>>>>    [3] BSgenome_1.28.0
>>>> BSgenome.Hsapiens.UCSC.hg19_1.**3.19
>>>>    [5] graph_1.38.2                            Matrix_1.0-12
>>>>    [7] org.Hs.eg.db_2.9.0                      RBGL_1.36.2
>>>>    [9] RCurl_1.95-4.1                          rtracklayer_1.20.2
>>>> [11] stats4_3.0.1                            tools_3.0.1
>>>> [13] TxDb.Hsapiens.UCSC.hg19.**knownGene_2.9.2 XML_3.96-1.1
>>>> [15] zlibbioc_1.6.0
>>>>
>>>>
>>>>     [[alternative HTML version deleted]]
>>>>
>>>> ______________________________**_________________
>>>> 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>
>>>>
>>>>
>>>
>> ______________________________**_________________
>> 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]]

