[BioC] .wig files for strand-specific paired-end RNA-Seq
Igor Ulitsky
ulitskyi at gmail.com
Wed Jun 20 15:06:05 CEST 2012
Yes, I'm trying to get the strand of stranscription. If the XS flags
are correct (I assume they will be if I run tophat with the
appropriate library-type?) then the export will also be on the right
strand?
Thanks!
Igor.
On Wed, Jun 20, 2012 at 3:59 PM, Michael Lawrence
<lawrence.michael at gene.com> wrote:
>
>
> On Wed, Jun 20, 2012 at 4:53 AM, Igor Ulitsky [guest]
> <guest at bioconductor.org> wrote:
>>
>>
>> Hi,
>>
>> Is there a simple way to make strand-specific .wig file (i.e., a separate
>> track for + and - strand) from paired-end data (where the second read maps
>> to the other strand)? I've tried using this:
>>
>> library(Rsamtools)
>> library(rtracklayer)
>> myReads <- readGappedAlignments("RNAseqMapping.bam")
>> coveragePlus <- coverage(myReads[strand(myReads) == '+'])
>> export(coveragePlus, "RNAplus.wig")
>> coverageMinus <- coverage(myReads[strand(myReads) == '-'])
>> export(coverageMinus, "RNAminus.wig")
>>
>> But it appears that the second read in the pair contributes to the other
>> strand, generating similar tracks for the + and the - strands.
>> Is there a way to deal with this better?
>>
>
> Are you trying to generate coverages for the actual strand of transcription?
> If so, you would probably get that information from the XS tag and set it as
> your strand prior to export, but unless you used a special protocol the XS
> information would be incomplete. Btw, I would recommend BigWig export of
> those coverage tracks.
>
> Michael
>
>>
>> Thanks!
>>
>> Igor.
>>
>> -- output of sessionInfo():
>>
>> R version 2.13.1 (2011-07-08)
>> Platform: i386-pc-mingw32/i386 (32-bit)
>>
>> locale:
>> [1] LC_COLLATE=English_United States.1252
>> [2] LC_CTYPE=English_United States.1252
>> [3] LC_MONETARY=English_United States.1252
>> [4] LC_NUMERIC=C
>> [5] LC_TIME=English_United States.1252
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> loaded via a namespace (and not attached):
>> [1] tools_2.13.1
>>
>> --
>> Sent via the guest posting facility at bioconductor.org.
>
>
More information about the Bioconductor
mailing list