Hi Michael,
Thank you very much for prompt reply.

With our ssRNA-Seq protocol, strand info of the second read represents the
whole fragment. Also, even the code can be applied for calculating
coverage, it may not work effectively with large BAM file because of
requiring to load the whole file at once. Is there any trick to avoid this?

Regarding countOverlaps function, as you suggested to convert
GappedAlignmentPairs into a GRangeList corresponding to two reads, it
happens that one fragment may be counted twice. Does it cause any bias on
the output, assuming that the read counts are used for DEG analysis?

Looking forward to your reply.
Best regards,


On Thu, Nov 7, 2013 at 11:32 PM, Michael Lawrence <lawrence.michael@gene.com
> wrote:

>
>
>
> On Thu, Nov 7, 2013 at 5:43 AM, bach le <le.bach.07@gmail.com> wrote:
>
>> Hi,
>> I would like to ask if any of you have experience with generating *WIG
>> coverage files from PE ssRNA-Seq* data. I got some discussion around this
>> (
>>
>> http://seqanswers.com/forums/showthread.php?t=20545) but it seems not to
>> be
>> solved yet.
>>
>
> I've no experience with strand-specific protocols, but it sounds like the
> strand of the first read in the pair indicates the strand. If that is true,
> please confirm as it would inform our design. As suggested in the thread,
> just flip the strand of the second ends and use the code earlier in that
> thread. GAlignmentPairs makes this easy, the accessors first() and last()
> get you the #1 and #2 reads, respectively, as GAlignments.
>
>
>
>> Also, does *countOverlaps* function work with GappedAlignmentPairs object?
>>
>> In that case, will the first or last alignment of each pair be used for
>> calculating the overlap?
>>
>
> GappedAlignmentPairs is now called GAlignmentPairs. Whenever you want to
> know whether there is a method for a certain signature, use showMethods()
> to see the list or selectMethod() to retrieve the method object for the
> signature. In this case:
>
> > showMethods(countOverlaps)
> [snip]
> query="GAlignmentPairs", subject="GAlignmentPairs"
> query="GAlignmentPairs", subject="Vector"
> [snip]
>
> This means we can call countOverlaps with a GAlignmentPairs and any Vector
> (of ranges), including a GRangesList of transcripts. The method works by
> first converting GAlignmentPairs into a GRangesList, where each element
> corresponds to a pair and contains the exonic ranges from both ends. The
> problem for you is that the strand will not be correct for the second end.
> Thus, you will need to call grglist() to get the GRangesList manually and
> modify the strand before passing it to countOverlaps. The strand
> modification would look like:
>
> strand(grl) <- relist(rep(strand(first(reads)), elementLengths(grl)), grl)
>
> There is some experimental functionality in GenomicRanges for computing
> isoform-specific overlaps: findSpliceOverlaps. These functions look for the
> BAM XS tag and if found use it to determine the strand in a way similar to
> above. If you were to set the XS mcol on your GAlignmentPairs to
> strand(first(reads)), then these functions should work the way you want.
>
> Michael
>
> Thanks much in advance.
>> Best,
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>
>

	[[alternative HTML version deleted]]

