Shaun Jackman — Feb 27, 2013, 11:59 AM
# STAT 540 Seminar 7
# Shaun Jackman
library(biomaRt)
library(easyRNASeq)
Loading required package: parallel```
Loading required package: genomeIntervals“
Loading required package: intervals```
Loading required package: BiocGenerics”
Attaching package: 'BiocGenerics'```
The following object(s) are masked from 'package:stats':
xtabs“
The following object(s) are masked from 'package:base':
anyDuplicated, cbind, colnames, duplicated, eval, Filter, Find, get,
intersect, lapply, Map, mapply, mget, order, paste, pmax, pmax.int, pmin,
pmin.int, Position, rbind, Reduce, rep.int, rownames, sapply, setdiff,
table, tapply, union, unique```
Warning: replacing previous import 'coerce' when loading 'intervals'”
Warning: replacing previous import 'initialize' when loading 'intervals'```
Loading required package: Biobase“
Welcome to Bioconductor
Vignettes contain introductory material; view with 'browseVignettes()'. To
cite Bioconductor, see 'citation("Biobase")', and for packages
'citation("pkgname")'.```
Loading required package: edgeR”
Loading required package: limma```
Loading required package: Biostrings“
Loading required package: IRanges```
Attaching package: 'IRanges'”
The following object(s) are masked from 'package:intervals':
reduce```
Attaching package: 'Biostrings'“
The following object(s) are masked from 'package:intervals':
type```
Loading required package: BSgenome”
Loading required package: GenomicRanges```
Loading required package: DESeq“
Loading required package: locfit```
locfit 1.5-8 2012-04-25”
Attaching package: 'locfit'```
The following object(s) are masked from 'package:GenomicRanges':
left, right“
Loading required package: Rsamtools```
Loading required package: ShortRead”
Loading required package: lattice```
Loading required package: latticeExtra“
Loading required package: RColorBrewer```
```r
library(Rsamtools)
library(ShortRead)
library(BSgenome.Dmelanogaster.UCSC.dm3)
# Reading BAM Files
bamDat <- readAligned(
"../data/drosophilaRnaSeq/drosophilaMelanogasterSubset.bam",
type = "BAM")
str(bamDat)
Formal class 'AlignedRead' [package "ShortRead"] with 8 slots
..@ chromosome : Factor w/ 15 levels "chrYHet","chrM",..: 2 2 2 2 2 2 2 2 2 2 ...
..@ position : int [1:64206] 548 1497 1506 1528 1540 1552 1552 1555 1559 1566 ...
..@ strand : Factor w/ 3 levels "+","-","*": 2 1 1 1 1 1 1 1 2 2 ...
..@ alignQuality:Formal class 'NumericQuality' [package "ShortRead"] with 1 slots
.. .. ..@ quality: int [1:64206] 132 132 127 130 130 122 132 132 132 132 ...
..@ alignData :Formal class 'AlignedDataFrame' [package "ShortRead"] with 4 slots
.. .. ..@ varMetadata :'data.frame': 1 obs. of 1 variable:
.. .. .. ..$ labelDescription: chr "Type of read; see ?scanBam"
.. .. ..@ data :'data.frame': 64206 obs. of 1 variable:
.. .. .. ..$ flag: int [1:64206] 16 0 0 0 0 0 0 0 16 16 ...
.. .. ..@ dimLabels : chr [1:2] "readName" "alignColumn"
.. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slots
.. .. .. .. ..@ .Data:List of 1
.. .. .. .. .. ..$ : int [1:3] 1 1 0
..@ quality :Formal class 'FastqQuality' [package "ShortRead"] with 1 slots
.. .. ..@ quality:Formal class 'BStringSet' [package "Biostrings"] with 5 slots
.. .. .. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. .. .. ..$ :<externalptr>
.. .. .. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. .. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. .. .. ..@ group : int [1:64206] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. .. .. ..@ start : int [1:64206] 1 37 73 109 145 181 217 253 289 325 ...
.. .. .. .. .. .. ..@ width : int [1:64206] 36 36 36 36 36 36 36 36 36 36 ...
.. .. .. .. .. .. ..@ NAMES : NULL
.. .. .. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. .. .. ..@ metadata : list()
.. .. .. .. ..@ elementType : chr "BString"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
..@ sread :Formal class 'DNAStringSet' [package "Biostrings"] with 5 slots
.. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. ..$ :<externalptr>
.. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. ..@ group : int [1:64206] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. ..@ start : int [1:64206] 1 37 73 109 145 181 217 253 289 325 ...
.. .. .. .. ..@ width : int [1:64206] 36 36 36 36 36 36 36 36 36 36 ...
.. .. .. .. ..@ NAMES : NULL
.. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
.. .. ..@ elementType : chr "DNAString"
.. .. ..@ elementMetadata: NULL
.. .. ..@ metadata : list()
..@ id :Formal class 'BStringSet' [package "Biostrings"] with 5 slots
.. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. ..$ :<externalptr>
.. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. ..@ group : int [1:64206] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. ..@ start : int [1:64206] 1 29 57 85 114 144 173 202 230 260 ...
.. .. .. .. ..@ width : int [1:64206] 28 28 28 29 30 29 29 28 30 30 ...
.. .. .. .. ..@ NAMES : NULL
.. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
.. .. ..@ elementType : chr "BString"
.. .. ..@ elementMetadata: NULL
.. .. ..@ metadata : list()
indexFile <- indexBam("../data/drosophilaRnaSeq/drosophilaMelanogasterSubset.bam")
# Filtering BAM Files
nFilt <- nFilter(2)
chrFilt <- chromosomeFilter(regex = "chr")
filt <- compose(nFilt, chrFilt)
bamDatFiltered <- bamDat[filt(bamDat)]
# Examining BAM Data
str(bamDatFiltered)
Formal class 'AlignedRead' [package "ShortRead"] with 8 slots
..@ chromosome : Factor w/ 7 levels "chrM","chr2L",..: 1 1 1 1 1 1 1 1 1 1 ...
..@ position : int [1:56883] 548 1497 1506 1528 1540 1552 1552 1555 1559 1566 ...
..@ strand : Factor w/ 3 levels "+","-","*": 2 1 1 1 1 1 1 1 2 2 ...
..@ alignQuality:Formal class 'NumericQuality' [package "ShortRead"] with 1 slots
.. .. ..@ quality: int [1:56883] 132 132 127 130 130 122 132 132 132 132 ...
..@ alignData :Formal class 'AlignedDataFrame' [package "ShortRead"] with 4 slots
.. .. ..@ varMetadata :'data.frame': 1 obs. of 1 variable:
.. .. .. ..$ labelDescription: chr "Type of read; see ?scanBam"
.. .. ..@ data :'data.frame': 56883 obs. of 1 variable:
.. .. .. ..$ flag: int [1:56883] 16 0 0 0 0 0 0 0 16 16 ...
.. .. ..@ dimLabels : chr [1:2] "readName" "alignColumn"
.. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slots
.. .. .. .. ..@ .Data:List of 1
.. .. .. .. .. ..$ : int [1:3] 1 1 0
..@ quality :Formal class 'FastqQuality' [package "ShortRead"] with 1 slots
.. .. ..@ quality:Formal class 'BStringSet' [package "Biostrings"] with 5 slots
.. .. .. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. .. .. ..$ :<externalptr>
.. .. .. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. .. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. .. .. ..@ group : int [1:56883] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. .. .. ..@ start : int [1:56883] 1 37 73 109 145 181 217 253 289 325 ...
.. .. .. .. .. .. ..@ width : int [1:56883] 36 36 36 36 36 36 36 36 36 36 ...
.. .. .. .. .. .. ..@ NAMES : NULL
.. .. .. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. .. .. ..@ metadata : list()
.. .. .. .. ..@ elementType : chr "BString"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
..@ sread :Formal class 'DNAStringSet' [package "Biostrings"] with 5 slots
.. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. ..$ :<externalptr>
.. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. ..@ group : int [1:56883] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. ..@ start : int [1:56883] 1 37 73 109 145 181 217 253 289 325 ...
.. .. .. .. ..@ width : int [1:56883] 36 36 36 36 36 36 36 36 36 36 ...
.. .. .. .. ..@ NAMES : NULL
.. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
.. .. ..@ elementType : chr "DNAString"
.. .. ..@ elementMetadata: NULL
.. .. ..@ metadata : list()
..@ id :Formal class 'BStringSet' [package "Biostrings"] with 5 slots
.. .. ..@ pool :Formal class 'SharedRaw_Pool' [package "IRanges"] with 2 slots
.. .. .. .. ..@ xp_list :List of 1
.. .. .. .. .. ..$ :<externalptr>
.. .. .. .. ..@ .link_to_cached_object_list:List of 1
.. .. .. .. .. ..$ :<environment: 0x10649faa8>
.. .. ..@ ranges :Formal class 'GroupedIRanges' [package "IRanges"] with 7 slots
.. .. .. .. ..@ group : int [1:56883] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. ..@ start : int [1:56883] 1 29 57 85 114 144 173 202 230 260 ...
.. .. .. .. ..@ width : int [1:56883] 28 28 28 29 30 29 29 28 30 30 ...
.. .. .. .. ..@ NAMES : NULL
.. .. .. .. ..@ elementType : chr "integer"
.. .. .. .. ..@ elementMetadata: NULL
.. .. .. .. ..@ metadata : list()
.. .. ..@ elementType : chr "BString"
.. .. ..@ elementMetadata: NULL
.. .. ..@ metadata : list()
levels(chromosome(bamDatFiltered))
[1] "chrM" "chr2L" "chrX" "chr3L" "chr4" "chr2R" "chr3R"
id(bamDatFiltered)[1:10]
A BStringSet instance of length 10
width seq
[1] 28 HWI-EAS225_90320:3:1:141:680
[2] 28 HWI-EAS225_90320:3:1:660:332
[3] 28 HWI-EAS225_90320:3:1:164:226
[4] 29 HWI-EAS225_90320:3:1:1088:176
[5] 30 HWI-EAS225_90320:3:1:1038:1484
[6] 29 HWI-EAS225_90320:3:1:850:1742
[7] 29 HWI-EAS225_90320:3:1:1319:586
[8] 28 HWI-EAS225_90320:3:1:103:631
[9] 30 HWI-EAS225_90320:3:1:1353:1498
[10] 30 HWI-EAS225_90320:3:1:1092:1016
sread(bamDatFiltered)[1:10]
A DNAStringSet instance of length 10
width seq
[1] 36 GGAAATCAAAAATGGAAAGGAGCGGCTCCACTTTTT
[2] 36 AAATCATAAAGATATTGGAACTTTATATTTTATTTT
[3] 36 AGATATTGGAACTTTATATTTTATTTTTGGAGCTTG
[4] 36 ATTTTTGGAGCTTGAGCTGGAATAGTTGGAACATCT
[5] 36 TGAGCTGGAATAGTTGGAACATCTTTAAGAATTTTA
[6] 36 GTTGGAACATCTTTAAGAATTTTAATTAGAGCTGAA
[7] 36 GTTGGAACATCTTTAAGAATTTTAATTCGAGCTGAA
[8] 36 GGAACATCTTTAAGAATTTTAATTCGAGCTGAATTA
[9] 36 GTCCTAATTCAGCTCGAATTAAAATTCTTAAAGATG
[10] 36 CCAGGATGTCCTAATTCAGCTCGAATTAAAATTCTT
quality(bamDatFiltered)[1:10]
class: FastqQuality
quality:
A BStringSet instance of length 10
width seq
[1] 36 BACBCCABBBBBA>8@B@@==>;5-9A<;=7A@@B@
[2] 36 BCBBABBA@@B;B>AB@@<>:AAA9?>?A@A<?A@@
[3] 36 @?8AB>A?=)A=@*8>6/@3>A)/@4>?BA'(-1B=
[4] 36 BBCACCA@-4ABC62?*;A?BBA?B@.8B9?33;+=
[5] 36 ?5@A4::@@55;;89<'6?A8@A=4@=>54>76);A
[6] 36 A8=B;462>;7BCBAA>1;=</?BA94%<:?(7@9=
[7] 36 BBB>ABB@@BBBBCBCC@7ABBBAABB@B?AAA@=@
[8] 36 =B=ACCBBC8ACCCBBBCCCBB=CAB9=BBBB@2?:
[9] 36 @7ABBBBABBB?BAAB=@CBB;7ABAABBA?@;2=A
[10] 36 BB>B?:ABABBBBABCB@@@BB@:@BA;>@;@B?AB
position(bamDatFiltered)[1:10]
[1] 548 1497 1506 1528 1540 1552 1552 1555 1559 1566
strand(bamDatFiltered)[1:10]
[1] - + + + + + + + - -
Levels: + - *
# What are the differences between the filtered and unfiltered BAM files?
length(position(bamDat))
[1] 64206
length(position(bamDatFiltered))
[1] 56883
# The unfiltered data has 64206 reads. The filtered data has 56883 reads.
# What are the chromosomes with aligned reads from the BAM file?
summary(chromosome(bamDat))
chrYHet chrM chr2L chrX chr3L chr4 chr2R
0 1921 9839 8892 9231 739 12370
chr3R chrUextra chr2RHet chr2LHet chr3LHet chr3RHet chrU
13891 0 0 0 0 0 0
chrXHet NA's
0 7323
summary(chromosome(bamDatFiltered))
chrM chr2L chrX chr3L chr4 chr2R chr3R
1921 9839 8892 9231 739 12370 13891
# Accessing Genome Annotations
chrSizes <- seqlengths(Dmelanogaster)
ensembl <- useMart("ensembl", dataset = "dmelanogaster_gene_ensembl")
annotation.fields <- c("ensembl_gene_id", "strand",
"chromosome_name", "start_position", "end_position")
gene.annotation <- getBM(annotation.fields, mart = ensembl,
filters = "chromosome_name", values = c("2L"))
str(gene.annotation)
'data.frame': 2856 obs. of 5 variables:
$ ensembl_gene_id: chr "FBgn0031208" "FBgn0002121" "FBgn0031209" "FBgn0051973" ...
$ strand : int 1 -1 -1 -1 1 1 1 -1 1 1 ...
$ chromosome_name: chr "2L" "2L" "2L" "2L" ...
$ start_position : int 7529 9839 21919 25402 66721 71757 76326 82456 94752 102382 ...
$ end_position : int 9484 21372 25151 59242 71390 76211 77785 87382 102086 104142 ...
levels(as.factor(gene.annotation$chromosome))
[1] "2L"
gene.annotation$chromosome <- paste("chr", gene.annotation$chromosome_name,
sep = "")
levels(as.factor(gene.annotation$chromosome))
[1] "chr2L"
gene.range <- RangedData(
IRanges(
start = gene.annotation$start_position,
end = gene.annotation$end_position),
space = gene.annotation$chromosome, strand = gene.annotation$strand,
gene = gene.annotation$ensembl_gene_id,
universe = "Dm3")
show(gene.range)
RangedData with 2856 rows and 2 value columns across 1 space
space ranges | strand gene
<factor> <IRanges> | <integer> <character>
1 chr2L [ 7529, 9484] | 1 FBgn0031208
2 chr2L [ 9839, 21372] | -1 FBgn0002121
3 chr2L [21919, 25151] | -1 FBgn0031209
4 chr2L [25402, 59242] | -1 FBgn0051973
5 chr2L [66721, 71390] | 1 FBgn0067779
6 chr2L [71757, 76211] | 1 FBgn0031213
7 chr2L [76326, 77785] | 1 FBgn0031214
8 chr2L [82456, 87382] | -1 FBgn0002931
9 chr2L [94752, 102086] | 1 FBgn0031216
... ... ... ... ... ...
2848 chr2L [18471434, 18471533] | 1 FBgn0262463
2849 chr2L [18472034, 18472111] | 1 FBgn0262406
2850 chr2L [18472315, 18472424] | 1 FBgn0262456
2851 chr2L [19569490, 19569561] | -1 FBgn0262460
2852 chr2L [19569897, 19569972] | -1 FBgn0262377
2853 chr2L [19570182, 19570264] | -1 FBgn0262449
2854 chr2L [20487441, 20487531] | 1 FBgn0262455
2855 chr2L [20616617, 20616714] | 1 FBgn0262445
2856 chr2L [20618366, 20618462] | -1 FBgn0262425
# Calculating Coverage
cover <- coverage(bamDatFiltered, width = chrSizes)
gene.coverage <- aggregate(cover[match(names(gene.range), names(cover))],
ranges(gene.range), sum)
gene.coverage <- ceiling(gene.coverage/unique(width(bamDat)))
length(gene.coverage[["chr2L"]])
[1] 2856
length(ranges(gene.range)$chr2L)
[1] 2856
countTable <- data.frame(chromosome = gene.range$space, gene_start = start(gene.range$ranges),
gene_end = end(gene.range$ranges), strand = gene.range$strand, gene = gene.range$gene,
count = as.vector(gene.coverage[["chr2L"]]))
dim(countTable)
[1] 2856 6
head(countTable)
chromosome gene_start gene_end strand gene count
1 chr2L 7529 9484 1 FBgn0031208 1
2 chr2L 9839 21372 -1 FBgn0002121 47
3 chr2L 21919 25151 -1 FBgn0031209 0
4 chr2L 25402 59242 -1 FBgn0051973 1
5 chr2L 66721 71390 1 FBgn0067779 6
6 chr2L 71757 76211 1 FBgn0031213 0
countTable <- data.frame(chromosome = gene.range$space, gene_start = start(gene.range$ranges),
gene_end = end(gene.range$ranges), strand = gene.range$strand, gene = gene.range$gene,
count = as.vector(gene.coverage[["chr2L"]]), RPKM = (as.vector(gene.coverage[["chr2L"]])
/ (end(gene.range$ranges) - start(gene.range$ranges)))
* (1e+09/length(bamDat)))
head(countTable)
chromosome gene_start gene_end strand gene count RPKM
1 chr2L 7529 9484 1 FBgn0031208 1 7.9667
2 chr2L 9839 21372 -1 FBgn0002121 47 63.4717
3 chr2L 21919 25151 -1 FBgn0031209 0 0.0000
4 chr2L 25402 59242 -1 FBgn0051973 1 0.4603
5 chr2L 66721 71390 1 FBgn0067779 6 20.0148
6 chr2L 71757 76211 1 FBgn0031213 0 0.0000
# Take Home Problem
# Create a similar count table for all the exons located on chr2L.
exon.annotation <- getBM(
c("ensembl_exon_id", "strand", "chromosome_name",
"exon_chrom_start", "exon_chrom_end"),
mart = ensembl, filters = "chromosome_name", values = c("2L"))
str(exon.annotation)
'data.frame': 13210 obs. of 5 variables:
$ ensembl_exon_id : chr "FBgn0031208:1" "FBgn0031208:3" "FBgn0031208:2" "FBgn0031208:4" ...
$ strand : int 1 1 1 1 -1 -1 -1 -1 -1 -1 ...
$ chromosome_name : chr "2L" "2L" "2L" "2L" ...
$ exon_chrom_start: int 7529 8193 8193 8668 21066 19880 17053 14933 13683 13520 ...
$ exon_chrom_end : int 8116 9484 8589 9484 21372 20015 17212 15711 14874 13625 ...
exon.annotation$chromosome <- paste("chr", exon.annotation$chromosome_name,
sep = "")
exon.range <- RangedData(
IRanges(
start = exon.annotation$exon_chrom_start,
end = exon.annotation$exon_chrom_end),
space = exon.annotation$chromosome, strand = exon.annotation$strand,
exon = exon.annotation$ensembl_exon_id,
universe = "Dm3")
show(exon.range)
RangedData with 13210 rows and 2 value columns across 1 space
space ranges | strand exon
<factor> <IRanges> | <integer> <character>
1 chr2L [ 7529, 8116] | 1 FBgn0031208:1
2 chr2L [ 8193, 9484] | 1 FBgn0031208:3
3 chr2L [ 8193, 8589] | 1 FBgn0031208:2
4 chr2L [ 8668, 9484] | 1 FBgn0031208:4
5 chr2L [21066, 21372] | -1 FBgn0002121:15
6 chr2L [19880, 20015] | -1 FBgn0002121:11
7 chr2L [17053, 17212] | -1 FBgn0002121:8
8 chr2L [14933, 15711] | -1 FBgn0002121:7
9 chr2L [13683, 14874] | -1 FBgn0002121:6
... ... ... ... ... ...
13202 chr2L [19569897, 19569972] | -1 FBtr0304215:1
13203 chr2L [19570189, 19570211] | -1 FBtr0304465:1
13204 chr2L [19570182, 19570264] | -1 FBtr0304464:1
13205 chr2L [20487496, 20487517] | 1 FBtr0304488:1
13206 chr2L [20487441, 20487531] | 1 FBtr0304487:1
13207 chr2L [20616673, 20616694] | 1 FBtr0304455:1
13208 chr2L [20616617, 20616714] | 1 FBtr0304454:1
13209 chr2L [20618375, 20618397] | -1 FBtr0304384:1
13210 chr2L [20618366, 20618462] | -1 FBtr0304383:1
cover <- coverage(bamDatFiltered, width = chrSizes)
exon.coverage <- aggregate(cover[match(names(exon.range), names(cover))],
ranges(exon.range), sum)
exon.coverage <- ceiling(exon.coverage/unique(width(bamDat)))
length(exon.coverage[["chr2L"]])
[1] 13210
length(ranges(exon.range)$chr2L)
[1] 13210
countTable <- data.frame(chromosome = exon.range$space, exon_start = start(exon.range$ranges),
exon_end = end(exon.range$ranges), strand = exon.range$strand, exon = exon.range$exon,
count = as.vector(exon.coverage[["chr2L"]]))
dim(countTable)
[1] 13210 6
head(countTable)
chromosome exon_start exon_end strand exon count
1 chr2L 7529 8116 1 FBgn0031208:1 0
2 chr2L 8193 9484 1 FBgn0031208:3 1
3 chr2L 8193 8589 1 FBgn0031208:2 0
4 chr2L 8668 9484 1 FBgn0031208:4 1
5 chr2L 21066 21372 -1 FBgn0002121:15 0
6 chr2L 19880 20015 -1 FBgn0002121:11 0
countTable <- data.frame(chromosome = exon.range$space, exon_start = start(exon.range$ranges),
exon_end = end(exon.range$ranges), strand = exon.range$strand, exon = exon.range$exon,
count = as.vector(exon.coverage[["chr2L"]]), RPKM = (as.vector(exon.coverage[["chr2L"]])
/ (end(exon.range$ranges) - start(exon.range$ranges)))
* (1e+09/length(bamDat)))
head(countTable)
chromosome exon_start exon_end strand exon count RPKM
1 chr2L 7529 8116 1 FBgn0031208:1 0 0.00
2 chr2L 8193 9484 1 FBgn0031208:3 1 12.06
3 chr2L 8193 8589 1 FBgn0031208:2 0 0.00
4 chr2L 8668 9484 1 FBgn0031208:4 1 19.09
5 chr2L 21066 21372 -1 FBgn0002121:15 0 0.00
6 chr2L 19880 20015 -1 FBgn0002121:11 0 0.00