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? 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.

On Wed, Jun 20, 2012 at 4:53 AM, Igor Ulitsky [guest] < guest@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. > [[alternative HTML version deleted]]

As long as you change the strand assignment to the XS values, it should be what you want. Michael On Wed, Jun 20, 2012 at 6:06 AM, Igor Ulitsky <ulitskyi@gmail.com> wrote: > 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@gene.com> wrote: > > > > > > On Wed, Jun 20, 2012 at 4:53 AM, Igor Ulitsky [guest] > > <guest@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. > > > > > [[alternative HTML version deleted]]