junctions-methods: Extract junctions from genomic alignments
In GenomicAlignments: Representation and manipulation of short genomic alignments
Description
Usage
Arguments
Details
Value
Author(s)
References
See Also
Examples
Description
Given an object x containing genomic alignments (e.g. a
GAlignments, GAlignmentPairs, or GAlignmentsList
object), junctions(x) extracts the junctions from it and
summarizeJunctions(x) extracts and summarizes them.
readTopHatJunctions and readSTARJunctions are utilities
for importing the junction file generated by the TopHat and STAR
aligners, respectively.
Usage
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
## junctions() generic and methods
## -------------------------------
junctions(x, use.mcols=FALSE, ...)
## S4 method for signature 'GAlignments'
junctions(x, use.mcols=FALSE)
## S4 method for signature 'GAlignmentPairs'
junctions(x, use.mcols=FALSE)
## S4 method for signature 'GAlignmentsList'
junctions(x, use.mcols=FALSE, ignore.strand=FALSE)
## summarizeJunctions() and NATURAL_INTRON_MOTIFS
## ----------------------------------------------
summarizeJunctions(x, with.revmap=FALSE, genome=NULL)
NATURAL_INTRON_MOTIFS
## Utilities for importing the junction file generated by some aligners
## --------------------------------------------------------------------
readTopHatJunctions(file, file.is.raw.juncs=FALSE)
readSTARJunctions(file)
Arguments
x
A GAlignments, GAlignmentPairs, or GAlignmentsList
object.
use.mcols
TRUE or FALSE (the default).
Whether the metadata columns on x (accessible with mcols(x))
should be propagated to the returned object or not.
...
Additional arguments, for use in specific methods.
ignore.strand
TRUE or FALSE (the default).
If set to TRUE, then the strand of x is set to "*"
prior to any computation.
with.revmap
TRUE or FALSE (the default).
If set to TRUE, then a revmap metadata column is added to
the output of summarizeJunctions. This metadata column is an
IntegerList object representing the mapping from each
element in the ouput (i.e. each junction) to the corresponding elements in
the input x.
genome
NULL (the default), or a BSgenome object
containing the sequences of the reference genome that was used to align
the reads, or the name of this reference genome specified in a way that
is accepted by the getBSgenome function defined
in the BSgenome software package. In that case the corresponding
BSgenome data package needs to be already installed (see
?getBSgenome in the BSgenome package for
the details).
If genome is supplied, then the intron_motif and
intron_strand metadata columns are computed (based on the
dinucleotides found at the intron boundaries) and added to the output
of summarizeJunctions. See the Value section below for a
description of these metadata columns.
file
The path (or a connection) to the junction file generated by the aligner.
This file should be the junctions.bed or new_list.juncs
file for readTopHatJunctions, and the SJ.out.tab file for
readSTARJunctions.
file.is.raw.juncs
TRUE or FALSE (the default).
If set to TRUE, then the input file is assumed to be a TopHat
.juncs file instead of the junctions.bed file generated
by TopHat. A TopHat .juncs file can be obtained by passing the
junctions.bed file thru TopHat's bed_to_juncs script.
See the TopHat manual at http://tophat.cbcb.umd.edu/manual.shtml
for more information.
Details
An N operation in the CIGAR of a genomic alignment is interpreted as a
junction. junctions(x) will return the genomic ranges of all
junctions found in x.
More precisely, on a GAlignments object x,
junctions(x) is equivalent to:
1
psetdiff(granges(x), grglist(x, order.as.in.query=TRUE))
On a GAlignmentPairs object x, it's equivalent to (but
faster than):
1
2
Note that starting with BioC 3.2, the behavior of junctions on a
GAlignmentPairs object has been slightly modified so that the
returned ranges now have the real strand set on them. See the
documentation of the real.strand argument in the man page of
GAlignmentPairs objects for more information.
NATURAL_INTRON_MOTIFS is a predefined character vector containing
the 5 natural intron motifs described at
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC84117/.
Value
junctions(x) returns the genomic ranges of the junctions in a
GRangesList object parallel to x
(i.e. with 1 list element per element in x).
If x has names on it, they're propagated to the returned object.
If use.mcols is TRUE and x has metadata columns on it
(accessible with mcols(x)), they're propagated to the returned object.
summarizeJunctions returns the genomic ranges of the unique
junctions in x in an unstranded GRanges
object with the following metadata columns:
-
score: The total number of alignments crossing each
junction, i.e. that have the junction encoded in their CIGAR.
-
plus_score and minus_score: The strand-specific
number of alignments crossing each junction.
-
revmap: [Only if with.revmap was set to TRUE.]
An IntegerList object representing the mapping from
each element in the ouput (i.e. each junction) to the corresponding
elements in input x.
-
intron_motif and intron_strand: [Only if genome
was supplied.] The intron motif and strand for each junction,
based on the dinucleotides found in the genome sequences at the
intron boundaries.
The intron_motif metadata column is a factor whose levels
are the 5 natural intron motifs stored in predefined character
vector NATURAL_INTRON_MOTIFS.
If the dinucleotides found at the intron boundaries don't match
any of these natural intron motifs, then intron_motif and
intron_strand are set to NA and *,
respectively.
readTopHatJunctions and readSTARJunctions return the
junctions reported in the input file in a stranded
GRanges object.
With the following metadata columns for readTopHatJunctions (when
reading in the junctions.bed file):
-
name: An id assigned by TopHat to each junction. This id is
of the form JUNC00000017 and is unique within the
junctions.bed file.
-
score: The total number of alignments crossing each
junction.
With the following metadata columns for readSTARJunctions:
-
intron_motif and intron_strand:
The intron motif and strand for each junction, based on the code
found in the input file (0: non-canonical, 1: GT/AG, 2: CT/AC,
3: GC/AG, 4: CT/GC, 5: AT/AC, 6: GT/AT).
Note that of the 5 natural intron motifs stored in predefined
character vector NATURAL_INTRON_MOTIFS, only the first 3 are
assigned codes by the STAR software (2 codes per motif, one if the
intron is on the plus strand and one if it's on the minus strand).
Thus the intron_motif metadata column is a factor with only
3 levels. If code is 0, then intron_motif and
intron_strand are set to NA and *, respectively.
-
um_reads: The number of uniquely mapping reads crossing
the junction (a pair where the 2 alignments cross the same junction
is counted only once).
-
mm_reads: The number of multi-mapping reads crossing the
junction (a pair where the 2 alignments cross the same junction is
counted only once).
See STAR manual at https://code.google.com/p/rna-star/ for more
information.
Author(s)
Hervé Pagès
References
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC84117/ for the 5 natural
intron motifs stored in predefined character vector
NATURAL_INTRON_MOTIFS.
TopHat2: accurate alignment of transcriptomes in the presence of insertions,
deletions and gene fusions
TopHat2 paper: http://genomebiology.com/2013/14/4/r36
TopHat2 software and manual: http://tophat.cbcb.umd.edu/
STAR: ultrafast universal RNA-seq aligner
STAR paper: http://bioinformatics.oxfordjournals.org/content/early/2012/10/25/bioinformatics.bts635
STAR software and manual: https://code.google.com/p/rna-star/
See Also
The readGAlignments and
readGAlignmentPairs functions for reading genomic
alignments from a BAM file.
-
GAlignments, GAlignmentPairs, and
GAlignmentsList objects.
-
GRanges and GRangesList
objects implemented and documented in the GenomicRanges
package.
-
IntegerList objects implemented and documented
in the IRanges package.
The getBSgenome function in the
BSgenome package, for searching the installed BSgenome
data packages for the specified genome and returning it as a
BSgenome object.
The extractList function in the IRanges
package, for extracting groups of elements from a vector-like object
and returning them into a List object.
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
library(RNAseqData.HNRNPC.bam.chr14)
bamfile <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
## ---------------------------------------------------------------------
## A. junctions()
## ---------------------------------------------------------------------
gal <- readGAlignments(bamfile)
table(njunc(gal)) # some alignments have 3 junctions!
juncs <- junctions(gal)
juncs
stopifnot(identical(unname(elementNROWS(juncs)), njunc(gal)))
galp <- readGAlignmentPairs(bamfile)
juncs <- junctions(galp)
juncs
stopifnot(identical(unname(elementNROWS(juncs)), njunc(galp)))
## ---------------------------------------------------------------------
## B. summarizeJunctions()
## ---------------------------------------------------------------------
## By default, only the "score", "plus_score", and "minus_score"
## metadata columns are returned:
junc_summary <- summarizeJunctions(gal)
junc_summary
## The "score" metadata column reports the total number of alignments
## crossing each junction, i.e., that have the junction encoded in their
## CIGAR:
median(mcols(junc_summary)$score)
## The "plus_score" and "minus_score" metadata columns report the
## strand-specific number of alignments crossing each junction:
stopifnot(identical(mcols(junc_summary)$score,
mcols(junc_summary)$plus_score +
mcols(junc_summary)$minus_score))
## If 'with.revmap' is TRUE, the "revmap" metadata column is added to
## the output. This metadata column is an IntegerList object represen-
## ting the mapping from each element in the ouput (i.e. a junction) to
## the corresponding elements in the input 'x'. Here we're going to use
## this to compute a 'score2' for each junction. We obtain this score
## by summing the mapping qualities of the alignments crossing the
## junction:
gal <- readGAlignments(bamfile, param=ScanBamParam(what="mapq"))
junc_summary <- summarizeJunctions(gal, with.revmap=TRUE)
junc_score2 <- sum(extractList(mcols(gal)$mapq,
mcols(junc_summary)$revmap))
mcols(junc_summary)$score2 <- junc_score2
## If the name of the reference genome is specified thru the 'genome'
## argument (in which case the corresponding BSgenome data package needs
## to be installed), then summarizeJunctions() returns the intron motif
## and strand for each junction.
## Since the reads in RNAseqData.HNRNPC.bam.chr14 were aligned to
## the hg19 genome, the following requires that you have
## BSgenome.Hsapiens.UCSC.hg19 installed:
junc_summary <- summarizeJunctions(gal, with.revmap=TRUE, genome="hg19")
mcols(junc_summary)$score2 <- junc_score2 # putting 'score2' back
## The "intron_motif" metadata column is a factor whose levels are the
## 5 natural intron motifs stored in predefined character vector
## 'NATURAL_INTRON_MOTIFS':
table(mcols(junc_summary)$intron_motif)
## ---------------------------------------------------------------------
## C. STRANDED RNA-seq PROTOCOL
## ---------------------------------------------------------------------
## Here is a simple test for checking whether the RNA-seq protocol was
## stranded or not:
strandedTest <- function(plus_score, minus_score)
(sum(plus_score ^ 2) + sum(minus_score ^ 2)) /
sum((plus_score + minus_score) ^ 2)
## The result of this test is guaranteed to be >= 0.5 and <= 1.
## If, for each junction, the strand of the crossing alignments looks
## random (i.e. "plus_score" and "minus_score" are close), then
## strandedTest() will return a value close to 0.5. If it doesn't look
## random (i.e. for each junction, one of "plus_score" and "minus_score"
## is much bigger than the other), then strandedTest() will return a
## value close to 1.
## If the reads are single-end, the test is meaningful when applied
## directly on 'junc_summary'. However, for the test to be meaningful
## on paired-end reads, it needs to be applied on the first and last
## alignments separately:
junc_summary1 <- summarizeJunctions(first(galp))
junc_summary2 <- summarizeJunctions(last(galp))
strandedTest(mcols(junc_summary1)$plus_score,
mcols(junc_summary1)$minus_score)
strandedTest(mcols(junc_summary2)$plus_score,
mcols(junc_summary2)$minus_score)
## Both values are close to 0.5 which suggests that the RNA-seq protocol
## used for this experiment was not stranded.
## ---------------------------------------------------------------------
## D. UTILITIES FOR IMPORTING THE JUNCTION FILE GENERATED BY SOME
## ALIGNERS
## ---------------------------------------------------------------------
## The TopHat aligner generates a junctions.bed file where it reports
## all the junctions satisfying some "quality" criteria (see the TopHat
## manual at http://tophat.cbcb.umd.edu/manual.shtml for more
## information). This file can be loaded with readTopHatJunctions():
runname <- names(RNAseqData.HNRNPC.bam.chr14_BAMFILES)[1]
junctions_file <- system.file("extdata", "tophat2_out", runname,
"junctions.bed",
package="RNAseqData.HNRNPC.bam.chr14")
th_junctions <- readTopHatJunctions(junctions_file)
## Comparing the "TopHat junctions" with the result of
## summarizeJunctions():
th_junctions14 <- th_junctions
seqlevels(th_junctions14, pruning.mode="coarse") <- "chr14"
mcols(th_junctions14)$intron_strand <- strand(th_junctions14)
strand(th_junctions14) <- "*"
## All the "TopHat junctions" are in 'junc_summary':
stopifnot(all(th_junctions14 %in% junc_summary))
## But not all the junctions in 'junc_summary' are reported by TopHat
## (that's because TopHat reports only junctions that satisfy some
## "quality" criteria):
is_in_th_junctions14 <- junc_summary %in% th_junctions14
table(is_in_th_junctions14) # 32 junctions are not in TopHat's
# junctions.bed file
junc_summary2 <- junc_summary[is_in_th_junctions14]
## 'junc_summary2' and 'th_junctions14' contain the same junctions in
## the same order:
stopifnot(all(junc_summary2 == th_junctions14))
## Let's merge their metadata columns. We use our own version of
## merge() for this, which is stricter (it checks that the common
## columns are the same in the 2 data frames to merge) and also
## simpler:
merge2 <- function(df1, df2)
{
common_colnames <- intersect(colnames(df1), colnames(df2))
lapply(common_colnames,
function(colname)
stopifnot(all(df1[ , colname] == df2[ , colname])))
extra_mcolnames <- setdiff(colnames(df2), colnames(df1))
cbind(df1, df2[ , extra_mcolnames, drop=FALSE])
}
mcols(th_junctions14) <- merge2(mcols(th_junctions14),
mcols(junc_summary2))
## Here is a peculiar junction reported by TopHat:
idx0 <- which(mcols(th_junctions14)$score2 == 0L)
th_junctions14[idx0]
gal[mcols(th_junctions14)$revmap[[idx0]]]
## The junction is crossed by 5 alignments (score is 5), all of which
## have a mapping quality of 0!
GenomicAlignments documentation built on Nov. 8, 2020, 8:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) References See Also Examples
Description
Given an object x containing genomic alignments (e.g. a
GAlignments, GAlignmentPairs, or GAlignmentsList
object), junctions(x) extracts the junctions from it and
summarizeJunctions(x) extracts and summarizes them.
readTopHatJunctions and readSTARJunctions are utilities
for importing the junction file generated by the TopHat and STAR
aligners, respectively.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | ## junctions() generic and methods
## -------------------------------
junctions(x, use.mcols=FALSE, ...)
## S4 method for signature 'GAlignments'
junctions(x, use.mcols=FALSE)
## S4 method for signature 'GAlignmentPairs'
junctions(x, use.mcols=FALSE)
## S4 method for signature 'GAlignmentsList'
junctions(x, use.mcols=FALSE, ignore.strand=FALSE)
## summarizeJunctions() and NATURAL_INTRON_MOTIFS
## ----------------------------------------------
summarizeJunctions(x, with.revmap=FALSE, genome=NULL)
NATURAL_INTRON_MOTIFS
## Utilities for importing the junction file generated by some aligners
## --------------------------------------------------------------------
readTopHatJunctions(file, file.is.raw.juncs=FALSE)
readSTARJunctions(file)
|
Arguments
x |
A GAlignments, GAlignmentPairs, or GAlignmentsList object. |
use.mcols |
|
... |
Additional arguments, for use in specific methods. |
ignore.strand |
|
with.revmap |
|
genome |
If |
file |
The path (or a connection) to the junction file generated by the aligner.
This file should be the junctions.bed or new_list.juncs
file for |
file.is.raw.juncs |
|
Details
An N operation in the CIGAR of a genomic alignment is interpreted as a
junction. junctions(x) will return the genomic ranges of all
junctions found in x.
More precisely, on a GAlignments object x,
junctions(x) is equivalent to:
1 | psetdiff(granges(x), grglist(x, order.as.in.query=TRUE))
|
On a GAlignmentPairs object x, it's equivalent to (but
faster than):
1 2 |
Note that starting with BioC 3.2, the behavior of junctions on a
GAlignmentPairs object has been slightly modified so that the
returned ranges now have the real strand set on them. See the
documentation of the real.strand argument in the man page of
GAlignmentPairs objects for more information.
NATURAL_INTRON_MOTIFS is a predefined character vector containing
the 5 natural intron motifs described at
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC84117/.
Value
junctions(x) returns the genomic ranges of the junctions in a
GRangesList object parallel to x
(i.e. with 1 list element per element in x).
If x has names on it, they're propagated to the returned object.
If use.mcols is TRUE and x has metadata columns on it
(accessible with mcols(x)), they're propagated to the returned object.
summarizeJunctions returns the genomic ranges of the unique
junctions in x in an unstranded GRanges
object with the following metadata columns:
-
score: The total number of alignments crossing each junction, i.e. that have the junction encoded in their CIGAR. -
plus_scoreandminus_score: The strand-specific number of alignments crossing each junction. -
revmap: [Only ifwith.revmapwas set toTRUE.] An IntegerList object representing the mapping from each element in the ouput (i.e. each junction) to the corresponding elements in inputx. -
intron_motifandintron_strand: [Only ifgenomewas supplied.] The intron motif and strand for each junction, based on the dinucleotides found in the genome sequences at the intron boundaries. Theintron_motifmetadata column is a factor whose levels are the 5 natural intron motifs stored in predefined character vectorNATURAL_INTRON_MOTIFS. If the dinucleotides found at the intron boundaries don't match any of these natural intron motifs, thenintron_motifandintron_strandare set toNAand*, respectively.
readTopHatJunctions and readSTARJunctions return the
junctions reported in the input file in a stranded
GRanges object.
With the following metadata columns for readTopHatJunctions (when
reading in the junctions.bed file):
-
name: An id assigned by TopHat to each junction. This id is of the form JUNC00000017 and is unique within the junctions.bed file. -
score: The total number of alignments crossing each junction.
With the following metadata columns for readSTARJunctions:
-
intron_motifandintron_strand: The intron motif and strand for each junction, based on the code found in the input file (0: non-canonical, 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5: AT/AC, 6: GT/AT). Note that of the 5 natural intron motifs stored in predefined character vectorNATURAL_INTRON_MOTIFS, only the first 3 are assigned codes by the STAR software (2 codes per motif, one if the intron is on the plus strand and one if it's on the minus strand). Thus theintron_motifmetadata column is a factor with only 3 levels. If code is 0, thenintron_motifandintron_strandare set toNAand*, respectively. -
um_reads: The number of uniquely mapping reads crossing the junction (a pair where the 2 alignments cross the same junction is counted only once). -
mm_reads: The number of multi-mapping reads crossing the junction (a pair where the 2 alignments cross the same junction is counted only once).
See STAR manual at https://code.google.com/p/rna-star/ for more information.
Author(s)
Hervé Pagès
References
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC84117/ for the 5 natural
intron motifs stored in predefined character vector
NATURAL_INTRON_MOTIFS.
TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions
TopHat2 paper: http://genomebiology.com/2013/14/4/r36
TopHat2 software and manual: http://tophat.cbcb.umd.edu/
STAR: ultrafast universal RNA-seq aligner
STAR paper: http://bioinformatics.oxfordjournals.org/content/early/2012/10/25/bioinformatics.bts635
STAR software and manual: https://code.google.com/p/rna-star/
See Also
The
readGAlignmentsandreadGAlignmentPairsfunctions for reading genomic alignments from a BAM file.-
GAlignments, GAlignmentPairs, and GAlignmentsList objects.
-
GRanges and GRangesList objects implemented and documented in the GenomicRanges package.
-
IntegerList objects implemented and documented in the IRanges package.
The
getBSgenomefunction in the BSgenome package, for searching the installed BSgenome data packages for the specified genome and returning it as a BSgenome object.The
extractListfunction in the IRanges package, for extracting groups of elements from a vector-like object and returning them into a List object.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | library(RNAseqData.HNRNPC.bam.chr14)
bamfile <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
## ---------------------------------------------------------------------
## A. junctions()
## ---------------------------------------------------------------------
gal <- readGAlignments(bamfile)
table(njunc(gal)) # some alignments have 3 junctions!
juncs <- junctions(gal)
juncs
stopifnot(identical(unname(elementNROWS(juncs)), njunc(gal)))
galp <- readGAlignmentPairs(bamfile)
juncs <- junctions(galp)
juncs
stopifnot(identical(unname(elementNROWS(juncs)), njunc(galp)))
## ---------------------------------------------------------------------
## B. summarizeJunctions()
## ---------------------------------------------------------------------
## By default, only the "score", "plus_score", and "minus_score"
## metadata columns are returned:
junc_summary <- summarizeJunctions(gal)
junc_summary
## The "score" metadata column reports the total number of alignments
## crossing each junction, i.e., that have the junction encoded in their
## CIGAR:
median(mcols(junc_summary)$score)
## The "plus_score" and "minus_score" metadata columns report the
## strand-specific number of alignments crossing each junction:
stopifnot(identical(mcols(junc_summary)$score,
mcols(junc_summary)$plus_score +
mcols(junc_summary)$minus_score))
## If 'with.revmap' is TRUE, the "revmap" metadata column is added to
## the output. This metadata column is an IntegerList object represen-
## ting the mapping from each element in the ouput (i.e. a junction) to
## the corresponding elements in the input 'x'. Here we're going to use
## this to compute a 'score2' for each junction. We obtain this score
## by summing the mapping qualities of the alignments crossing the
## junction:
gal <- readGAlignments(bamfile, param=ScanBamParam(what="mapq"))
junc_summary <- summarizeJunctions(gal, with.revmap=TRUE)
junc_score2 <- sum(extractList(mcols(gal)$mapq,
mcols(junc_summary)$revmap))
mcols(junc_summary)$score2 <- junc_score2
## If the name of the reference genome is specified thru the 'genome'
## argument (in which case the corresponding BSgenome data package needs
## to be installed), then summarizeJunctions() returns the intron motif
## and strand for each junction.
## Since the reads in RNAseqData.HNRNPC.bam.chr14 were aligned to
## the hg19 genome, the following requires that you have
## BSgenome.Hsapiens.UCSC.hg19 installed:
junc_summary <- summarizeJunctions(gal, with.revmap=TRUE, genome="hg19")
mcols(junc_summary)$score2 <- junc_score2 # putting 'score2' back
## The "intron_motif" metadata column is a factor whose levels are the
## 5 natural intron motifs stored in predefined character vector
## 'NATURAL_INTRON_MOTIFS':
table(mcols(junc_summary)$intron_motif)
## ---------------------------------------------------------------------
## C. STRANDED RNA-seq PROTOCOL
## ---------------------------------------------------------------------
## Here is a simple test for checking whether the RNA-seq protocol was
## stranded or not:
strandedTest <- function(plus_score, minus_score)
(sum(plus_score ^ 2) + sum(minus_score ^ 2)) /
sum((plus_score + minus_score) ^ 2)
## The result of this test is guaranteed to be >= 0.5 and <= 1.
## If, for each junction, the strand of the crossing alignments looks
## random (i.e. "plus_score" and "minus_score" are close), then
## strandedTest() will return a value close to 0.5. If it doesn't look
## random (i.e. for each junction, one of "plus_score" and "minus_score"
## is much bigger than the other), then strandedTest() will return a
## value close to 1.
## If the reads are single-end, the test is meaningful when applied
## directly on 'junc_summary'. However, for the test to be meaningful
## on paired-end reads, it needs to be applied on the first and last
## alignments separately:
junc_summary1 <- summarizeJunctions(first(galp))
junc_summary2 <- summarizeJunctions(last(galp))
strandedTest(mcols(junc_summary1)$plus_score,
mcols(junc_summary1)$minus_score)
strandedTest(mcols(junc_summary2)$plus_score,
mcols(junc_summary2)$minus_score)
## Both values are close to 0.5 which suggests that the RNA-seq protocol
## used for this experiment was not stranded.
## ---------------------------------------------------------------------
## D. UTILITIES FOR IMPORTING THE JUNCTION FILE GENERATED BY SOME
## ALIGNERS
## ---------------------------------------------------------------------
## The TopHat aligner generates a junctions.bed file where it reports
## all the junctions satisfying some "quality" criteria (see the TopHat
## manual at http://tophat.cbcb.umd.edu/manual.shtml for more
## information). This file can be loaded with readTopHatJunctions():
runname <- names(RNAseqData.HNRNPC.bam.chr14_BAMFILES)[1]
junctions_file <- system.file("extdata", "tophat2_out", runname,
"junctions.bed",
package="RNAseqData.HNRNPC.bam.chr14")
th_junctions <- readTopHatJunctions(junctions_file)
## Comparing the "TopHat junctions" with the result of
## summarizeJunctions():
th_junctions14 <- th_junctions
seqlevels(th_junctions14, pruning.mode="coarse") <- "chr14"
mcols(th_junctions14)$intron_strand <- strand(th_junctions14)
strand(th_junctions14) <- "*"
## All the "TopHat junctions" are in 'junc_summary':
stopifnot(all(th_junctions14 %in% junc_summary))
## But not all the junctions in 'junc_summary' are reported by TopHat
## (that's because TopHat reports only junctions that satisfy some
## "quality" criteria):
is_in_th_junctions14 <- junc_summary %in% th_junctions14
table(is_in_th_junctions14) # 32 junctions are not in TopHat's
# junctions.bed file
junc_summary2 <- junc_summary[is_in_th_junctions14]
## 'junc_summary2' and 'th_junctions14' contain the same junctions in
## the same order:
stopifnot(all(junc_summary2 == th_junctions14))
## Let's merge their metadata columns. We use our own version of
## merge() for this, which is stricter (it checks that the common
## columns are the same in the 2 data frames to merge) and also
## simpler:
merge2 <- function(df1, df2)
{
common_colnames <- intersect(colnames(df1), colnames(df2))
lapply(common_colnames,
function(colname)
stopifnot(all(df1[ , colname] == df2[ , colname])))
extra_mcolnames <- setdiff(colnames(df2), colnames(df1))
cbind(df1, df2[ , extra_mcolnames, drop=FALSE])
}
mcols(th_junctions14) <- merge2(mcols(th_junctions14),
mcols(junc_summary2))
## Here is a peculiar junction reported by TopHat:
idx0 <- which(mcols(th_junctions14)$score2 == 0L)
th_junctions14[idx0]
gal[mcols(th_junctions14)$revmap[[idx0]]]
## The junction is crossed by 5 alignments (score is 5), all of which
## have a mapping quality of 0!
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.