coordinate-mapping-methods: Map range coordinates between transcripts and genome space
In jmacdon/GenomicFeatures: Conveniently import and query gene models
Description
Usage
Arguments
Details
Value
Author(s)
See Also
Examples
Description
Map range coordinates between features in the transcriptome and
genome (reference) space.
See ?mapToAlignments in the
GenomicAlignments package for mapping coordinates between
reads (local) and genome (reference) space using a CIGAR alignment.
Usage
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## mapping to transcripts
## S4 method for signature 'GenomicRanges,GenomicRanges'
mapToTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'GenomicRanges,GRangesList'
mapToTranscripts(x, transcripts,
ignore.strand = FALSE, intronJunctions=FALSE)
## S4 method for signature 'ANY,TxDb'
mapToTranscripts(x, transcripts, ignore.strand = FALSE,
extractor.fun = GenomicFeatures::transcripts, ...)
## S4 method for signature 'GenomicRanges,GRangesList'
pmapToTranscripts(x, transcripts,
ignore.strand = FALSE)
## mapping from transcripts
## S4 method for signature 'GenomicRanges,GRangesList'
mapFromTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'GenomicRanges,GRangesList'
pmapFromTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'IntegerRanges,GRangesList'
pmapFromTranscripts(x, transcripts)
Arguments
x
GenomicRanges object of positions to be mapped.
The seqnames of x are used in mapFromTranscripts, i.e.,
when mapping from transcripts to the genome. In the case of
pmapFromTranscripts, x can be an
IntegerRanges object.
transcripts
A named GenomicRanges or
GRangesList object used to map between x
and the result. The ranges can be any feature in the transcriptome
extracted from a TxDb (e.g., introns, exons, cds regions).
See ?transcripts and ?transcriptsBy for a list of
extractor functions.
The transcripts object must have names. When mapping from
transcripts to the genome, they are used to determine mapping pairs; in
the reverse direction they become the seqlevels of the output object.
ignore.strand
When ignore.strand is TRUE, strand is ignored in overlaps operations
(i.e., all strands are considered "+") and the strand in the output is '*'.
When ignore.strand is FALSE strand in the output is taken from the
transcripts argument. When transcripts is a
GRangesList, all inner list elements of a common list element must
have the same strand or an error is thrown.
Mapped position is computed by counting from the transcription start site
(TSS) and is not affected by the value of ignore.strand.
intronJunctions
Logical to indicate if intronic ranges in x should be reported.
This argument is only supported in mapToTranscripts when
transcripts is a GRangesList. When transcripts is a
GRangesList, individual ranges can be thought of as exons and the spaces
between the ranges as introns.
When intronJunctions=TRUE, ranges that fall completely "within" an
intron are reported as a zero-width range (start and end are taken
from the ranges they fall between). A metadata column called "intronic" is
returned with the GRanges and marked as TRUE for these ranges. By
default, intronJunctions=FALSE and these ranges are not mapped.
Ranges that have either the start or end in an intron are considered "non
hits" and are never mapped. Ranges that span introns are always mapped.
Neither of these range types are controlled by the intronJunctions
argument.
extractor.fun
Function to extract genomic features from a TxDb object.
This argument is only applicable to mapToTranscripts when
transcripts is a TxDb object. The extractor should
be the name of a function (not a character()) described on the
?transcripts, ?transcriptsBy, or ?microRNAs man page.
Valid extractor functions:
transcripts ## default
exons
cds
genes
promoters
exonicParts
intronicParts
transcriptsBy
exonsBy
cdsBy
intronsByTranscript
fiveUTRsByTranscript
threeUTRsByTranscript
microRNAs
tRNAs
...
Additional arguments passed to extractor.fun functions.
Details
In GenomicFeatures >= 1.21.10, the default for ignore.strand was
changed to FALSE for consistency with other methods in the
GenomicRanges and GenomicAlignments packages. Additionally,
the mapped position is computed from the TSS and does not depend on the
ignore.strand argument.
See the section on ignore.strand for details.
mapToTranscripts, pmapToTranscripts
The genomic range in x is mapped to the local position in the
transcripts ranges. A successful mapping occurs when x
is completely within the transcripts range, equivalent to:
findOverlaps(..., type="within")
Transcriptome-based coordinates start counting at 1 at the beginning
of the transcripts range and return positions where x
was aligned. The seqlevels of the return object are taken from the
transcripts object and should be transcript names. In this
direction, mapping is attempted between all elements of x and
all elements of transcripts.
mapToTranscripts uses findOverlaps to map ranges in
x to ranges in transcripts. This method does not return
unmapped ranges.
pmapToTranscripts maps the i-th range in x to the
i-th range in transcripts. Recycling is supported for both
x and transcripts when either is length == 1L; otherwise
the lengths must match. Ranges in x that do not map (out of bounds
or strand mismatch) are returned as zero-width ranges starting at 0.
These ranges are given the seqname of "UNMAPPED".
mapFromTranscripts, pmapFromTranscripts
The transcript-based position in x is mapped to genomic coordinates
using the ranges in transcripts. A successful mapping occurs when
the following is TRUE:
width(transcripts) >= start(x) + width(x)
x is aligned to transcripts by moving in start(x)
positions in from the beginning of the transcripts range. The
seqlevels of the return object are chromosome names.
mapFromTranscripts uses the seqname of x and the names
of transcripts to determine mapping pairs (vs attempting to match
all possible pairs). Name matching is motivated by use cases such as
differentially expressed regions where the expressed regions in x
would only be related to a subset of regions in transcripts.
This method does not return unmapped ranges.
pmapFromTranscripts maps the i-th range in x to the i-th
range in transcripts and therefore does not use name matching.
Recycling is supported in pmapFromTranscripts when either
x or transcripts is length == 1L; otherwise the lengths
must match. Ranges in x that do not map (out of bounds or strand
mismatch) are returned as zero-width ranges starting at 0. These ranges
are given the seqname of "UNMAPPED".
Value
pmapToTranscripts returns a GRanges the same length as
x.
pmapFromTranscripts returns a GRanges when transcripts
is a GRanges and a GRangesList when transcripts
is a GRangesList. In both cases the return object is the same
length as x. The rational for returning the GRangesList is
to preserve exon structure; ranges in a list element that are not overlapped
by x are returned as a zero-width range. The GRangesList
return object will have no seqlevels called "UNMAPPED"; those will only
occur when a GRanges is returned.
mapToTranscripts and mapFromTranscripts return GRanges
objects that vary in length similar to a Hits object. The result
contains mapped records only; strand mismatch and out of bound ranges are
not returned. xHits and transcriptsHits metadata columns
(similar to the queryHits and subjectHits of a Hits
object) indicate elements of x and transcripts used in
the mapping.
When intronJunctions is TRUE, mapToTranscripts returns an
extra metdata column named intronic to identify the intron ranges.
When mapping to transcript coordinates, seqlevels of the output are the names
on the transcripts object and most often these will be transcript
names. When mapping to the genome, seqlevels of the output are the seqlevels
of transcripts which are usually chromosome names.
Author(s)
V. Obenchain, M. Lawrence and H. Pagès
See Also
-
?mapToAlignments in the
GenomicAlignments package for methods mapping between
reads and genome space using a CIGAR alignment.
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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
## ---------------------------------------------------------------------
## A. Basic Use: Conversion between CDS and Exon coordinates and the
## genome
## ---------------------------------------------------------------------
## Gene "Dgkb" has ENTREZID "217480":
library(org.Mm.eg.db)
Dgkb_geneid <- get("Dgkb", org.Mm.egSYMBOL2EG)
## The gene is on the positive strand, chromosome 12:
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
tx_by_gene <- transcriptsBy(txdb, by="gene")
Dgkb_transcripts <- tx_by_gene[[Dgkb_geneid]]
Dgkb_transcripts # all 7 Dgkb transcripts are on chr12, positive strand
## To map coordinates from local CDS or exon space to genome
## space use mapFromTranscripts().
## When mapping CDS coordinates to genome space the 'transcripts'
## argument is the collection of CDS regions by transcript.
coord <- GRanges("chr12", IRanges(4, width=1))
## Get the names of the transcripts in the gene:
Dgkb_tx_names <- mcols(Dgkb_transcripts)$tx_name
Dgkb_tx_names
## Use these names to isolate the region of interest:
cds_by_tx <- cdsBy(txdb, "tx", use.names=TRUE)
Dgkb_cds_by_tx <- cds_by_tx[intersect(Dgkb_tx_names, names(cds_by_tx))]
## Dgkb CDS grouped by transcript (no-CDS transcripts omitted):
Dgkb_cds_by_tx
lengths(Dgkb_cds_by_tx) # nb of CDS per transcript
## A requirement for mapping from transcript space to genome space
## is that seqnames in 'x' match the names in 'transcripts'.
names(Dgkb_cds_by_tx) <- rep(seqnames(coord), length(Dgkb_cds_by_tx))
## There are 6 results, one for each transcript.
mapFromTranscripts(coord, Dgkb_cds_by_tx)
## To map exon coordinates to genome space the 'transcripts'
## argument is the collection of exon regions by transcript.
coord <- GRanges("chr12", IRanges(100, width=1))
ex_by_tx <- exonsBy(txdb, "tx", use.names=TRUE)
Dgkb_ex_by_tx <- ex_by_tx[Dgkb_tx_names]
names(Dgkb_ex_by_tx) <- rep(seqnames(coord), length(Dgkb_ex_by_tx))
## Again the output has 6 results, one for each transcript.
mapFromTranscripts(coord, Dgkb_ex_by_tx)
## To go the reverse direction and map from genome space to
## local CDS or exon space, use mapToTranscripts().
## Genomic position 37981944 maps to CDS position 4:
coord <- GRanges("chr12", IRanges(37981944, width=1))
mapToTranscripts(coord, Dgkb_cds_by_tx)
## Genomic position 37880273 maps to exon position 100:
coord <- GRanges("chr12", IRanges(37880273, width=1))
mapToTranscripts(coord, Dgkb_ex_by_tx)
## The following examples use more than 2GB of memory, which is more
## than what 32-bit Windows can handle:
is_32bit_windows <- .Platform$OS.type == "windows" &&
.Platform$r_arch == "i386"
if (!is_32bit_windows) {
## ---------------------------------------------------------------------
## B. Map sequence locations in exons to the genome
## ---------------------------------------------------------------------
## NAGNAG alternative splicing plays an essential role in biological
## processes and represents a highly adaptable system for
## posttranslational regulation of gene function. The majority of
## NAGNAG studies largely focus on messenger RNA. A study by Sun,
## Lin, and Yan (http://www.hindawi.com/journals/bmri/2014/736798/)
## demonstrated that NAGNAG splicing is also operative in large
## intergenic noncoding RNA (lincRNA). One finding of interest was
## that linc-POLR3G-10 exhibited two NAGNAG acceptors located in two
## distinct transcripts: TCONS_00010012 and TCONS_00010010.
## Extract the exon coordinates of TCONS_00010012 and TCONS_00010010:
lincrna <- c("TCONS_00010012", "TCONS_00010010")
library(TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts)
txdb <- TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts
exons <- exonsBy(txdb, by="tx", use.names=TRUE)[lincrna]
exons
## The two NAGNAG acceptors were identified in the upstream region of
## the fourth and fifth exons located in TCONS_00010012.
## Extract the sequences for transcript TCONS_00010012:
library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
exons_seq <- getSeq(genome, exons[[1]])
## TCONS_00010012 has 4 exons:
exons_seq
## The most common triplet among the lincRNA sequences was CAG. Identify
## the location of this pattern in all exons.
cag_loc <- vmatchPattern("CAG", exons_seq)
## Convert the first occurance of CAG in each exon back to genome
## coordinates.
first_loc <- do.call(c, sapply(cag_loc, "[", 1, simplify=TRUE))
pmapFromTranscripts(first_loc, exons[[1]])
## ---------------------------------------------------------------------
## C. Map dbSNP variants to CDS or cDNA coordinates
## ---------------------------------------------------------------------
## The GIPR gene encodes a G-protein coupled receptor for gastric
## inhibitory polypeptide (GIP). Originally GIP was identified to
## inhibited gastric acid secretion and gastrin release but was later
## demonstrated to stimulate insulin release in the presence of elevated
## glucose.
## In this example 5 SNPs located in the GIPR gene are mapped to cDNA
## coordinates. A list of SNPs in GIPR can be downloaded from dbSNP or
## NCBI.
rsids <- c("rs4803846", "rs139322374", "rs7250736", "rs7250754",
"rs9749185")
## Extract genomic coordinates with a SNPlocs package.
library(SNPlocs.Hsapiens.dbSNP144.GRCh38)
snps <- snpsById(SNPlocs.Hsapiens.dbSNP144.GRCh38, rsids)
## Gene regions of GIPR can be extracted from a TxDb package of
## compatible build. The TxDb package uses Entrez gene identifiers
## and GIPR is a gene symbol. Let's first lookup its Entrez gene ID.
library(org.Hs.eg.db)
GIPR_geneid <- get("GIPR", org.Hs.egSYMBOL2EG)
## The transcriptsBy() extractor returns a range for each transcript that
## includes the UTR and exon regions (i.e., cDNA).
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
tx_by_gene <- transcriptsBy(txdb, "gene")
GIPR_transcripts <- tx_by_gene[GIPR_geneid]
GIPR_transcripts # all 8 GIPR transcripts are on chr19, positive strand
## Before mapping, the chromosome names (seqlevels) in the two
## objects must be harmonized. The style is NCBI for 'snps' and
## UCSC for 'GIPR_transcripts'.
seqlevelsStyle(snps)
seqlevelsStyle(GIPR_transcripts)
## Modify the style (and genome) in 'snps' to match 'GIPR_transcripts'.
seqlevelsStyle(snps) <- seqlevelsStyle(GIPR_transcripts)
## The 'GIPR_transcripts' object is a GRangesList of length 1. This single
## list element contains the cDNA range for 8 different transcripts. To
## map to each transcript individually 'GIPR_transcripts' must be unlisted
## before mapping.
## Map all 5 SNPS to all 8 transcripts:
mapToTranscripts(snps, unlist(GIPR_transcripts))
## Map the first SNP to transcript "ENST00000590918.5" and the second to
## "ENST00000263281.7".
pmapToTranscripts(snps[1:2], unlist(GIPR_transcripts)[1:2])
## The cdsBy() extractor returns coding regions by gene or by transcript.
## Extract the coding regions for transcript "ENST00000263281.7".
cds <- cdsBy(txdb, "tx", use.names=TRUE)["ENST00000263281.7"]
cds
## The 'cds' object is a GRangesList of length 1 containing all CDS ranges
## for the single transcript "ENST00000263281.7".
## To map to the concatenated group of ranges leave 'cds' as a GRangesList.
mapToTranscripts(snps, cds)
## Only the second SNP could be mapped. Unlisting the 'cds' object maps
## the SNPs to the individual cds ranges (vs the concatenated range).
mapToTranscripts(snps[2], unlist(cds))
## The location is the same because the SNP hit the first CDS range. If
## the transcript were on the "-" strand the difference in concatenated
## vs non-concatenated position would be more obvious.
## Change strand:
strand(cds) <- strand(snps) <- "-"
mapToTranscripts(snps[2], unlist(cds))
}
jmacdon/GenomicFeatures documentation built on Jan. 2, 2022, 7:40 a.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) See Also Examples
Description
Map range coordinates between features in the transcriptome and genome (reference) space.
See ?mapToAlignments in the
GenomicAlignments package for mapping coordinates between
reads (local) and genome (reference) space using a CIGAR alignment.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | ## mapping to transcripts
## S4 method for signature 'GenomicRanges,GenomicRanges'
mapToTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'GenomicRanges,GRangesList'
mapToTranscripts(x, transcripts,
ignore.strand = FALSE, intronJunctions=FALSE)
## S4 method for signature 'ANY,TxDb'
mapToTranscripts(x, transcripts, ignore.strand = FALSE,
extractor.fun = GenomicFeatures::transcripts, ...)
## S4 method for signature 'GenomicRanges,GRangesList'
pmapToTranscripts(x, transcripts,
ignore.strand = FALSE)
## mapping from transcripts
## S4 method for signature 'GenomicRanges,GRangesList'
mapFromTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'GenomicRanges,GRangesList'
pmapFromTranscripts(x, transcripts,
ignore.strand = FALSE)
## S4 method for signature 'IntegerRanges,GRangesList'
pmapFromTranscripts(x, transcripts)
|
Arguments
x |
GenomicRanges object of positions to be mapped.
The seqnames of |
transcripts |
A named GenomicRanges or
GRangesList object used to map between The |
ignore.strand |
When When Mapped position is computed by counting from the transcription start site
(TSS) and is not affected by the value of |
intronJunctions |
Logical to indicate if intronic ranges in This argument is only supported in When Ranges that have either the start or end in an intron are considered "non
hits" and are never mapped. Ranges that span introns are always mapped.
Neither of these range types are controlled by the |
extractor.fun |
Function to extract genomic features from a This argument is only applicable to Valid
|
... |
Additional arguments passed to |
Details
In GenomicFeatures >= 1.21.10, the default for ignore.strand was
changed to FALSE for consistency with other methods in the
GenomicRanges and GenomicAlignments packages. Additionally,
the mapped position is computed from the TSS and does not depend on the
ignore.strand argument.
See the section on ignore.strand for details.
mapToTranscripts,pmapToTranscriptsThe genomic range inxis mapped to the local position in thetranscriptsranges. A successful mapping occurs whenxis completely within thetranscriptsrange, equivalent to:findOverlaps(..., type="within")Transcriptome-based coordinates start counting at 1 at the beginning of the
transcriptsrange and return positions wherexwas aligned. The seqlevels of the return object are taken from thetranscriptsobject and should be transcript names. In this direction, mapping is attempted between all elements ofxand all elements oftranscripts.mapToTranscriptsusesfindOverlapsto map ranges inxto ranges intranscripts. This method does not return unmapped ranges.pmapToTranscriptsmaps the i-th range inxto the i-th range intranscripts. Recycling is supported for bothxandtranscriptswhen either is length == 1L; otherwise the lengths must match. Ranges inxthat do not map (out of bounds or strand mismatch) are returned as zero-width ranges starting at 0. These ranges are given the seqname of "UNMAPPED".mapFromTranscripts,pmapFromTranscriptsThe transcript-based position inxis mapped to genomic coordinates using the ranges intranscripts. A successful mapping occurs when the following is TRUE:width(transcripts) >= start(x) + width(x)xis aligned totranscriptsby moving instart(x)positions in from the beginning of thetranscriptsrange. The seqlevels of the return object are chromosome names.mapFromTranscriptsuses the seqname ofxand the names oftranscriptsto determine mapping pairs (vs attempting to match all possible pairs). Name matching is motivated by use cases such as differentially expressed regions where the expressed regions inxwould only be related to a subset of regions intranscripts. This method does not return unmapped ranges.pmapFromTranscriptsmaps the i-th range inxto the i-th range intranscriptsand therefore does not use name matching. Recycling is supported inpmapFromTranscriptswhen eitherxortranscriptsis length == 1L; otherwise the lengths must match. Ranges inxthat do not map (out of bounds or strand mismatch) are returned as zero-width ranges starting at 0. These ranges are given the seqname of "UNMAPPED".
Value
pmapToTranscripts returns a GRanges the same length as
x.
pmapFromTranscripts returns a GRanges when transcripts
is a GRanges and a GRangesList when transcripts
is a GRangesList. In both cases the return object is the same
length as x. The rational for returning the GRangesList is
to preserve exon structure; ranges in a list element that are not overlapped
by x are returned as a zero-width range. The GRangesList
return object will have no seqlevels called "UNMAPPED"; those will only
occur when a GRanges is returned.
mapToTranscripts and mapFromTranscripts return GRanges
objects that vary in length similar to a Hits object. The result
contains mapped records only; strand mismatch and out of bound ranges are
not returned. xHits and transcriptsHits metadata columns
(similar to the queryHits and subjectHits of a Hits
object) indicate elements of x and transcripts used in
the mapping.
When intronJunctions is TRUE, mapToTranscripts returns an
extra metdata column named intronic to identify the intron ranges.
When mapping to transcript coordinates, seqlevels of the output are the names
on the transcripts object and most often these will be transcript
names. When mapping to the genome, seqlevels of the output are the seqlevels
of transcripts which are usually chromosome names.
Author(s)
V. Obenchain, M. Lawrence and H. Pagès
See Also
-
?mapToAlignmentsin the GenomicAlignments package for methods mapping between reads and genome space using a CIGAR alignment.
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 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 | ## ---------------------------------------------------------------------
## A. Basic Use: Conversion between CDS and Exon coordinates and the
## genome
## ---------------------------------------------------------------------
## Gene "Dgkb" has ENTREZID "217480":
library(org.Mm.eg.db)
Dgkb_geneid <- get("Dgkb", org.Mm.egSYMBOL2EG)
## The gene is on the positive strand, chromosome 12:
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
tx_by_gene <- transcriptsBy(txdb, by="gene")
Dgkb_transcripts <- tx_by_gene[[Dgkb_geneid]]
Dgkb_transcripts # all 7 Dgkb transcripts are on chr12, positive strand
## To map coordinates from local CDS or exon space to genome
## space use mapFromTranscripts().
## When mapping CDS coordinates to genome space the 'transcripts'
## argument is the collection of CDS regions by transcript.
coord <- GRanges("chr12", IRanges(4, width=1))
## Get the names of the transcripts in the gene:
Dgkb_tx_names <- mcols(Dgkb_transcripts)$tx_name
Dgkb_tx_names
## Use these names to isolate the region of interest:
cds_by_tx <- cdsBy(txdb, "tx", use.names=TRUE)
Dgkb_cds_by_tx <- cds_by_tx[intersect(Dgkb_tx_names, names(cds_by_tx))]
## Dgkb CDS grouped by transcript (no-CDS transcripts omitted):
Dgkb_cds_by_tx
lengths(Dgkb_cds_by_tx) # nb of CDS per transcript
## A requirement for mapping from transcript space to genome space
## is that seqnames in 'x' match the names in 'transcripts'.
names(Dgkb_cds_by_tx) <- rep(seqnames(coord), length(Dgkb_cds_by_tx))
## There are 6 results, one for each transcript.
mapFromTranscripts(coord, Dgkb_cds_by_tx)
## To map exon coordinates to genome space the 'transcripts'
## argument is the collection of exon regions by transcript.
coord <- GRanges("chr12", IRanges(100, width=1))
ex_by_tx <- exonsBy(txdb, "tx", use.names=TRUE)
Dgkb_ex_by_tx <- ex_by_tx[Dgkb_tx_names]
names(Dgkb_ex_by_tx) <- rep(seqnames(coord), length(Dgkb_ex_by_tx))
## Again the output has 6 results, one for each transcript.
mapFromTranscripts(coord, Dgkb_ex_by_tx)
## To go the reverse direction and map from genome space to
## local CDS or exon space, use mapToTranscripts().
## Genomic position 37981944 maps to CDS position 4:
coord <- GRanges("chr12", IRanges(37981944, width=1))
mapToTranscripts(coord, Dgkb_cds_by_tx)
## Genomic position 37880273 maps to exon position 100:
coord <- GRanges("chr12", IRanges(37880273, width=1))
mapToTranscripts(coord, Dgkb_ex_by_tx)
## The following examples use more than 2GB of memory, which is more
## than what 32-bit Windows can handle:
is_32bit_windows <- .Platform$OS.type == "windows" &&
.Platform$r_arch == "i386"
if (!is_32bit_windows) {
## ---------------------------------------------------------------------
## B. Map sequence locations in exons to the genome
## ---------------------------------------------------------------------
## NAGNAG alternative splicing plays an essential role in biological
## processes and represents a highly adaptable system for
## posttranslational regulation of gene function. The majority of
## NAGNAG studies largely focus on messenger RNA. A study by Sun,
## Lin, and Yan (http://www.hindawi.com/journals/bmri/2014/736798/)
## demonstrated that NAGNAG splicing is also operative in large
## intergenic noncoding RNA (lincRNA). One finding of interest was
## that linc-POLR3G-10 exhibited two NAGNAG acceptors located in two
## distinct transcripts: TCONS_00010012 and TCONS_00010010.
## Extract the exon coordinates of TCONS_00010012 and TCONS_00010010:
lincrna <- c("TCONS_00010012", "TCONS_00010010")
library(TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts)
txdb <- TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts
exons <- exonsBy(txdb, by="tx", use.names=TRUE)[lincrna]
exons
## The two NAGNAG acceptors were identified in the upstream region of
## the fourth and fifth exons located in TCONS_00010012.
## Extract the sequences for transcript TCONS_00010012:
library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
exons_seq <- getSeq(genome, exons[[1]])
## TCONS_00010012 has 4 exons:
exons_seq
## The most common triplet among the lincRNA sequences was CAG. Identify
## the location of this pattern in all exons.
cag_loc <- vmatchPattern("CAG", exons_seq)
## Convert the first occurance of CAG in each exon back to genome
## coordinates.
first_loc <- do.call(c, sapply(cag_loc, "[", 1, simplify=TRUE))
pmapFromTranscripts(first_loc, exons[[1]])
## ---------------------------------------------------------------------
## C. Map dbSNP variants to CDS or cDNA coordinates
## ---------------------------------------------------------------------
## The GIPR gene encodes a G-protein coupled receptor for gastric
## inhibitory polypeptide (GIP). Originally GIP was identified to
## inhibited gastric acid secretion and gastrin release but was later
## demonstrated to stimulate insulin release in the presence of elevated
## glucose.
## In this example 5 SNPs located in the GIPR gene are mapped to cDNA
## coordinates. A list of SNPs in GIPR can be downloaded from dbSNP or
## NCBI.
rsids <- c("rs4803846", "rs139322374", "rs7250736", "rs7250754",
"rs9749185")
## Extract genomic coordinates with a SNPlocs package.
library(SNPlocs.Hsapiens.dbSNP144.GRCh38)
snps <- snpsById(SNPlocs.Hsapiens.dbSNP144.GRCh38, rsids)
## Gene regions of GIPR can be extracted from a TxDb package of
## compatible build. The TxDb package uses Entrez gene identifiers
## and GIPR is a gene symbol. Let's first lookup its Entrez gene ID.
library(org.Hs.eg.db)
GIPR_geneid <- get("GIPR", org.Hs.egSYMBOL2EG)
## The transcriptsBy() extractor returns a range for each transcript that
## includes the UTR and exon regions (i.e., cDNA).
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
tx_by_gene <- transcriptsBy(txdb, "gene")
GIPR_transcripts <- tx_by_gene[GIPR_geneid]
GIPR_transcripts # all 8 GIPR transcripts are on chr19, positive strand
## Before mapping, the chromosome names (seqlevels) in the two
## objects must be harmonized. The style is NCBI for 'snps' and
## UCSC for 'GIPR_transcripts'.
seqlevelsStyle(snps)
seqlevelsStyle(GIPR_transcripts)
## Modify the style (and genome) in 'snps' to match 'GIPR_transcripts'.
seqlevelsStyle(snps) <- seqlevelsStyle(GIPR_transcripts)
## The 'GIPR_transcripts' object is a GRangesList of length 1. This single
## list element contains the cDNA range for 8 different transcripts. To
## map to each transcript individually 'GIPR_transcripts' must be unlisted
## before mapping.
## Map all 5 SNPS to all 8 transcripts:
mapToTranscripts(snps, unlist(GIPR_transcripts))
## Map the first SNP to transcript "ENST00000590918.5" and the second to
## "ENST00000263281.7".
pmapToTranscripts(snps[1:2], unlist(GIPR_transcripts)[1:2])
## The cdsBy() extractor returns coding regions by gene or by transcript.
## Extract the coding regions for transcript "ENST00000263281.7".
cds <- cdsBy(txdb, "tx", use.names=TRUE)["ENST00000263281.7"]
cds
## The 'cds' object is a GRangesList of length 1 containing all CDS ranges
## for the single transcript "ENST00000263281.7".
## To map to the concatenated group of ranges leave 'cds' as a GRangesList.
mapToTranscripts(snps, cds)
## Only the second SNP could be mapped. Unlisting the 'cds' object maps
## the SNPs to the individual cds ranges (vs the concatenated range).
mapToTranscripts(snps[2], unlist(cds))
## The location is the same because the SNP hit the first CDS range. If
## the transcript were on the "-" strand the difference in concatenated
## vs non-concatenated position would be more obvious.
## Change strand:
strand(cds) <- strand(snps) <- "-"
mapToTranscripts(snps[2], unlist(cds))
}
|
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.