.wig files for strand-specific paired-end RNA-Seq
2
0
Entering edit mode
Guest User ★ 13k
@guest-user-4897
Last seen 9.6 years ago
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.
• 1.9k views
ADD COMMENT
0
Entering edit mode
@michael-lawrence-3846
Last seen 2.3 years ago
United States
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]]
ADD COMMENT
0
Entering edit mode
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. > >
ADD REPLY
0
Entering edit mode
@michael-lawrence-3846
Last seen 2.3 years ago
United States
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]]
ADD COMMENT

Login before adding your answer.

Traffic: 509 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6