Thanks Martin,
that gives me a new lead to explore further.

Marcus


On Wed, Nov 28, 2012 at 11:26 AM, Martin Morgan <mtmorgan@fhcrc.org> wrote:

> On 11/27/2012 01:17 PM, Marcus Davy wrote:
>
>> Hi,
>> I would like to export certain reads filtered using Rsamtools, for example
>> aligning against contaminant references
>> to exclude a set of sequences filtering on the bitwise flag=0x4 (similar
>> to
>> SamToFastq.jar in the picard http://picard.sourceforge.net/ tools),
>> but couldn't find much on exporting/writing to file in the
>> Rsamtools-Overview.pdf vignette or on google.
>>
>> Is there a function/utility which converts an Rsamtools bam list region to
>> ShortReadQ object(s) for export using writeFastq()
>> with single end *and* paired end alignments which have been read in via
>> the
>> Rsamtools interface?
>>
>> ## Load a bam
>> bam <- scanBam(...)
>>
>> ## Filter using R's subsetting,
>> bamRegion <- bam[[1]][ ... ] ## list of list
>>
>> ## Should work for single end alignments, but not paired end alignments
>> ## Read pairs may have to index on mpos = PNEXT field
>> rfq <- ShortReadQ(bamRegion[["seq"]], bamRegion[["qual"]],
>> BStringSet(bamRegion[["qname"]**]))
>>
>> If not is the area of exporting from Rsamtools worth developing futher?
>>
>
> I think I would use readBamGappedAlignments, specifying what=c("seq",
> "qual") and maybe "qname", flag=scanBamFlag(<...>) and perhaps the tag
> argument to filter / load relevant records, and then create a fastq object
> ShortReadQ(<...>) and writeFastq.
>
> This could be made space efficient using BamFile(, yieldSize=) and
> iterating through
>
>   bf = open(BamFile(...))
>   out = FastqFile("..."
>   while (nrow(ga <- readBamGappedAlignments(bf, ...))) {
>      ## filter, create fq = ShortReadQ
>      writeFastq(fq, out, mode ="a")
>   }
>   close(bf)
>   close(out)
>
>
>> cheers,
>>
>> Marcus
>>
>>         [[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>
>>
>>
>
> --
> 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
>

	[[alternative HTML version deleted]]

