Libraries
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
## Loading required package: GenomicFeatures
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind,
## colMeans, colnames, colSums, dirname, do.call, duplicated,
## eval, evalq, Filter, Find, get, grep, grepl, intersect,
## is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
## paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
## Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which, which.max,
## which.min
## Loading required package: S4Vectors
## Loading required package: stats4
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
##
## expand.grid
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:grDevices':
##
## windows
## Loading required package: GenomeInfoDb
## Loading required package: GenomicRanges
## Loading required package: AnnotationDbi
## 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")'.
library(GenomicFeatures)
library(AnnotationDbi)
Homo sapiens hg38 transcript database
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
txdb
## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: UCSC
## # Genome: hg38
## # Organism: Homo sapiens
## # Taxonomy ID: 9606
## # UCSC Table: knownGene
## # UCSC Track: GENCODE v24
## # Resource URL: http://genome.ucsc.edu/
## # Type of Gene ID: Entrez Gene ID
## # Full dataset: yes
## # miRBase build ID: NA
## # transcript_nrow: 197782
## # exon_nrow: 581036
## # cds_nrow: 293052
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time: 2016-09-29 13:02:09 +0000 (Thu, 29 Sep 2016)
## # GenomicFeatures version at creation time: 1.25.18
## # RSQLite version at creation time: 1.0.0
## # DBSCHEMAVERSION: 1.1
columns(txdb)
## [1] "CDSCHROM" "CDSEND" "CDSID" "CDSNAME" "CDSSTART"
## [6] "CDSSTRAND" "EXONCHROM" "EXONEND" "EXONID" "EXONNAME"
## [11] "EXONRANK" "EXONSTART" "EXONSTRAND" "GENEID" "TXCHROM"
## [16] "TXEND" "TXID" "TXNAME" "TXSTART" "TXSTRAND"
## [21] "TXTYPE"
keytypes(txdb)
## [1] "CDSID" "CDSNAME" "EXONID" "EXONNAME" "GENEID" "TXID"
## [7] "TXNAME"
genes(txdb) # GRanges object with 24,183 ranges
## GRanges object with 24183 ranges and 1 metadata column:
## seqnames ranges strand | gene_id
## <Rle> <IRanges> <Rle> | <character>
## 1 chr19 58345178-58362751 - | 1
## 10 chr8 18391245-18401218 + | 10
## 100 chr20 44619522-44651742 - | 100
## 1000 chr18 27950966-28177446 - | 1000
## 100009613 chr11 70072434-70075433 - | 100009613
## ... ... ... ... . ...
## 9991 chr9 112217716-112333667 - | 9991
## 9992 chr21 34364024-34371389 + | 9992
## 9993 chr22 19036282-19122454 - | 9993
## 9994 chr6 89829894-89874436 + | 9994
## 9997 chr22 50523568-50526439 - | 9997
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
exons(txdb) # GRanges object with 581,036 ranges
## GRanges object with 581036 ranges and 1 metadata column:
## seqnames ranges strand | exon_id
## <Rle> <IRanges> <Rle> | <integer>
## [1] chr1 29554-30039 + | 1
## [2] chr1 30267-30667 + | 2
## [3] chr1 30366-30503 + | 3
## [4] chr1 30564-30667 + | 4
## [5] chr1 30976-31097 + | 5
## ... ... ... ... . ...
## [581032] chrUn_KI270750v1 148668-148843 + | 581032
## [581033] chrUn_KI270752v1 144-268 + | 581033
## [581034] chrUn_KI270752v1 21813-21944 + | 581034
## [581035] chrUn_KI270752v1 3497-3623 - | 581035
## [581036] chrUn_KI270752v1 9943-10067 - | 581036
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
transcripts(txdb) #GRanges object with 197,782 ranges
## GRanges object with 197782 ranges and 2 metadata columns:
## seqnames ranges strand | tx_id tx_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] chr1 29554-31097 + | 1 uc057aty.1
## [2] chr1 30267-31109 + | 2 uc057atz.1
## [3] chr1 30366-30503 + | 3 uc031tlb.1
## [4] chr1 69091-70008 + | 4 uc001aal.1
## [5] chr1 160446-161525 + | 5 uc057aum.1
## ... ... ... ... . ... ...
## [197778] chrUn_KI270750v1 148668-148843 + | 197778 uc064xrp.1
## [197779] chrUn_KI270752v1 144-268 + | 197779 uc064xrq.1
## [197780] chrUn_KI270752v1 21813-21944 + | 197780 uc064xrt.1
## [197781] chrUn_KI270752v1 3497-3623 - | 197781 uc064xrr.1
## [197782] chrUn_KI270752v1 9943-10067 - | 197782 uc064xrs.1
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
A sliced TxDb dataframe
keys <- keys(txdb, keytype="GENEID")
cols <- c("TXID", "TXSTART")
selected <- select(txdb, keys=keys, columns=cols, keytype="GENEID")
## 'select()' returned 1:many mapping between keys and columns
selected[1:5,]
## GENEID TXID TXSTART
## 1 1 166436 58345178
## 2 1 166437 58346850
## 3 1 166438 58346854
## 4 1 166439 58346858
## 5 1 166440 58346860
Transcripts in CHRNA4 gene
hschrna4_txs <- transcriptsBy(txdb, by = "gene")[["1137"]]
hschrna4_txs # GRanges object with 7 ranges
## GRanges object with 7 ranges and 2 metadata columns:
## seqnames ranges strand | tx_id tx_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] chr20 63343223-63361343 - | 170358 uc002yes.4
## [2] chr20 63343310-63361396 - | 170359 uc010gke.2
## [3] chr20 63344069-63377993 - | 170360 uc061ymb.1
## [4] chr20 63346594-63359699 - | 170362 uc061ymd.1
## [5] chr20 63346594-63361304 - | 170363 uc061yme.1
## [6] chr20 63346594-63361304 - | 170364 uc061ymf.1
## [7] chr20 63350908-63361190 - | 170366 uc061ymh.1
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
Exons in CHRNA4 gene/transcrpt(s)
hschrna4_exbytx <- exonsBy(txdb, "tx", use.names=TRUE)[hschrna4_txs$tx_name]
hschrna4_exbytx
## GRangesList object of length 7:
## $uc002yes.4
## GRanges object with 6 ranges and 3 metadata columns:
## seqnames ranges strand | exon_id exon_name exon_rank
## <Rle> <IRanges> <Rle> | <integer> <character> <integer>
## [1] chr20 63361090-63361343 - | 503108 <NA> 1
## [2] chr20 63359548-63359699 - | 503104 <NA> 2
## [3] chr20 63356371-63356415 - | 503102 <NA> 3
## [4] chr20 63355975-63356084 - | 503100 <NA> 4
## [5] chr20 63349653-63351027 - | 503095 <NA> 5
## [6] chr20 63343223-63346863 - | 503089 <NA> 6
##
## $uc010gke.2
## GRanges object with 7 ranges and 3 metadata columns:
## seqnames ranges strand | exon_id exon_name exon_rank
## [1] chr20 63361090-63361396 - | 503109 <NA> 1
## [2] chr20 63359548-63359699 - | 503104 <NA> 2
## [3] chr20 63356371-63356415 - | 503102 <NA> 3
## [4] chr20 63355975-63356084 - | 503100 <NA> 4
## [5] chr20 63355496-63355595 - | 503097 <NA> 5
## [6] chr20 63349653-63351027 - | 503095 <NA> 6
## [7] chr20 63343310-63346863 - | 503090 <NA> 7
##
## $uc061ymb.1
## GRanges object with 5 ranges and 3 metadata columns:
## seqnames ranges strand | exon_id exon_name exon_rank
## [1] chr20 63377755-63377993 - | 503117 <NA> 1
## [2] chr20 63375253-63375392 - | 503115 <NA> 2
## [3] chr20 63373881-63374532 - | 503113 <NA> 3
## [4] chr20 63349653-63351027 - | 503095 <NA> 4
## [5] chr20 63344069-63346863 - | 503091 <NA> 5
##
## ...
## <4 more elements>
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
elementNROWS(hschrna4_exbytx) #lists number of exons per transcripts/gene
## uc002yes.4 uc010gke.2 uc061ymb.1 uc061ymd.1 uc061yme.1 uc061ymf.1
## 6 7 5 6 6 6
## uc061ymh.1
## 5
barplot(elementNROWS(hschrna4_exbytx))
Introns in CHRNA4 gene/transcript(s)
hschrna4_inbytx <- intronsByTranscript(txdb, use.names=TRUE)[hschrna4_txs$tx_name]
hschrna4_inbytx
## GRangesList object of length 7:
## $uc002yes.4
## GRanges object with 5 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr20 63346864-63349652 -
## [2] chr20 63351028-63355974 -
## [3] chr20 63356085-63356370 -
## [4] chr20 63356416-63359547 -
## [5] chr20 63359700-63361089 -
##
## $uc010gke.2
## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## [1] chr20 63346864-63349652 -
## [2] chr20 63351028-63355495 -
## [3] chr20 63355596-63355974 -
## [4] chr20 63356085-63356370 -
## [5] chr20 63356416-63359547 -
## [6] chr20 63359700-63361089 -
##
## $uc061ymb.1
## GRanges object with 4 ranges and 0 metadata columns:
## seqnames ranges strand
## [1] chr20 63346864-63349652 -
## [2] chr20 63351028-63373880 -
## [3] chr20 63374533-63375252 -
## [4] chr20 63375393-63377754 -
##
## ...
## <4 more elements>
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome
elementNROWS(hschrna4_inbytx)
## uc002yes.4 uc010gke.2 uc061ymb.1 uc061ymd.1 uc061yme.1 uc061ymf.1
## 5 6 4 5 5 5
## uc061ymh.1
## 4
barplot(elementNROWS(hschrna4_exbytx))
Extract exon and intron sequences of transcript/gene(s)
library(BSgenome.Hsapiens.UCSC.hg38)
## Loading required package: BSgenome
## Loading required package: Biostrings
## Loading required package: XVector
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
##
## strsplit
## Loading required package: rtracklayer
BSgenome.Hsapiens.UCSC.hg38
## Human genome:
## # organism: Homo sapiens (Human)
## # provider: UCSC
## # provider version: hg38
## # release date: Dec. 2013
## # release name: Genome Reference Consortium GRCh38
## # 455 sequences:
## # chr1 chr2 chr3
## # chr4 chr5 chr6
## # chr7 chr8 chr9
## # chr10 chr11 chr12
## # chr13 chr14 chr15
## # ... ... ...
## # chrUn_KI270744v1 chrUn_KI270745v1 chrUn_KI270746v1
## # chrUn_KI270747v1 chrUn_KI270748v1 chrUn_KI270749v1
## # chrUn_KI270750v1 chrUn_KI270751v1 chrUn_KI270752v1
## # chrUn_KI270753v1 chrUn_KI270754v1 chrUn_KI270755v1
## # chrUn_KI270756v1 chrUn_KI270757v1
## # (use 'seqnames()' to see all the sequence names, use the '$' or '[['
## # operator to access a given sequence)
Hsapiens <- BSgenome.Hsapiens.UCSC.hg38
hschrna4_ex_seqs <- getSeq(Hsapiens, hschrna4_exbytx) ## extract exon sequences of CHRNA4 transcripts
hschrna4_ex_seqs
## DNAStringSetList of length 7
## [["uc002yes.4"]] AGCCCGGCGCTCCCTGCCGCGCCGCCGCCGCACCGCGCCCCACAGGAGAAGACGA...
## [["uc010gke.2"]] TCCCAGCCGGCTGAGGCGGGCAGGGCCGGGCGGGGCCGCGCCACGGAGCCCACAG...
## [["uc061ymb.1"]] GGCCCAGGGTTTGAAGCAGGTGAAAGTCTGAGCTACTTCTGCAAGTGCAGCCTTT...
## [["uc061ymd.1"]] CCAGCAGCCATGTGGAGACCCGGGCCCACGCCGAGGAGCGGCTCCTGAAGAAACT...
## [["uc061yme.1"]] CCACAGGAGAAGACGAACCGGGCCCGGCGGCCGAAGCGGCCCGCGAGGCGCGGGA...
## [["uc061ymf.1"]] CCACAGGAGAAGACGAACCGGGCCCGGCGGCCGAAGCGGCCCGCGAGGCGCGGGA...
## [["uc061ymh.1"]] CTAGAGCCCGCGAGGTGCGTGCGCCATGGAGCTAGGGGGCCCCGGAGCGCCGCGG...
hschrna4_ex_seqs[["uc002yes.4"]] # extract exon sequences of one CHRNA4 transcript
## A DNAStringSet instance of length 6
## width seq
## [1] 254 AGCCCGGCGCTCCCTGCCGCGCCGCCGCCGC...CTGCTGCTTCTGGGGACCGGCCTCCTGCGCG
## [2] 152 CCAGCAGCCATGTGGAGACCCGGGCCCACGC...CGGCCTGTCCATCGCTCAGCTCATTGACGTG
## [3] 45 GATGAGAAGAACCAGATGATGACCACGAACGTATGGGTGAAGCAG
## [4] 110 GAGTGGCACGACTACAAGCTGCGCTGGGACC...TCTGGCGGCCGGACATCGTCCTCTACAACAA
## [5] 1375 TGCTGACGGGGACTTCGCGGTCACCCACCTG...CCACCTGAAGGCCGAAGACACAGACTTCTCG
## [6] 3641 GTGAAGGAGGACTGGAAGTACGTGGCCATGG...AGACTTCAATGAATAAAACGCACCGTCACCA
hschrna4_in_seqs <- getSeq(Hsapiens, hschrna4_inbytx) # intron sequences of CHRNA4 gene/transcripts
hschrna4_in_seqs
## DNAStringSetList of length 7
## [["uc002yes.4"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc010gke.2"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc061ymb.1"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc061ymd.1"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc061yme.1"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc061ymf.1"]] GTAAGTCCCGCCCATGGCTGTGTTGGGCGCCTCTGACCCAGACGGAACCGGCGGA...
## [["uc061ymh.1"]] GTGAGTCCTCTACAGGGCCTTGGTGGGCAGGAGTGGCCAGGGCCTGCCCCTGCCC...
hschrna4_in_seqs[["uc002yes.4"]] # extract intron sequences of one isoform of CHRNA4 gene
## A DNAStringSet instance of length 5
## width seq
## [1] 2789 GTAAGTCCCGCCCATGGCTGTGTTGGGCGCC...CGTGCTGGAGTGACGGCGTCTGTGCTTGCAG
## [2] 4947 GTGAGTCCTCTACAACAAGTGAGTCCTCTAC...GTCTCACACCCTTCGCTCTCTTCCTGCCCAG
## [3] 286 GTGAGTGGGAGGGCACTGCCCTGCCCCACCC...CTGCAGCCTCCCCCCGCCCCACCTCATCCAG
## [4] 3132 GTAGGTGAGGGCGTGGCCATCGTGCACTGTG...CACTGGCCCCCTTCCTCTCTCGGTGCCCTAG
## [5] 1390 GTAAGTTACGGGGCGCGGGATGGGCCCCCTT...ACCTGAGCCACTGGCCTGCCTCTCCCCGCAG