Genomic annotation and analysis of human voltage-gated sodium channel subunits

Libraries

library(Biobase)

## 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

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

library(BiocGenerics)
library(GenomicFeatures)

## 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

library(IRanges)
library(GenomicRanges)
library(GenomeInfoDb)
library(AnnotationDbi)
library(MultiAssayExperiment)

## Loading required package: SummarizedExperiment

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## Loading required package: BiocParallel

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following objects are masked from 'package:base':
## 
##     aperm, apply

library(org.Hs.eg.db)

## 

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
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:DelayedArray':
## 
##     type

## The following object is masked from 'package:base':
## 
##     strsplit

## Loading required package: rtracklayer

Select organism database: Homo sapiens

org.Hs.eg.db

## OrgDb object:
## | DBSCHEMAVERSION: 2.1
## | Db type: OrgDb
## | Supporting package: AnnotationDbi
## | DBSCHEMA: HUMAN_DB
## | ORGANISM: Homo sapiens
## | SPECIES: Human
## | EGSOURCEDATE: 2018-Oct11
## | EGSOURCENAME: Entrez Gene
## | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | CENTRALID: EG
## | TAXID: 9606
## | GOSOURCENAME: Gene Ontology
## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
## | GOSOURCEDATE: 2018-Oct10
## | GOEGSOURCEDATE: 2018-Oct11
## | GOEGSOURCENAME: Entrez Gene
## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | KEGGSOURCENAME: KEGG GENOME
## | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
## | KEGGSOURCEDATE: 2011-Mar15
## | GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
## | GPSOURCEURL: 
## | GPSOURCEDATE: 2018-Oct2
## | ENSOURCEDATE: 2018-Oct05
## | ENSOURCENAME: Ensembl
## | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
## | UPSOURCENAME: Uniprot
## | UPSOURCEURL: http://www.UniProt.org/
## | UPSOURCEDATE: Thu Oct 18 05:22:10 2018

## 
## Please see: help('select') for usage information

columns(org.Hs.eg.db)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"

HsSCN <- select(org.Hs.eg.db,
                keys=c("SCN1A", "SCN2A", "SCN3A", "SCN4A", "SCN5A", "SCN7A", "SCN8A", "SCN9A",
                       "SCN10A", "SCN11A", "SCN1B", "SCN2B", "SCN3B", "SCN4B"),
                columns=columns(org.Hs.eg.db)[c(3, 4, 6, 10, 25)], 
                keytype="SYMBOL")

## 'select()' returned 1:many mapping between keys and columns

HsSCN

##    SYMBOL         ENSEMBL     ENSEMBLPROT ENTREZID
## 1   SCN1A ENSG00000144285            <NA>     6323
## 2   SCN1A ENSG00000144285            <NA>     6323
## 3   SCN1A ENSG00000144285            <NA>     6323
## 4   SCN1A ENSG00000144285            <NA>     6323
## 5   SCN2A ENSG00000136531            <NA>     6326
## 6   SCN3A ENSG00000153253            <NA>     6328
## 7   SCN4A ENSG00000007314            <NA>     6329
## 8   SCN5A ENSG00000183873            <NA>     6331
## 9   SCN7A ENSG00000136546            <NA>     6332
## 10  SCN7A ENSG00000136546            <NA>     6332
## 11  SCN8A ENSG00000196876            <NA>     6334
## 12  SCN8A ENSG00000196876            <NA>     6334
## 13  SCN9A ENSG00000169432            <NA>     6335
## 14 SCN10A ENSG00000185313 ENSP00000390600     6336
## 15 SCN10A ENSG00000185313 ENSP00000496110     6336
## 16 SCN10A ENSG00000185313 ENSP00000495595     6336
## 17 SCN10A ENSG00000185313 ENSP00000494789     6336
## 18 SCN11A ENSG00000168356            <NA>    11280
## 19  SCN1B ENSG00000105711            <NA>     6324
## 20  SCN2B ENSG00000149575 ENSP00000278947     6327
## 21  SCN3B ENSG00000166257 ENSP00000376523    55800
## 22  SCN3B ENSG00000166257 ENSP00000432785    55800
## 23  SCN3B ENSG00000166257 ENSP00000435554    55800
## 24  SCN3B ENSG00000166257 ENSP00000434363    55800
## 25  SCN3B ENSG00000166257 ENSP00000299333    55800
## 26  SCN4B ENSG00000177098 ENSP00000322460     6330
## 27  SCN4B ENSG00000177098 ENSP00000436343     6330
##                                         GENENAME   UNIGENE
## 1   sodium voltage-gated channel alpha subunit 1  Hs.22654
## 2   sodium voltage-gated channel alpha subunit 1 Hs.629873
## 3   sodium voltage-gated channel alpha subunit 1 Hs.693440
## 4   sodium voltage-gated channel alpha subunit 1 Hs.740081
## 5   sodium voltage-gated channel alpha subunit 2  Hs.93485
## 6   sodium voltage-gated channel alpha subunit 3 Hs.435274
## 7   sodium voltage-gated channel alpha subunit 4  Hs.46038
## 8   sodium voltage-gated channel alpha subunit 5 Hs.517898
## 9   sodium voltage-gated channel alpha subunit 7 Hs.596087
## 10  sodium voltage-gated channel alpha subunit 7 Hs.644853
## 11  sodium voltage-gated channel alpha subunit 8 Hs.436550
## 12  sodium voltage-gated channel alpha subunit 8 Hs.710638
## 13  sodium voltage-gated channel alpha subunit 9 Hs.439145
## 14 sodium voltage-gated channel alpha subunit 10 Hs.250443
## 15 sodium voltage-gated channel alpha subunit 10 Hs.250443
## 16 sodium voltage-gated channel alpha subunit 10 Hs.250443
## 17 sodium voltage-gated channel alpha subunit 10 Hs.250443
## 18 sodium voltage-gated channel alpha subunit 11 Hs.591657
## 19   sodium voltage-gated channel beta subunit 1 Hs.436646
## 20   sodium voltage-gated channel beta subunit 2 Hs.129783
## 21   sodium voltage-gated channel beta subunit 3   Hs.4865
## 22   sodium voltage-gated channel beta subunit 3   Hs.4865
## 23   sodium voltage-gated channel beta subunit 3   Hs.4865
## 24   sodium voltage-gated channel beta subunit 3   Hs.4865
## 25   sodium voltage-gated channel beta subunit 3   Hs.4865
## 26   sodium voltage-gated channel beta subunit 4  Hs.65239
## 27   sodium voltage-gated channel beta subunit 4  Hs.65239

table(HsSCN$ENTREZID)

## 
## 11280 55800  6323  6324  6326  6327  6328  6329  6330  6331  6332  6334 
##     1     5     4     1     1     1     1     1     2     1     2     2 
##  6335  6336 
##     1     4

unique(HsSCN$ENTREZID)

##  [1] "6323"  "6326"  "6328"  "6329"  "6331"  "6332"  "6334"  "6335" 
##  [9] "6336"  "11280" "6324"  "6327"  "55800" "6330"

According to ENTREZ/NCBI

“6323” = SCN1A “6326” = SCN2A “6328” = SCN3A “6329” = SCN4A “6331” = SCN5A “6332” = SCN7A “6334” = SCN8A “6335” = SCN9A “6336” = SCN10A “11280” = SCN11A “6324” = SCN1B “6327” = SCN2B “55800” = SCN3B “6330” = SCN4B

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

Extract the cordinates of genes, transcripts, exons and introns of SCN genes using hg38

HsSCNgenes = c("6323",  "6326",  "6328",  "6329",  "6331",  "6332",  "6334",  "6335",
                "6336",  "11280", "6324",  "6327", "55800", "6330")

hsscn_txs = transcriptsBy(txdb, by = "gene")[HsSCNgenes]
hsscn_txs #GRangesList object of length 14

## GRangesList object of length 14:
## $6323 
## GRanges object with 8 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 165989160-166073697      - |     27086  uc061pes.1
##   [2]     chr2 165989160-166128047      - |     27087  uc061pet.1
##   [3]     chr2 165989163-166073639      - |     27088  uc061peu.1
##   [4]     chr2 165991245-166073621      - |     27089  uc061pev.1
##   [5]     chr2 165994167-165996322      - |     27090  uc061pew.1
##   [6]     chr2 165996053-166000063      - |     27091  uc061pex.1
##   [7]     chr2 165998125-166002908      - |     27092  uc061pey.1
##   [8]     chr2 166073302-166128020      - |     27093  uc061pfd.1
## 
## $6326 
## GRanges object with 7 ranges and 2 metadata columns:
##       seqnames              ranges strand | tx_id    tx_name
##   [1]     chr2 165239402-165315604      + | 20181 uc061pdn.1
##   [2]     chr2 165239402-165392308      + | 20182 uc002udc.4
##   [3]     chr2 165239488-165392032      + | 20183 uc061pdo.1
##   [4]     chr2 165294031-165392310      + | 20184 uc002udd.4
##   [5]     chr2 165294044-165392051      + | 20185 uc061pdp.1
##   [6]     chr2 165295773-165392308      + | 20186 uc002ude.4
##   [7]     chr2 165307921-165310745      + | 20187 uc061pdq.1
## 
## ...
## <12 more elements>
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome

We can extract information about transcripts of each gene separately

Nav1.1 is coded by SCN1A/6323

hsscn_txs$`6323`

## GRanges object with 8 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 165989160-166073697      - |     27086  uc061pes.1
##   [2]     chr2 165989160-166128047      - |     27087  uc061pet.1
##   [3]     chr2 165989163-166073639      - |     27088  uc061peu.1
##   [4]     chr2 165991245-166073621      - |     27089  uc061pev.1
##   [5]     chr2 165994167-165996322      - |     27090  uc061pew.1
##   [6]     chr2 165996053-166000063      - |     27091  uc061pex.1
##   [7]     chr2 165998125-166002908      - |     27092  uc061pey.1
##   [8]     chr2 166073302-166128020      - |     27093  uc061pfd.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.2 is coded by SCN2A/6326

hsscn_txs$`6326`

## GRanges object with 7 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 165239402-165315604      + |     20181  uc061pdn.1
##   [2]     chr2 165239402-165392308      + |     20182  uc002udc.4
##   [3]     chr2 165239488-165392032      + |     20183  uc061pdo.1
##   [4]     chr2 165294031-165392310      + |     20184  uc002udd.4
##   [5]     chr2 165294044-165392051      + |     20185  uc061pdp.1
##   [6]     chr2 165295773-165392308      + |     20186  uc002ude.4
##   [7]     chr2 165307921-165310745      + |     20187  uc061pdq.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.3 is coded by SCN3A/6328

hsscn_txs$`6328`

## GRanges object with 8 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 165087522-165204067      - |     27059  uc061pdg.1
##   [2]     chr2 165087530-165204043      - |     27060  uc002ucx.4
##   [3]     chr2 165089718-165203886      - |     27061  uc002ucz.4
##   [4]     chr2 165092256-165093828      - |     27062  uc061pdi.1
##   [5]     chr2 165095511-165097614      - |     27063  uc061pdj.1
##   [6]     chr2 165097252-165186554      - |     27064  uc061pdk.1
##   [7]     chr2 165122800-165127792      - |     27065  uc061pdl.1
##   [8]     chr2 165164393-165204050      - |     27066  uc061pdm.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.4 is coded by SCN4A/6329

hsscn_txs$`6329`

## GRanges object with 4 ranges and 2 metadata columns:
##       seqnames            ranges strand |     tx_id     tx_name
##          <Rle>         <IRanges>  <Rle> | <integer> <character>
##   [1]    chr17 63938554-63972918      - |    149348  uc002jds.1
##   [2]    chr17 63938554-63972918      - |    149349  uc060iti.1
##   [3]    chr17 63947954-63949704      - |    149350  uc060itj.1
##   [4]    chr17 63959373-63961693      - |    149351  uc060itk.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.5 is coded by SCN5A/6331

hsscn_txs$`6331`

## GRanges object with 14 ranges and 2 metadata columns:
##        seqnames            ranges strand |     tx_id     tx_name
##           <Rle>         <IRanges>  <Rle> | <integer> <character>
##    [1]     chr3 38548057-38633349      - |     35902  uc062ihd.1
##    [2]     chr3 38548062-38633362      - |     35903  uc062ihe.1
##    [3]     chr3 38548062-38649672      - |     35904  uc062ihf.1
##    [4]     chr3 38548062-38649672      - |     35905  uc062ihg.1
##    [5]     chr3 38548062-38649673      - |     35906  uc062ihh.1
##    ...      ...               ...    ... .       ...         ...
##   [10]     chr3 38552981-38554549      - |     35911  uc062ihm.1
##   [11]     chr3 38614028-38649673      - |     35912  uc032rjy.2
##   [12]     chr3 38614543-38633349      - |     35913  uc062ihn.1
##   [13]     chr3 38629842-38633220      - |     35914  uc062iho.1
##   [14]     chr3 38633111-38645776      - |     35915  uc062ihp.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nax is coded by SCN7A/6323

hsscn_txs$`6332`

## GRanges object with 7 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 166403573-166486968      - |     27103  uc002udu.3
##   [2]     chr2 166403573-166486971      - |     27104  uc061pfp.1
##   [3]     chr2 166403573-166486971      - |     27105  uc061pfq.1
##   [4]     chr2 166405513-166486968      - |     27106  uc061pfr.1
##   [5]     chr2 166432321-166494247      - |     27107  uc061pfs.1
##   [6]     chr2 166456953-166494205      - |     27108  uc061pft.1
##   [7]     chr2 166461755-166465526      - |     27109  uc061pfu.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.6 is coded by SCN8A/6334

hsscn_txs$`6334`

## GRanges object with 10 ranges and 2 metadata columns:
##        seqnames            ranges strand |     tx_id     tx_name
##           <Rle>         <IRanges>  <Rle> | <integer> <character>
##    [1]    chr12 51590266-51662820      + |    104501  uc058ody.1
##    [2]    chr12 51590926-51766251      + |    104502  uc001ryy.3
##    [3]    chr12 51591236-51812864      + |    104503  uc001ryw.4
##    [4]    chr12 51591236-51812864      + |    104504  uc010snl.3
##    [5]    chr12 51662818-51807429      + |    104505  uc001ryx.2
##    [6]    chr12 51662818-51807429      + |    104506  uc058odz.1
##    [7]    chr12 51662818-51807429      + |    104507  uc058oea.1
##    [8]    chr12 51686423-51745948      + |    104508  uc058oeb.1
##    [9]    chr12 51768972-51770009      + |    104509  uc058oec.1
##   [10]    chr12 51774240-51786771      + |    104510  uc001rza.2
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.7 is coded by SCN9A/6335

hsscn_txs$`6335`

## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]     chr2 166195185-166375993      - |     27097  uc010fpl.4
##   [2]     chr2 166195194-166375987      - |     27098  uc002udr.2
##   [3]     chr2 166198672-166311756      - |     27099  uc061pfl.1
##   [4]     chr2 166281760-166306571      - |     27100  uc061pfm.1
##   [5]     chr2 166281760-166306571      - |     27101  uc061pfn.1
##   [6]     chr2 166302619-166305920      - |     27102  uc061pfo.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.8 is coded by SCN10A/6336

hsscn_txs$`6336`

## GRanges object with 1 range and 2 metadata columns:
##       seqnames            ranges strand |     tx_id     tx_name
##          <Rle>         <IRanges>  <Rle> | <integer> <character>
##   [1]     chr3 38696802-38794010      - |     35916  uc003ciq.4
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Nav1.9 is coded by SCN10A/11280

hsscn_txs$`11280`

## GRanges object with 3 ranges and 2 metadata columns:
##       seqnames            ranges strand |     tx_id     tx_name
##          <Rle>         <IRanges>  <Rle> | <integer> <character>
##   [1]     chr3 38845769-38950561      - |     35918  uc062ihs.1
##   [2]     chr3 38846694-38950362      - |     35919  uc010hhn.2
##   [3]     chr3 38849278-38950561      - |     35920  uc003cis.2
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

NavB1 is coded by SCN1B/6324

hsscn_txs$`6324`

## GRanges object with 4 ranges and 2 metadata columns:
##       seqnames            ranges strand |     tx_id     tx_name
##          <Rle>         <IRanges>  <Rle> | <integer> <character>
##   [1]    chr19 35030684-35040448      + |    157234  uc002nxp.4
##   [2]    chr19 35030733-35040073      + |    157235  uc010xsg.3
##   [3]    chr19 35030821-35034270      + |    157236  uc002nxo.3
##   [4]    chr19 35036899-35040448      + |    157239  uc060wye.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

NavB2 is coded by SCN2B/6327

hsscn_txs$`6327`

## GRanges object with 1 range and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]    chr11 118161951-118176673      - |    101863  uc001psf.3
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

NavB3 is coded by SCN3B/55800

hsscn_txs$`55800`

## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]    chr11 123629187-123641124      - |    102152  uc058ipv.1
##   [2]    chr11 123629189-123654604      - |    102153  uc001pzb.2
##   [3]    chr11 123629189-123654607      - |    102154  uc001pza.2
##   [4]    chr11 123631214-123654607      - |    102155  uc058ipw.1
##   [5]    chr11 123642467-123655006      - |    102156  uc058ipx.1
##   [6]    chr11 123642610-123655244      - |    102157  uc058ipy.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

NavB4 is coded by SCN4B/6330

hsscn_txs$`6330`

## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames              ranges strand |     tx_id     tx_name
##          <Rle>           <IRanges>  <Rle> | <integer> <character>
##   [1]    chr11 118133377-118145433      - |    101857  uc058hwk.1
##   [2]    chr11 118133377-118152820      - |    101858  uc001pse.4
##   [3]    chr11 118133379-118142873      - |    101859  uc058hwl.1
##   [4]    chr11 118136573-118141864      - |    101860  uc058hwm.1
##   [5]    chr11 118136961-118152888      - |    101861  uc010rxv.3
##   [6]    chr11 118141206-118145700      - |    101862  uc058hwn.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

Make a dataframe for SCN transcripts/genes

hsscn_txs_df <- data.frame(hsscn_txs)
hsscn_txs_df

##    group group_name seqnames     start       end  width strand  tx_id
## 1      1       6323     chr2 165989160 166073697  84538      -  27086
## 2      1       6323     chr2 165989160 166128047 138888      -  27087
## 3      1       6323     chr2 165989163 166073639  84477      -  27088
## 4      1       6323     chr2 165991245 166073621  82377      -  27089
## 5      1       6323     chr2 165994167 165996322   2156      -  27090
## 6      1       6323     chr2 165996053 166000063   4011      -  27091
## 7      1       6323     chr2 165998125 166002908   4784      -  27092
## 8      1       6323     chr2 166073302 166128020  54719      -  27093
## 9      2       6326     chr2 165239402 165315604  76203      +  20181
## 10     2       6326     chr2 165239402 165392308 152907      +  20182
## 11     2       6326     chr2 165239488 165392032 152545      +  20183
## 12     2       6326     chr2 165294031 165392310  98280      +  20184
## 13     2       6326     chr2 165294044 165392051  98008      +  20185
## 14     2       6326     chr2 165295773 165392308  96536      +  20186
## 15     2       6326     chr2 165307921 165310745   2825      +  20187
## 16     3       6328     chr2 165087522 165204067 116546      -  27059
## 17     3       6328     chr2 165087530 165204043 116514      -  27060
## 18     3       6328     chr2 165089718 165203886 114169      -  27061
## 19     3       6328     chr2 165092256 165093828   1573      -  27062
## 20     3       6328     chr2 165095511 165097614   2104      -  27063
## 21     3       6328     chr2 165097252 165186554  89303      -  27064
## 22     3       6328     chr2 165122800 165127792   4993      -  27065
## 23     3       6328     chr2 165164393 165204050  39658      -  27066
## 24     4       6329    chr17  63938554  63972918  34365      - 149348
## 25     4       6329    chr17  63938554  63972918  34365      - 149349
## 26     4       6329    chr17  63947954  63949704   1751      - 149350
## 27     4       6329    chr17  63959373  63961693   2321      - 149351
## 28     5       6331     chr3  38548057  38633349  85293      -  35902
## 29     5       6331     chr3  38548062  38633362  85301      -  35903
## 30     5       6331     chr3  38548062  38649672 101611      -  35904
## 31     5       6331     chr3  38548062  38649672 101611      -  35905
## 32     5       6331     chr3  38548062  38649673 101612      -  35906
## 33     5       6331     chr3  38548066  38649628 101563      -  35907
## 34     5       6331     chr3  38549128  38633332  84205      -  35908
## 35     5       6331     chr3  38549968  38633349  83382      -  35909
## 36     5       6331     chr3  38550321  38633316  82996      -  35910
## 37     5       6331     chr3  38552981  38554549   1569      -  35911
## 38     5       6331     chr3  38614028  38649673  35646      -  35912
## 39     5       6331     chr3  38614543  38633349  18807      -  35913
## 40     5       6331     chr3  38629842  38633220   3379      -  35914
## 41     5       6331     chr3  38633111  38645776  12666      -  35915
## 42     6       6332     chr2 166403573 166486968  83396      -  27103
## 43     6       6332     chr2 166403573 166486971  83399      -  27104
## 44     6       6332     chr2 166403573 166486971  83399      -  27105
## 45     6       6332     chr2 166405513 166486968  81456      -  27106
## 46     6       6332     chr2 166432321 166494247  61927      -  27107
## 47     6       6332     chr2 166456953 166494205  37253      -  27108
## 48     6       6332     chr2 166461755 166465526   3772      -  27109
## 49     7       6334    chr12  51590266  51662820  72555      + 104501
## 50     7       6334    chr12  51590926  51766251 175326      + 104502
## 51     7       6334    chr12  51591236  51812864 221629      + 104503
## 52     7       6334    chr12  51591236  51812864 221629      + 104504
## 53     7       6334    chr12  51662818  51807429 144612      + 104505
## 54     7       6334    chr12  51662818  51807429 144612      + 104506
## 55     7       6334    chr12  51662818  51807429 144612      + 104507
## 56     7       6334    chr12  51686423  51745948  59526      + 104508
## 57     7       6334    chr12  51768972  51770009   1038      + 104509
## 58     7       6334    chr12  51774240  51786771  12532      + 104510
## 59     8       6335     chr2 166195185 166375993 180809      -  27097
## 60     8       6335     chr2 166195194 166375987 180794      -  27098
## 61     8       6335     chr2 166198672 166311756 113085      -  27099
## 62     8       6335     chr2 166281760 166306571  24812      -  27100
## 63     8       6335     chr2 166281760 166306571  24812      -  27101
## 64     8       6335     chr2 166302619 166305920   3302      -  27102
## 65     9       6336     chr3  38696802  38794010  97209      -  35916
## 66    10      11280     chr3  38845769  38950561 104793      -  35918
## 67    10      11280     chr3  38846694  38950362 103669      -  35919
## 68    10      11280     chr3  38849278  38950561 101284      -  35920
## 69    11       6324    chr19  35030684  35040448   9765      + 157234
## 70    11       6324    chr19  35030733  35040073   9341      + 157235
## 71    11       6324    chr19  35030821  35034270   3450      + 157236
## 72    11       6324    chr19  35036899  35040448   3550      + 157239
## 73    12       6327    chr11 118161951 118176673  14723      - 101863
## 74    13      55800    chr11 123629187 123641124  11938      - 102152
## 75    13      55800    chr11 123629189 123654604  25416      - 102153
## 76    13      55800    chr11 123629189 123654607  25419      - 102154
## 77    13      55800    chr11 123631214 123654607  23394      - 102155
## 78    13      55800    chr11 123642467 123655006  12540      - 102156
## 79    13      55800    chr11 123642610 123655244  12635      - 102157
## 80    14       6330    chr11 118133377 118145433  12057      - 101857
## 81    14       6330    chr11 118133377 118152820  19444      - 101858
## 82    14       6330    chr11 118133379 118142873   9495      - 101859
## 83    14       6330    chr11 118136573 118141864   5292      - 101860
## 84    14       6330    chr11 118136961 118152888  15928      - 101861
## 85    14       6330    chr11 118141206 118145700   4495      - 101862
##       tx_name
## 1  uc061pes.1
## 2  uc061pet.1
## 3  uc061peu.1
## 4  uc061pev.1
## 5  uc061pew.1
## 6  uc061pex.1
## 7  uc061pey.1
## 8  uc061pfd.1
## 9  uc061pdn.1
## 10 uc002udc.4
## 11 uc061pdo.1
## 12 uc002udd.4
## 13 uc061pdp.1
## 14 uc002ude.4
## 15 uc061pdq.1
## 16 uc061pdg.1
## 17 uc002ucx.4
## 18 uc002ucz.4
## 19 uc061pdi.1
## 20 uc061pdj.1
## 21 uc061pdk.1
## 22 uc061pdl.1
## 23 uc061pdm.1
## 24 uc002jds.1
## 25 uc060iti.1
## 26 uc060itj.1
## 27 uc060itk.1
## 28 uc062ihd.1
## 29 uc062ihe.1
## 30 uc062ihf.1
## 31 uc062ihg.1
## 32 uc062ihh.1
## 33 uc062ihi.1
## 34 uc062ihj.1
## 35 uc062ihk.1
## 36 uc062ihl.1
## 37 uc062ihm.1
## 38 uc032rjy.2
## 39 uc062ihn.1
## 40 uc062iho.1
## 41 uc062ihp.1
## 42 uc002udu.3
## 43 uc061pfp.1
## 44 uc061pfq.1
## 45 uc061pfr.1
## 46 uc061pfs.1
## 47 uc061pft.1
## 48 uc061pfu.1
## 49 uc058ody.1
## 50 uc001ryy.3
## 51 uc001ryw.4
## 52 uc010snl.3
## 53 uc001ryx.2
## 54 uc058odz.1
## 55 uc058oea.1
## 56 uc058oeb.1
## 57 uc058oec.1
## 58 uc001rza.2
## 59 uc010fpl.4
## 60 uc002udr.2
## 61 uc061pfl.1
## 62 uc061pfm.1
## 63 uc061pfn.1
## 64 uc061pfo.1
## 65 uc003ciq.4
## 66 uc062ihs.1
## 67 uc010hhn.2
## 68 uc003cis.2
## 69 uc002nxp.4
## 70 uc010xsg.3
## 71 uc002nxo.3
## 72 uc060wye.1
## 73 uc001psf.3
## 74 uc058ipv.1
## 75 uc001pzb.2
## 76 uc001pza.2
## 77 uc058ipw.1
## 78 uc058ipx.1
## 79 uc058ipy.1
## 80 uc058hwk.1
## 81 uc001pse.4
## 82 uc058hwl.1
## 83 uc058hwm.1
## 84 uc010rxv.3
## 85 uc058hwn.1

dim(hsscn_txs_df)

## [1] 85  9

write.csv(hsscn_txs_df, file = 'hsscn_txs')

Exons in SCN transcripts/genes: We can do it for all genes one by one

Exons in SCN1A gene/transcrpt(s)

hsscn1a_exbytx <- exonsBy(txdb, "tx", use.names=TRUE)[hsscn_txs$`6323`$tx_name]
hsscn1a_exbytx #GRangesList object of length 8

## GRangesList object of length 8:
## $uc061pes.1 
## GRanges object with 26 ranges and 3 metadata columns:
##        seqnames              ranges strand |   exon_id   exon_name
##           <Rle>           <IRanges>  <Rle> | <integer> <character>
##    [1]     chr2 166073358-166073697      - |     82817        <NA>
##    [2]     chr2 166058568-166058688      - |     82811        <NA>
##    [3]     chr2 166056411-166056498      - |     82809        <NA>
##    [4]     chr2 166054638-166054766      - |     82808        <NA>
##    [5]     chr2 166052852-166052943      - |     82807        <NA>
##    ...      ...                 ...    ... .       ...         ...
##   [22]     chr2 165999723-165999776      - |     82786        <NA>
##   [23]     chr2 165998038-165998175      - |     82784        <NA>
##   [24]     chr2 165996013-165996117      - |     82781        <NA>
##   [25]     chr2 165994146-165994416      - |     82779        <NA>
##   [26]     chr2 165989160-165992422      - |     82776        <NA>
##        exon_rank
##        <integer>
##    [1]         1
##    [2]         2
##    [3]         3
##    [4]         4
##    [5]         5
##    ...       ...
##   [22]        22
##   [23]        23
##   [24]        24
##   [25]        25
##   [26]        26
## 
## ...
## <7 more elements>
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome

elementNROWS(hsscn1a_exbytx) #lists number of exons per transcripts/in SCN1A gene

## uc061pes.1 uc061pet.1 uc061peu.1 uc061pev.1 uc061pew.1 uc061pex.1 
##         26         28         26         26          2          3 
## uc061pey.1 uc061pfd.1 
##          3          3

barplot(elementNROWS(hsscn1a_exbytx))

Introns in SCN transcripts/genes: We can do it for all the SCN genes one by one

Introns in SCN1A gene/transcrpt(s)

hsscn1a_inbytx <- intronsByTranscript(txdb, use.names=TRUE)[hsscn_txs$`6323`$tx_name]
hsscn1a_inbytx #GRangesList object of length 8

## GRangesList object of length 8:
## $uc061pes.1 
## GRanges object with 25 ranges and 0 metadata columns:
##        seqnames              ranges strand
##           <Rle>           <IRanges>  <Rle>
##    [1]     chr2 165992423-165994145      -
##    [2]     chr2 165994417-165996012      -
##    [3]     chr2 165996118-165998037      -
##    [4]     chr2 165998176-165999722      -
##    [5]     chr2 165999777-166002471      -
##    ...      ...                 ...    ...
##   [21]     chr2 166051989-166052851      -
##   [22]     chr2 166052944-166054637      -
##   [23]     chr2 166054767-166056410      -
##   [24]     chr2 166056499-166058567      -
##   [25]     chr2 166058689-166073357      -
## 
## ...
## <7 more elements>
## -------
## seqinfo: 455 sequences (1 circular) from hg38 genome

elementNROWS(hsscn1a_inbytx)

## uc061pes.1 uc061pet.1 uc061peu.1 uc061pev.1 uc061pew.1 uc061pex.1 
##         25         27         25         25          1          2 
## uc061pey.1 uc061pfd.1 
##          2          2

barplot(elementNROWS(hsscn1a_inbytx))

BSgenome.Hsapiens.UCSC.hg38

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)

Extract exon sequences of transcript/gene(s): We can do it for all SCN genes

Extract exon sequences of transcript/gene(s) of SCN1A

Hsapiens <- BSgenome.Hsapiens.UCSC.hg38
hsscn1a_ex_seqs <- getSeq(Hsapiens, hsscn1a_exbytx) ## extract exon sequences of SCN1A transcripts
hsscn1a_ex_seqs #DNAStringSetList of length 8

## DNAStringSetList of length 8
## [["uc061pes.1"]] AAAATATCTGTATTCCTTTTATTTTAGGAATTTCATATGCAGAATAAATGGTAAT...
## [["uc061pet.1"]] CAAGCTATTTGCTGATTTGTATTAGGTACCATAGAGTGAGGCGAGGATGAAGCCG...
## [["uc061peu.1"]] AATGTGCAGGATGACAAGATGGAGCAAACAGTGCTTGTACCACCAGGACCTGACA...
## [["uc061pev.1"]] ATGGAGCAAACAGTGCTTGTACCACCAGGACCTGACAGCTTCAACTTCTTCACCA...
## [["uc061pew.1"]] TAAAGTCAATTTTGAATTATTTAGGGAATTAAATATTATCATACCTAAAGAGTAC...
## [["uc061pex.1"]] ATACTGCTTTAAATATCATGACACCATGAACTCTAAAGAATTGAAAAGTGCTTGT...
## [["uc061pey.1"]] TTTTAAGATTTTTTTGAACCTTGCTTTTACATATCCTAGAATAAATAGCATTGAT...
## [["uc061pfd.1"]] ACCATAGAGTGAGGCGAGGATGAAGCCGAGAGGATACTGCAGAGGTCTCTGGTGC...

hsscn1a_ex_seqs$uc061pes.1 # extract exon sequences of one SCN1A transcript

##   A DNAStringSet instance of length 26
##      width seq
##  [1]   340 AAAATATCTGTATTCCTTTTATTTTAGGAAT...GACCTGGACCCCTACTATATCAATAAGAAA
##  [2]   121 ACTTTTATAGTATTGAATAAAGGGAAGGCCA...AAATAGCTATTAAGATTTTGGTACATTCAT
##  [3]    88 TATTCAGCATGCTAATTATGTGCACTATTTT...TAACCCTCCTGATTGGACAAAGAATGTAGA
##  [4]   129 ATACACCTTCACAGGAATATATACTTTTGAA...CTGGCTCGATTTCACTGTCATTACATTTGC
##  [5]    92 GTACGTCACAGAGTTTGTGGACCTGGGCAAT...GAGCATTGAAGACGATTTCAGTCATTCCAG
##  ...   ... ...
## [22]    54 GCCACATTCAAAGGATGGATGGATATAATGTATGCAGCAGTTGATTCCAGAAAT
## [23]   138 GTGGAACTCCAGCCTAAGTATGAAGAAAGTC...ATAGATAATTTCAACCAGCAGAAAAAGAAG
## [24]   105 TTTGGAGGTCAAGACATCTTTATGACAGAAG...AAACCGCAAAAGCCTATACCTCGACCAGGA
## [25]   271 AACAAATTTCAAGGAATGGTCTTTGACTTCG...ATTTTGTGGTTGTCATTCTCTCCATTGTAG
## [26]  3263 GTATGTTTCTTGCCGAGCTGATAGAAAAGTA...AAACTAATAAAGATTACATTTTTTATTTTA

#A DNAStringSet instance of length 26

Extract intron sequences of transcript/gene(s): We can do it for all SCN genes

Extract intron sequences of transcript/gene(s) of SCN1A

hsscn1a_in_seqs <- getSeq(Hsapiens, hsscn1a_inbytx) # intron sequences of SCN1A gene/transcripts
hsscn1a_in_seqs #DNAStringSetList of length 8

## DNAStringSetList of length 8
## [["uc061pes.1"]] GTAAGAAATATTTAAAGTTCTTAAATTCAGTTAAATAAAAGTGAAAGCTGAAACA...
## [["uc061pet.1"]] GTAAGAAATATTTAAAGTTCTTAAATTCAGTTAAATAAAAGTGAAAGCTGAAACA...
## [["uc061peu.1"]] GTAAGAAATATTTAAAGTTCTTAAATTCAGTTAAATAAAAGTGAAAGCTGAAACA...
## [["uc061pev.1"]] GTAAGAAATATTTAAAGTTCTTAAATTCAGTTAAATAAAAGTGAAAGCTGAAACA...
## [["uc061pew.1"]] GTAAGAAGTATCAAATGATATGGGGGAAAAATACAAAAACAAAAACTTGCATGCT...
## [["uc061pex.1"]] ATAAGTATTTCTAATATTTTCTCTCCCACTGAGATAGAAAATTATTCCTTGGAGT...
## [["uc061pey.1"]] GTAAGTATTCCTTGTATTCTAAGTCTTTTTACAATATTGATCAGGTGGTAAAATT...
## [["uc061pfd.1"]] GTATGCAATTGCTGGATCATATGGTAAGAGTATGTTTGGTTTTTAAGAAACTGCT...

hsscn1a_in_seqs$uc061pes.1 ## extract intron sequences of one isoforms of SCN1A

##   A DNAStringSet instance of length 25
##      width seq
##  [1]  1723 GTAAGAAATATTTAAAGTTCTTAAATTCAGT...CTGGTTGGTATTCTCATTGTTTATTCATAG
##  [2]  1596 GTAAGAAGTATCAAATGATATGGGGGAAAAA...GTTGTTAAAAGCATTTCTATTTCTCTACAG
##  [3]  1920 ATAAGTATTTCTAATATTTTCTCTCCCACTG...TTTAACCAGTTTGATTTTTCTTTTCTATAC
##  [4]  1547 GTAAGTATTCCTTGTATTCTAAGTCTTTTTA...AATATTTTACAAAATATTCCCCTTTGGTAG
##  [5]  2695 GTAAGTGAACACTATTTTCTCTGAATATTTT...AAAGTAAAACATGCATGTCCTTCTTAATAG
##  ...   ... ...
## [21]   863 GTGAGAGCAAGGTTAGATAATGAGACGGACC...TCTTCATAGTAATAGACTTCTTTCTTTCAG
## [22]  1694 GTAAGTGCCTTTTTTGAAACTTTAAGAGAGA...TTGTGTTTGTGTGTGAACTCCCTATTACAG
## [23]  1644 GTAAGTTCAACTTATATTTTTAATAACATAT...TTATCAAATTTTTTGGATGCTTGTTTTCAG
## [24]  2069 ATCCTTTTTCAAGTGATTAATATTAACTATT...TTTTGTGTGTTCCTTAACAATAACCTACAT
## [25] 14669 GTGAGTGTTTTTTATCAGGCATATTTTGCTG...ATACTTCTATGTTGTGTTCCTGTCTTACAG

# A DNAStringSet instance of length 25