NGS Analysis Basics

34 minute read

Source code downloads:     [ .Rmd ]     [ .R ]

Overview

Sequence Analysis in R and Bioconductor

R Base

  • Some basic string handling utilities. Wide spectrum of numeric data analysis tools.

Bioconductor

Bioconductor packages provide much more sophisticated string handling utilities for sequence analysis (Lawrence et al. 2013; Huber et al. 2015).

Package Requirements

Several Bioconductor packages are required for this tutorial. To install them, execute the following lines in the R console. Please also make sure that you have a recent R version installed on your system. R versions 4.0.x or higher are recommended.

Please do not run this install on the HPCC unless you want to reinstall some of these packages in your own user account.

source("https://bioconductor.org/biocLite.R")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("Biostrings", "GenomicRanges", "rtracklayer", "systemPipeR", "seqLogo", "ShortRead"))

Strings in R Base

Basic String Matching and Parsing

String matching

Generate sample sequence data set

myseq <- c("ATGCAGACATAGTG", "ATGAACATAGATCC", "GTACAGATCAC")

String searching with regular expression support

myseq[grep("ATG", myseq)] 
## [1] "ATGCAGACATAGTG" "ATGAACATAGATCC"

Searches myseq for first match of pattern “AT”

pos1 <- regexpr("AT", myseq) 
as.numeric(pos1); attributes(pos1)$match.length # Returns position information of matches
## [1] 1 1 7

## [1] 2 2 2

Searches myseq for all matches of pattern “AT”

pos2 <- gregexpr("AT", myseq) 
as.numeric(pos2[[1]]); attributes(pos2[[1]])$match.length # Returns positions of matches in first sequence
## [1] 1 9

## [1] 2 2

String substitution with regular expression support

gsub("^ATG", "atg", myseq) 
## [1] "atgCAGACATAGTG" "atgAACATAGATCC" "GTACAGATCAC"

Positional parsing

nchar(myseq) # Computes length of strings
## [1] 14 14 11
substring(myseq[1], c(1,3), c(2,5)) # Positional parsing of several fragments from one string
## [1] "AT"  "GCA"
substring(myseq, c(1,4,7), c(2,6,10)) # Positional parsing of many strings
## [1] "AT"   "AAC"  "ATCA"

Random Sequence Generation

Random DNA sequences of any length

rand <- sapply(1:100, function(x) paste(sample(c("A","T","G","C"), sample(10:20), replace=TRUE), collapse=""))
rand[1:3]
## [1] "GCACACCATGTT"         "CCAAGGGGTCACGCAAGTAA" "GAATACGGAAAAG"

Count identical sequences

table(c(rand[1:4], rand[1]))
## 
## CCAAGGGGTCACGCAAGTAA        GAATACGGAAAAG         GCACACCATGTT       TAATCATCGAAGTT 
##                    1                    1                    2                    1

Extract reads from reference

Note: this requires the Biostrings package.

library(Biostrings)
ref <- DNAString(paste(sample(c("A","T","G","C"), 100000, replace=T), collapse=""))
randstart <- sample(1:(length(ref)-15), 1000)
randreads <- Views(ref, randstart, width=15)
rand_set <- DNAStringSet(randreads)
unlist(rand_set)
## 15000-letter DNAString object
## seq: AATCGCGGCACTTAGCCAGCTTGTCTTATCGGGTATTTTGAATCTT...ACTGTCGCCAATGGGTCCAACGGATTCTACGCCACGACTCGTAACT

Sequences in Bioconductor

Important Data Objects of Biostrings

XString for single sequence

  • DNAString: for DNA
  • RNAString: for RNA
  • AAString: for amino acid
  • BString: for any string

XStringSet for many sequences

  • `DNAStringSet``: for DNA
  • RNAStringSet: for RNA
  • AAStringSet: for amino acid
  • BStringSet: for any string

QualityScaleXStringSet for sequences with quality data

  • QualityScaledDNAStringSet: for DNA
  • QualityScaledRNAStringSet: for RNA
  • QualityScaledAAStringSet: for amino acid
  • QualityScaledBStringSet: for any string

Sequence Import and Export

Download the following sequences to your current working directory and then import them into R: https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_genbank/Bacteria/Halobacterium_sp_uid217/AE004437.ffn

dir.create("data", showWarnings = FALSE)
# system("wget https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_genbank/Bacteria/Halobacterium_sp_uid217/AE004437.ffn")
download.file("https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_genbank/Bacteria/Halobacterium_sp_uid217/AE004437.ffn", "data/AE004437.ffn")

Import FASTA file with readDNAStringSet

myseq <- readDNAStringSet("data/AE004437.ffn")
myseq[1:3]
## DNAStringSet object of length 3:
##     width seq                                                                   names               
## [1]  1206 ATGACTCGGCGGTCTCGTGTCGGTGCCGGCCTC...GTCGTCGTTGTTCGACGCTGGCGGAACCCATGA gi|12057215|gb|AE...
## [2]   666 ATGAGCATCATCGAACTCGAAGGCGTGGTCAAA...GTCAACCTCGTCGATGGGGTGTTACACACGTGA gi|12057215|gb|AE...
## [3]  1110 ATGGCGTGGCGGAACCTCGGGCGGAACCGCGTG...AACGATCCGCCCGTCGAGGCGCTCGGCGAATGA gi|12057215|gb|AE...

Subset sequences with regular expression on sequence name field

sub <- myseq[grep("99.*", names(myseq))]
length(sub)
## [1] 170

Export subsetted sequences to FASTA file

writeXStringSet(sub, file="./data/AE004437sub.ffn", width=80)

Now inspect exported sequence file AE004437sub.ffn in a text editor

Working with XString Containers

The XString stores the different types of biosequences in dedicated containers

library(Biostrings)
d <- DNAString("GCATAT-TAC")
d
## 10-letter DNAString object
## seq: GCATAT-TAC
d[1:4]
## 4-letter DNAString object
## seq: GCAT

RNA sequences

r <- RNAString("GCAUAU-UAC") 
r <- RNAString(d) # Converts d to RNAString object
r
## 10-letter RNAString object
## seq: GCAUAU-UAC

Protein sequences

p <- AAString("HCWYHH")
p
## 6-letter AAString object
## seq: HCWYHH

Any type of character strings

b <- BString("I store any set of characters. Other XString objects store only the IUPAC characters.")
b
## 85-letter BString object
## seq: I store any set of characters. Other XString objects store only the IUPAC characters.

Working with XStringSet Containers

XStringSet containers allow to store many biosequences in one object

dset <- DNAStringSet(c("GCATATTAC", "AATCGATCC", "GCATATTAC")) 
names(dset) <- c("seq1", "seq2", "seq3") # Assigns names
dset[1:2]
## DNAStringSet object of length 2:
##     width seq                                                                   names               
## [1]     9 GCATATTAC                                                             seq1
## [2]     9 AATCGATCC                                                             seq2

Important utilities for XStringSet containers

width(dset) # Returns the length of each sequences
## [1] 9 9 9
d <- dset[[1]] # The [[ subsetting operator returns a single entry as XString object
dset2 <- c(dset, dset) # Appends/concatenates two XStringSet objects
dsetchar <- as.character(dset) # Converts XStringSet to named vector 
dsetone <- unlist(dset) # Collapses many sequences to a single one stored in a DNAString container

Sequence subsetting by positions:

DNAStringSet(dset, start=c(1,2,3), end=c(4,8,5)) 
## DNAStringSet object of length 3:
##     width seq                                                                   names               
## [1]     4 GCAT                                                                  seq1
## [2]     7 ATCGATC                                                               seq2
## [3]     3 ATA                                                                   seq3

Multiple Alignment Class

The XMultipleAlignment class stores the different types of multiple sequence alignments:

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata",
              "msx2_mRNA.aln", package = "Biostrings"), format = "clustal")
origMAlign
## DNAMultipleAlignment with 8 rows and 2343 columns
##      aln                                                                        names               
## [1] -----TCCCGTCTCCGCAGCAAAAAAGTTTGAGTCG...TTGTCCAAACTCACAATTAAAAAAAAAAAAAAAAA gi|84452153|ref|N...
## [2] ------------------------------------...----------------------------------- gi|208431713|ref|...
## [3] ------------------------------------...----------------------------------- gi|118601823|ref|...
## [4] ----------------------AAAAGTTGGAGTCT...----------------------------------- gi|114326503|ref|...
## [5] ------------------------------------...----------------------------------- gi|119220589|ref|...
## [6] ------------------------------------...----------------------------------- gi|148540149|ref|...
## [7] --------------CGGCTCCGCAGCGCCTCACTCG...----------------------------------- gi|45383056|ref|N...
## [8] GGGGGAGACTTCAGAAGTTGTTGTCCTCTCCGCTGA...----------------------------------- gi|213515133|ref|...

Basic Sequence Manipulations

Reverse and Complement

randset <- DNAStringSet(rand)
complement(randset[1:2])
## DNAStringSet object of length 2:
##     width seq
## [1]    12 CGTGTGGTACAA
## [2]    20 GGTTCCCCAGTGCGTTCATT
reverse(randset[1:2])
## DNAStringSet object of length 2:
##     width seq
## [1]    12 TTGTACCACACG
## [2]    20 AATGAACGCACTGGGGAACC
reverseComplement(randset[1:2])
## DNAStringSet object of length 2:
##     width seq
## [1]    12 AACATGGTGTGC
## [2]    20 TTACTTGCGTGACCCCTTGG

Translate DNA into Protein

translate(randset[1:2])
## Warning in .Call2("DNAStringSet_translate", x, skip_code, dna_codes[codon_alphabet], : in 'x[[2]]':
## last 2 bases were ignored

## AAStringSet object of length 2:
##     width seq
## [1]     4 AHHV
## [2]     6 PRGHAS

Pattern Matching

Pattern matching with mismatches

Find pattern matches in reference

myseq1 <- readDNAStringSet("./data/AE004437.ffn") 
mypos <- matchPattern("ATGGTG", myseq1[[1]], max.mismatch=1) 

Count only the corresponding matches

countPattern("ATGGCT", myseq1[[1]], max.mismatch=1) 
## [1] 3

Count matches in many sequences

vcountPattern("ATGGCT", myseq1, max.mismatch=1)[1:20]
##  [1] 3 0 5 4 1 2 2 1 4 3 0 0 1 2 0 1 4 0 0 1

Results shown in DNAStringSet object

tmp <- c(DNAStringSet("ATGGTG"), DNAStringSet(mypos)) 

Return a consensus matrix for query and hits

consensusMatrix(tmp)[1:4,] 
##   [,1] [,2] [,3] [,4] [,5] [,6]
## A    3    0    0    0    0    0
## C    1    1    0    0    0    0
## G    0    0    4    4    1    4
## T    0    3    0    0    3    0

Find all pattern matches in reference

myvpos <- vmatchPattern("ATGGCT", myseq1, max.mismatch=1) 
myvpos # The results are stored as MIndex object.
## MIndex object of length 2058
## $`gi|12057215|gb|AE004437.1|:248-1453 Halobacterium sp. NRC-1, complete genome`
## IRanges object with 3 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         1         6         6
##   [2]       383       388         6
##   [3]       928       933         6
## 
## $`gi|12057215|gb|AE004437.1|:1450-2115 Halobacterium sp. NRC-1, complete genome`
## IRanges object with 0 ranges and 0 metadata columns:
##        start       end     width
##    <integer> <integer> <integer>
## 
## $`gi|12057215|gb|AE004437.1|:2145-3254 Halobacterium sp. NRC-1, complete genome`
## IRanges object with 5 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         1         6         6
##   [2]        94        99         6
##   [3]       221       226         6
##   [4]       535       540         6
##   [5]       601       606         6
## 
## ...
## <2055 more elements>
Views(myseq1[[1]], start(myvpos[[1]]), end(myvpos[[1]])) # Retrieves the result for single entry
## Views on a 1206-letter DNAString subject
## subject: ATGACTCGGCGGTCTCGTGTCGGTGCCGGCCTCGCAGCCATTGT...TTGCGATCGTCGTCGTCGTTGTTCGACGCTGGCGGAACCCATGA
## views:
##       start end width
##   [1]     1   6     6 [ATGACT]
##   [2]   383 388     6 [ATGGCA]
##   [3]   928 933     6 [ATGACT]

Return all matches

sapply(names(myseq1), function(x) 
       as.character(Views(myseq1[[x]], start(myvpos[[x]]), end(myvpos[[x]]))))[1:4] 

Pattern matching with regular expression support

myseq <- DNAStringSet(c("ATGCAGACATAGTG", "ATGAACATAGATCC", "GTACAGATCAC"))
myseq[grep("^ATG", myseq, perl=TRUE)] # String searching with regular expression support
## DNAStringSet object of length 2:
##     width seq
## [1]    14 ATGCAGACATAGTG
## [2]    14 ATGAACATAGATCC
pos1 <- regexpr("AT", myseq) # Searches 'myseq' for first match of pattern "AT"
as.numeric(pos1); attributes(pos1)$match.length # Returns position information of matches
## [1] 1 1 7

## [1] 2 2 2
pos2 <- gregexpr("AT", myseq) # Searches 'myseq' for all matches of pattern "AT"
as.numeric(pos2[[1]]); attributes(pos2[[1]])$match.length # Match positions in first sequence
## [1] 1 9

## [1] 2 2
DNAStringSet(gsub("^ATG", "NNN", myseq)) # String substitution with regular expression support
## DNAStringSet object of length 3:
##     width seq
## [1]    14 NNNCAGACATAGTG
## [2]    14 NNNAACATAGATCC
## [3]    11 GTACAGATCAC

PWM Viewing and Searching

library(seqLogo) 
pwm <- PWM(DNAStringSet(c("GCT", "GGT", "GCA"))) 
pwm
##        [,1]      [,2]      [,3]
## A 0.0000000 0.0000000 0.2312611
## C 0.0000000 0.3157205 0.0000000
## G 0.3685591 0.2312611 0.0000000
## T 0.0000000 0.0000000 0.3157205
seqLogo(t(t(pwm) * 1/colSums(pwm)))

The ggseqlogo package (manual) provides many customization options for plotting sequence logos. It also supports various alphabets including sequence logos for amino acid sequences.

library(ggplot2); library(ggseqlogo)
pwm <- PWM(DNAStringSet(c("GCT", "GGT", "GCA"))) 
ggseqlogo(pwm)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

Search sequence for PWM matches with score better than min.score

chr <- DNAString("AAAGCTAAAGGTAAAGCAAAA") 
matchPWM(pwm, chr, min.score=0.9) 
## Views on a 21-letter DNAString subject
## subject: AAAGCTAAAGGTAAAGCAAAA
## views:
##       start end width
##   [1]     4   6     3 [GCT]
##   [2]    10  12     3 [GGT]
##   [3]    16  18     3 [GCA]

NGS Sequences

Sequence and Quality Data: FASTQ Format

Four lines per sequence:

  1. ID
  2. Sequence
  3. ID
  4. Base call qualities (Phred scores) as ASCII characters

The following gives an example of 3 Illumina reads in a FASTQ file. The numbers at the beginning of each line are not part of the FASTQ format. They have been added solely for illustration purposes.

1. @SRR038845.3 HWI-EAS038:6:1:0:1938 length=36
2. CAACGAGTTCACACCTTGGCCGACAGGCCCGGGTAA
3. +SRR038845.3 HWI-EAS038:6:1:0:1938 length=36
4. BA@7>B=>:>>7@7@>>9=BAA?;>52;>:9=8.=A
1. @SRR038845.41 HWI-EAS038:6:1:0:1474 length=36
2. CCAATGATTTTTTTCCGTGTTTCAGAATACGGTTAA
3. +SRR038845.41 HWI-EAS038:6:1:0:1474 length=36
4. BCCBA@BB@BBBBAB@B9B@=BABA@A:@693:@B=
1. @SRR038845.53 HWI-EAS038:6:1:1:360 length=36
2. GTTCAAAAAGAACTAAATTGTGTCAATAGAAAACTC
3. +SRR038845.53 HWI-EAS038:6:1:1:360 length=36
4. BBCBBBBBB@@BAB?BBBBCBC>BBBAA8>BBBAA@

Sequence and Quality Data: QualityScaleXStringSet

Phred quality scores are integers from 0-50 that are stored as ASCII characters after adding 33. The basic R functions rawToChar and charToRaw can be used to interconvert among their representations.

Phred score interconversion

phred <- 1:9
phreda <- paste(sapply(as.raw((phred)+33), rawToChar), collapse="")
phreda
## [1] "\"#$%&'()*"
as.integer(charToRaw(phreda))-33 
## [1] 1 2 3 4 5 6 7 8 9

Construct QualityScaledDNAStringSet from scratch

dset <- DNAStringSet(sapply(1:100, function(x) paste(sample(c("A","T","G","C"), 20, replace=T), collapse=""))) # Creates random sample sequence.
myqlist <- lapply(1:100, function(x) sample(1:40, 20, replace=T)) # Creates random Phred score list.
myqual <- sapply(myqlist, function(x) toString(PhredQuality(x))) # Converts integer scores into ASCII characters.
myqual <- PhredQuality(myqual) # Converts to a PhredQuality object.
dsetq1 <- QualityScaledDNAStringSet(dset, myqual) # Combines DNAStringSet and quality data in QualityScaledDNAStringSet object.
dsetq1[1:2]
##   A QualityScaledDNAStringSet instance containing:
## 
## DNAStringSet object of length 2:
##     width seq
## [1]    20 TAAAGGGCCCTACGCCTTGC
## [2]    20 TGCCTACACACAGTTAGTAC
## 
## PhredQuality object of length 2:
##     width seq
## [1]    20 DE4<,-D;7:E94G724C1E
## [2]    20 H#)I"/+2)3>/I@?9'8/)

Processing FASTQ Files with ShortRead

The following explains the basic usage of ShortReadQ objects. To make the sample code work, download and unzip this file to your current working directory. The following code performs the download for you.

library(ShortRead)
download.file("http://cluster.hpcc.ucr.edu/~tgirke/HTML_Presentations/Manuals/testdata/samplefastq/data.zip", "data.zip")
unzip("data.zip")

Important utilities for accessing FASTQ files

fastq <- list.files("data", "*.fastq$"); fastq <- paste("data/", fastq, sep="")
names(fastq) <- paste("flowcell6_lane", 1:length(fastq), sep="_") 
(fq <- readFastq(fastq[1])) # Imports first FASTQ file
## class: ShortReadQ
## length: 1000 reads; width: 36 cycles
countLines(dirPath="./data", pattern=".fastq$")/4 # Counts numbers of reads in FASTQ files
## SRR038845.fastq SRR038846.fastq SRR038848.fastq SRR038850.fastq 
##            1000            1000            1000            1000
id(fq)[1] # Returns ID field
## BStringSet object of length 1:
##     width seq
## [1]    43 SRR038845.3 HWI-EAS038:6:1:0:1938 length=36
sread(fq)[1] # Returns sequence
## DNAStringSet object of length 1:
##     width seq
## [1]    36 CAACGAGTTCACACCTTGGCCGACAGGCCCGGGTAA
quality(fq)[1] # Returns Phred scores 
## class: FastqQuality
## quality:
## BStringSet object of length 1:
##     width seq
## [1]    36 BA@7>B=>:>>7@7@>>9=BAA?;>52;>:9=8.=A
as(quality(fq), "matrix")[1:4,1:12] # Coerces Phred scores to numeric matrix
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]   33   32   31   22   29   33   28   29   25    29    29    22
## [2,]   33   34   34   33   32   31   33   33   31    33    33    33
## [3,]   33   33   34   33   33   33   33   33   33    31    31    33
## [4,]   33   33   33   33   31   33   28   31   28    32    33    33
ShortReadQ(sread=sread(fq), quality=quality(fq), id=id(fq)) # Constructs a ShortReadQ from components
## class: ShortReadQ
## length: 1000 reads; width: 36 cycles

FASTQ Quality Reports

Using systemPipeR

The following seeFastq and seeFastqPlot functions generate and plot a series of useful quality statistics for a set of FASTQ files.

library(systemPipeR)
fqlist <- seeFastq(fastq=fastq, batchsize=800, klength=8) # For real data set batchsize to at least 10^5 
seeFastqPlot(fqlist)

Handles many samples in one PDF file. For more details see here

Using ShortRead

The ShortRead package contains several FASTQ quality reporting functions.

sp <- SolexaPath(system.file('extdata', package='ShortRead'))
fl <- file.path(analysisPath(sp), "s_1_sequence.txt") 
fls <- c(fl, fl) 
coll <- QACollate(QAFastqSource(fls), QAReadQuality(), QAAdapterContamination(), 
        QANucleotideUse(), QAQualityUse(), QASequenceUse(), QAFrequentSequence(n=10), 
        QANucleotideByCycle(), QAQualityByCycle())
x <- qa2(coll, verbose=TRUE)
res <- report(x)
if(interactive())
browseURL(res) 

Filtering and Trimming FASTQ Files with ShortRead

Adaptor trimming

fqtrim <- trimLRPatterns(Rpattern="GCCCGGGTAA", subject=fq)
sread(fq)[1:2] # Before trimming
## DNAStringSet object of length 2:
##     width seq
## [1]    36 CAACGAGTTCACACCTTGGCCGACAGGCCCGGGTAA
## [2]    36 CCAATGATTTTTTTCCGTGTTTCAGAATACGGTTAA
sread(fqtrim)[1:2] # After trimming
## DNAStringSet object of length 2:
##     width seq
## [1]    26 CAACGAGTTCACACCTTGGCCGACAG
## [2]    36 CCAATGATTTTTTTCCGTGTTTCAGAATACGGTTAA

Read counting and duplicate removal

tables(fq)$distribution # Counts read occurences
##   nOccurrences nReads
## 1            1    948
## 2            2     26
sum(srduplicated(fq)) # Identifies duplicated reads
## [1] 26
fq[!srduplicated(fq)]
## class: ShortReadQ
## length: 974 reads; width: 36 cycles

Trimming low quality tails

cutoff <- 30
cutoff <- rawToChar(as.raw(cutoff+33))
sread(trimTails(fq, k=2, a=cutoff, successive=FALSE))[1:2]
## DNAStringSet object of length 2:
##     width seq
## [1]     4 CAAC
## [2]    20 CCAATGATTTTTTTCCGTGT

Removal of reads with Phred scores below a threshold value

cutoff <- 30
qcount <- rowSums(as(quality(fq), "matrix") <= 20) 
fq[qcount == 0] # Number of reads where all Phred scores >= 20
## class: ShortReadQ
## length: 349 reads; width: 36 cycles

Removal of reads with x Ns and/or low complexity segments

filter1 <- nFilter(threshold=1) # Keeps only reads without Ns
filter2 <- polynFilter(threshold=20, nuc=c("A","T","G","C")) # Removes reads with nucleotide bias, >=20 of any base
filter <- compose(filter1, filter2)
fq[filter(fq)]
## class: ShortReadQ
## length: 989 reads; width: 36 cycles

Memory Efficient FASTQ Processing

Streaming through FASTQ files with FastqStreamer and random sampling reads with FastqSampler

fq <- yield(FastqStreamer(fastq[1], 50)) # Imports first 50 reads 
fq <- yield(FastqSampler(fastq[1], 50)) # Random samples 50 reads 

Streaming through a FASTQ file while applying filtering/trimming functions and writing the results to a new file here SRR038845.fastq_sub in data directory.

f <- FastqStreamer(fastq[1], 50) 
while(length(fq <- yield(f))) {
    fqsub <- fq[grepl("^TT", sread(fq))] 
    writeFastq(fqsub, paste(fastq[1], "sub", sep="_"), mode="a", compress=FALSE)
}
close(f)

Range Operations

Important Data Objects for Range Operations

  • IRanges: stores range data only (IRanges library)
  • GRanges: stores ranges and annotations (GenomicRanges library)
  • GRangesList: list version of GRanges container (GenomicRanges library)

Range Data Are Stored in IRanges and GRanges Containers

Construct GRanges Object

library(GenomicRanges); library(rtracklayer)
gr <- GRanges(seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)), ranges = IRanges(1:10, end = 7:16, names = head(letters, 10)), strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)), score = 1:10, GC = seq(1, 0, length = 10)) # Example of creating a GRanges object with its constructor function.

Import GFF into GRanges Object

gff <- import.gff("http://cluster.hpcc.ucr.edu/~tgirke/Documents/R_BioCond/Samples/gff3.gff") # Imports a simplified GFF3 genome annotation file.
seqlengths(gff) <- end(ranges(gff[which(values(gff)[,"type"]=="chromosome"),])) 
names(gff) <- 1:length(gff) # Assigns names to corresponding slot
gff[1:4,]
## GRanges object with 4 ranges and 10 metadata columns:
##     seqnames     ranges strand |   source       type     score     phase                  ID
##        <Rle>  <IRanges>  <Rle> | <factor>   <factor> <numeric> <integer>         <character>
##   1     Chr1 1-30427671      + |   TAIR10 chromosome        NA      <NA>                Chr1
##   2     Chr1  3631-5899      + |   TAIR10 gene              NA      <NA>           AT1G01010
##   3     Chr1  3631-5899      + |   TAIR10 mRNA              NA      <NA>         AT1G01010.1
##   4     Chr1  3760-5630      + |   TAIR10 protein           NA      <NA> AT1G01010.1-Protein
##            Name                Note          Parent       Index Derives_from
##     <character>     <CharacterList> <CharacterList> <character>  <character>
##   1        Chr1                                            <NA>         <NA>
##   2   AT1G01010 protein_coding_gene                        <NA>         <NA>
##   3 AT1G01010.1                           AT1G01010           1         <NA>
##   4 AT1G01010.1                                            <NA>  AT1G01010.1
##   -------
##   seqinfo: 7 sequences from an unspecified genome
seqinfo(gff)
## Seqinfo object with 7 sequences from an unspecified genome:
##   seqnames seqlengths isCircular genome
##   Chr1       30427671         NA   <NA>
##   Chr2       19698289         NA   <NA>
##   Chr3       23459830         NA   <NA>
##   Chr4       18585056         NA   <NA>
##   Chr5       26975502         NA   <NA>
##   ChrC         154478         NA   <NA>
##   ChrM         366924         NA   <NA>

Coerce GRanges object to data.frame

as.data.frame(gff)[1:4, 1:7]
##   seqnames start      end    width strand source       type
## 1     Chr1     1 30427671 30427671      + TAIR10 chromosome
## 2     Chr1  3631     5899     2269      + TAIR10       gene
## 3     Chr1  3631     5899     2269      + TAIR10       mRNA
## 4     Chr1  3760     5630     1871      + TAIR10    protein

Utilities for Range Containers

Accessor and subsetting methods for GRanges objects

Subsetting and replacement

gff[1:4]
## GRanges object with 4 ranges and 10 metadata columns:
##     seqnames     ranges strand |   source       type     score     phase                  ID
##        <Rle>  <IRanges>  <Rle> | <factor>   <factor> <numeric> <integer>         <character>
##   1     Chr1 1-30427671      + |   TAIR10 chromosome        NA      <NA>                Chr1
##   2     Chr1  3631-5899      + |   TAIR10 gene              NA      <NA>           AT1G01010
##   3     Chr1  3631-5899      + |   TAIR10 mRNA              NA      <NA>         AT1G01010.1
##   4     Chr1  3760-5630      + |   TAIR10 protein           NA      <NA> AT1G01010.1-Protein
##            Name                Note          Parent       Index Derives_from
##     <character>     <CharacterList> <CharacterList> <character>  <character>
##   1        Chr1                                            <NA>         <NA>
##   2   AT1G01010 protein_coding_gene                        <NA>         <NA>
##   3 AT1G01010.1                           AT1G01010           1         <NA>
##   4 AT1G01010.1                                            <NA>  AT1G01010.1
##   -------
##   seqinfo: 7 sequences from an unspecified genome
gff[1:4, c("type", "ID")] 
## GRanges object with 4 ranges and 2 metadata columns:
##     seqnames     ranges strand |       type                  ID
##        <Rle>  <IRanges>  <Rle> |   <factor>         <character>
##   1     Chr1 1-30427671      + | chromosome                Chr1
##   2     Chr1  3631-5899      + | gene                 AT1G01010
##   3     Chr1  3631-5899      + | mRNA               AT1G01010.1
##   4     Chr1  3760-5630      + | protein    AT1G01010.1-Protein
##   -------
##   seqinfo: 7 sequences from an unspecified genome
gff[2] <- gff[3] 

GRanges objects can be concatenated with the c function

c(gff[1:2], gff[401:402]) 
## GRanges object with 4 ranges and 10 metadata columns:
##       seqnames     ranges strand |   source           type     score     phase                  ID
##          <Rle>  <IRanges>  <Rle> | <factor>       <factor> <numeric> <integer>         <character>
##     1     Chr1 1-30427671      + |   TAIR10 chromosome            NA      <NA>                Chr1
##     2     Chr1  3631-5899      + |   TAIR10 mRNA                  NA      <NA>         AT1G01010.1
##   401     Chr5  5516-5769      - |   TAIR10 protein               NA      <NA> AT5G01015.2-Protein
##   402     Chr5  5770-5801      - |   TAIR10 five_prime_UTR        NA      <NA>                <NA>
##              Name            Note          Parent       Index Derives_from
##       <character> <CharacterList> <CharacterList> <character>  <character>
##     1        Chr1                                        <NA>         <NA>
##     2 AT1G01010.1                       AT1G01010           1         <NA>
##   401 AT5G01015.2                                        <NA>  AT5G01015.2
##   402        <NA>                     AT5G01015.2        <NA>         <NA>
##   -------
##   seqinfo: 7 sequences from an unspecified genome

Acessor functions

seqnames(gff)
## factor-Rle of length 449 with 7 runs
##   Lengths:   72   22   38  118  172   13   14
##   Values : Chr1 Chr2 Chr3 Chr4 Chr5 ChrC ChrM
## Levels(7): Chr1 Chr2 Chr3 Chr4 Chr5 ChrC ChrM
ranges(gff)
## IRanges object with 449 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##     1         1  30427671  30427671
##     2      3631      5899      2269
##     3      3631      5899      2269
##     4      3760      5630      1871
##     5      3631      3913       283
##   ...       ...       ...       ...
##   445     11918     12241       324
##   446     11918     12241       324
##   447     11918     12241       324
##   448     11918     12241       324
##   449     11918     12241       324
strand(gff)
## factor-Rle of length 449 with 13 runs
##   Lengths:  18  54  28  21  12 117   1 171   1  12   1   8   5
##   Values :   +   -   +   -   +   -   +   -   +   -   +   -   +
## Levels(3): + - *
seqlengths(gff) 
##     Chr1     Chr2     Chr3     Chr4     Chr5     ChrC     ChrM 
## 30427671 19698289 23459830 18585056 26975502   154478   366924
start(gff[1:4])
## [1]    1 3631 3631 3760
end(gff[1:4])
## [1] 30427671     5899     5899     5630
width(gff[1:4]) 
## [1] 30427671     2269     2269     1871

Accessing metadata component

values(gff) # or elementMetadata(gff)
## DataFrame with 449 rows and 10 columns
##       source       type     score     phase                  ID        Name                Note
##     <factor>   <factor> <numeric> <integer>         <character> <character>     <CharacterList>
## 1     TAIR10 chromosome        NA        NA                Chr1        Chr1                    
## 2     TAIR10 mRNA              NA        NA         AT1G01010.1 AT1G01010.1                    
## 3     TAIR10 mRNA              NA        NA         AT1G01010.1 AT1G01010.1                    
## 4     TAIR10 protein           NA        NA AT1G01010.1-Protein AT1G01010.1                    
## 5     TAIR10 exon              NA        NA                  NA          NA                    
## ...      ...        ...       ...       ...                 ...         ...                 ...
## 445   TAIR10    gene           NA        NA           ATMG00030   ATMG00030 protein_coding_gene
## 446   TAIR10    mRNA           NA        NA         ATMG00030.1 ATMG00030.1                    
## 447   TAIR10    protein        NA        NA ATMG00030.1-Protein ATMG00030.1                    
## 448   TAIR10    exon           NA        NA                  NA          NA                    
## 449   TAIR10    CDS            NA         0                  NA          NA                    
##                              Parent       Index Derives_from
##                     <CharacterList> <character>  <character>
## 1                                            NA           NA
## 2                         AT1G01010           1           NA
## 3                         AT1G01010           1           NA
## 4                                            NA  AT1G01010.1
## 5                       AT1G01010.1          NA           NA
## ...                             ...         ...          ...
## 445                                          NA           NA
## 446                       ATMG00030           1           NA
## 447                                          NA  ATMG00030.1
## 448                     ATMG00030.1          NA           NA
## 449 ATMG00030.1,ATMG00030.1-Protein          NA           NA
values(gff)[, "type"][1:20] 
##  [1] chromosome      mRNA            mRNA            protein         exon            five_prime_UTR 
##  [7] CDS             exon            CDS             exon            CDS             exon           
## [13] CDS             exon            CDS             exon            CDS             three_prime_UTR
## [19] gene            mRNA           
## Levels: chromosome gene mRNA protein exon five_prime_UTR CDS three_prime_UTR rRNA tRNA
gff[values(gff)[ ,"type"] == "gene"] 
## GRanges object with 21 ranges and 10 metadata columns:
##       seqnames      ranges strand |   source     type     score     phase          ID        Name
##          <Rle>   <IRanges>  <Rle> | <factor> <factor> <numeric> <integer> <character> <character>
##    19     Chr1   5928-8737      - |   TAIR10     gene        NA      <NA>   AT1G01020   AT1G01020
##    64     Chr1 11649-13714      - |   TAIR10     gene        NA      <NA>   AT1G01030   AT1G01030
##    74     Chr2   1025-2810      + |   TAIR10     gene        NA      <NA>   AT2G01008   AT2G01008
##    84     Chr2   3706-5513      + |   TAIR10     gene        NA      <NA>   AT2G01010   AT2G01010
##    87     Chr2   5782-5945      + |   TAIR10     gene        NA      <NA>   AT2G01020   AT2G01020
##   ...      ...         ...    ... .      ...      ...       ...       ...         ...         ...
##   427     ChrC    383-1444      - |   TAIR10     gene        NA      <NA>   ATCG00020   ATCG00020
##   432     ChrC   1717-4347      - |   TAIR10     gene        NA      <NA>   ATCG00030   ATCG00030
##   437     ChrM     273-734      - |   TAIR10     gene        NA      <NA>   ATMG00010   ATMG00010
##   442     ChrM  8848-11415      - |   TAIR10     gene        NA      <NA>   ATMG00020   ATMG00020
##   445     ChrM 11918-12241      + |   TAIR10     gene        NA      <NA>   ATMG00030   ATMG00030
##                      Note          Parent       Index Derives_from
##           <CharacterList> <CharacterList> <character>  <character>
##    19 protein_coding_gene                        <NA>         <NA>
##    64 protein_coding_gene                        <NA>         <NA>
##    74 protein_coding_gene                        <NA>         <NA>
##    84                rRNA                        <NA>         <NA>
##    87                rRNA                        <NA>         <NA>
##   ...                 ...             ...         ...          ...
##   427 protein_coding_gene                        <NA>         <NA>
##   432                tRNA                        <NA>         <NA>
##   437 protein_coding_gene                        <NA>         <NA>
##   442                rRNA                        <NA>         <NA>
##   445 protein_coding_gene                        <NA>         <NA>
##   -------
##   seqinfo: 7 sequences from an unspecified genome

Useful utilities for GRanges objects

Remove chromosome ranges

gff <- gff[values(gff)$type != "chromosome"] 

Erase the strand information

strand(gff) <- "*" 

Collapses overlapping ranges to continuous ranges.

reduce(gff) 
## GRanges object with 22 ranges and 0 metadata columns:
##        seqnames      ranges strand
##           <Rle>   <IRanges>  <Rle>
##    [1]     Chr1   3631-5899      *
##    [2]     Chr1   5928-8737      *
##    [3]     Chr1 11649-13714      *
##    [4]     Chr2   1025-2810      *
##    [5]     Chr2   3706-5513      *
##    ...      ...         ...    ...
##   [18]     ChrC    383-1444      *
##   [19]     ChrC   1717-4347      *
##   [20]     ChrM     273-734      *
##   [21]     ChrM  8848-11415      *
##   [22]     ChrM 11918-12241      *
##   -------
##   seqinfo: 7 sequences from an unspecified genome

Return uncovered regions

gaps(gff) 
## GRanges object with 43 ranges and 0 metadata columns:
##        seqnames       ranges strand
##           <Rle>    <IRanges>  <Rle>
##    [1]     Chr1   1-30427671      +
##    [2]     Chr1   1-30427671      -
##    [3]     Chr1       1-3630      *
##    [4]     Chr1    5900-5927      *
##    [5]     Chr1   8738-11648      *
##    ...      ...          ...    ...
##   [39]     ChrM     1-366924      -
##   [40]     ChrM        1-272      *
##   [41]     ChrM     735-8847      *
##   [42]     ChrM  11416-11917      *
##   [43]     ChrM 12242-366924      *
##   -------
##   seqinfo: 7 sequences from an unspecified genome

More intuitive way to get uncovered regions

setdiff(as(seqinfo(gff), "GRanges"), gff) 
## GRanges object with 29 ranges and 0 metadata columns:
##        seqnames         ranges strand
##           <Rle>      <IRanges>  <Rle>
##    [1]     Chr1         1-3630      *
##    [2]     Chr1      5900-5927      *
##    [3]     Chr1     8738-11648      *
##    [4]     Chr1 13715-30427671      *
##    [5]     Chr2         1-1024      *
##    ...      ...            ...    ...
##   [25]     ChrC    4348-154478      *
##   [26]     ChrM          1-272      *
##   [27]     ChrM       735-8847      *
##   [28]     ChrM    11416-11917      *
##   [29]     ChrM   12242-366924      *
##   -------
##   seqinfo: 7 sequences from an unspecified genome

Return disjoint ranges

disjoin(gff)
## GRanges object with 211 ranges and 0 metadata columns:
##         seqnames      ranges strand
##            <Rle>   <IRanges>  <Rle>
##     [1]     Chr1   3631-3759      *
##     [2]     Chr1   3760-3913      *
##     [3]     Chr1   3914-3995      *
##     [4]     Chr1   3996-4276      *
##     [5]     Chr1   4277-4485      *
##     ...      ...         ...    ...
##   [207]     ChrC   1752-4310      *
##   [208]     ChrC   4311-4347      *
##   [209]     ChrM     273-734      *
##   [210]     ChrM  8848-11415      *
##   [211]     ChrM 11918-12241      *
##   -------
##   seqinfo: 7 sequences from an unspecified genome

Returns coverage of ranges

coverage(gff)
## RleList of length 7
## $Chr1
## integer-Rle of length 30427671 with 45 runs
##   Lengths:     3630      129      154       82      281 ...      233      161      380 30413957
##   Values :        0        4        5        3        5 ...        4        2        4        0
## 
## $Chr2
## integer-Rle of length 19698289 with 14 runs
##   Lengths:     1024      248      185       53      362 ...      164      625      102 19691617
##   Values :        0        5        3        5        3 ...        3        0        5        0
## 
## $Chr3
## integer-Rle of length 23459830 with 29 runs
##   Lengths:     1652      145      139      111       95 ...      155      148      156 23453781
##   Values :        0        4        5        3        5 ...        3        5        4        0
## 
## $Chr4
## integer-Rle of length 18585056 with 72 runs
##   Lengths:     1179      357     1358      128      872 ...      212      114       74 18571697
##   Values :        0        5        0        5        3 ...        3        5        4        0
## 
## $Chr5
## integer-Rle of length 26975502 with 64 runs
##   Lengths:     1222       28       28      109       72 ...       76       55      174 26967058
##   Values :        0        4        7       13       16 ...        3        5        4        0
## 
## ...
## <2 more elements>

Return the index pairings for overlapping ranges

findOverlaps(gff, gff[1:4]) 
## Hits object with 55 hits and 0 metadata columns:
##        queryHits subjectHits
##        <integer>   <integer>
##    [1]         1           1
##    [2]         1           2
##    [3]         1           4
##    [4]         1           3
##    [5]         2           1
##    ...       ...         ...
##   [51]        16           1
##   [52]        16           2
##   [53]        16           3
##   [54]        17           1
##   [55]        17           2
##   -------
##   queryLength: 442 / subjectLength: 4

Counts overlapping ranges

countOverlaps(gff, gff[1:4])[1:40]
##  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 
##  4  4  4  4  3  4  3  3  3  3  3  3  3  3  3  3  2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 
## 35 36 37 38 39 40 41 
##  0  0  0  0  0  0  0

Return only overlapping ranges

subsetByOverlaps(gff, gff[1:4]) 
## GRanges object with 17 ranges and 10 metadata columns:
##      seqnames    ranges strand |   source            type     score     phase                  ID
##         <Rle> <IRanges>  <Rle> | <factor>        <factor> <numeric> <integer>         <character>
##    2     Chr1 3631-5899      * |   TAIR10  mRNA                  NA      <NA>         AT1G01010.1
##    3     Chr1 3631-5899      * |   TAIR10  mRNA                  NA      <NA>         AT1G01010.1
##    4     Chr1 3760-5630      * |   TAIR10  protein               NA      <NA> AT1G01010.1-Protein
##    5     Chr1 3631-3913      * |   TAIR10  exon                  NA      <NA>                <NA>
##    6     Chr1 3631-3759      * |   TAIR10  five_prime_UTR        NA      <NA>                <NA>
##   ..      ...       ...    ... .      ...             ...       ...       ...                 ...
##   14     Chr1 5174-5326      * |   TAIR10 exon                   NA      <NA>                <NA>
##   15     Chr1 5174-5326      * |   TAIR10 CDS                    NA         0                <NA>
##   16     Chr1 5439-5899      * |   TAIR10 exon                   NA      <NA>                <NA>
##   17     Chr1 5439-5630      * |   TAIR10 CDS                    NA         0                <NA>
##   18     Chr1 5631-5899      * |   TAIR10 three_prime_UTR        NA      <NA>                <NA>
##             Name            Note                          Parent       Index Derives_from
##      <character> <CharacterList>                 <CharacterList> <character>  <character>
##    2 AT1G01010.1                                       AT1G01010           1         <NA>
##    3 AT1G01010.1                                       AT1G01010           1         <NA>
##    4 AT1G01010.1                                                        <NA>  AT1G01010.1
##    5        <NA>                                     AT1G01010.1        <NA>         <NA>
##    6        <NA>                                     AT1G01010.1        <NA>         <NA>
##   ..         ...             ...                             ...         ...          ...
##   14        <NA>                                     AT1G01010.1        <NA>         <NA>
##   15        <NA>                 AT1G01010.1,AT1G01010.1-Protein        <NA>         <NA>
##   16        <NA>                                     AT1G01010.1        <NA>         <NA>
##   17        <NA>                 AT1G01010.1,AT1G01010.1-Protein        <NA>         <NA>
##   18        <NA>                                     AT1G01010.1        <NA>         <NA>
##   -------
##   seqinfo: 7 sequences from an unspecified genome

GRangesList Objects

sp <- split(gff, seq(along=gff)) # Stores every range in separate component of a GRangesList object
split(gff, seqnames(gff)) # Stores ranges of each chromosome in separate component.
## GRangesList object of length 7:
## $Chr1
## GRanges object with 71 ranges and 10 metadata columns:
##      seqnames      ranges strand |   source            type     score     phase                  ID
##         <Rle>   <IRanges>  <Rle> | <factor>        <factor> <numeric> <integer>         <character>
##    2     Chr1   3631-5899      * |   TAIR10  mRNA                  NA      <NA>         AT1G01010.1
##    3     Chr1   3631-5899      * |   TAIR10  mRNA                  NA      <NA>         AT1G01010.1
##    4     Chr1   3760-5630      * |   TAIR10  protein               NA      <NA> AT1G01010.1-Protein
##    5     Chr1   3631-3913      * |   TAIR10  exon                  NA      <NA>                <NA>
##    6     Chr1   3631-3759      * |   TAIR10  five_prime_UTR        NA      <NA>                <NA>
##   ..      ...         ...    ... .      ...             ...       ...       ...                 ...
##   68     Chr1 13335-13714      * |   TAIR10 exon                   NA      <NA>                <NA>
##   69     Chr1 12941-13173      * |   TAIR10 five_prime_UTR         NA      <NA>                <NA>
##   70     Chr1 11864-12940      * |   TAIR10 CDS                    NA         0                <NA>
##   71     Chr1 11649-11863      * |   TAIR10 three_prime_UTR        NA      <NA>                <NA>
##   72     Chr1 11649-13173      * |   TAIR10 exon                   NA      <NA>                <NA>
##             Name            Note                          Parent       Index Derives_from
##      <character> <CharacterList>                 <CharacterList> <character>  <character>
##    2 AT1G01010.1                                       AT1G01010           1         <NA>
##    3 AT1G01010.1                                       AT1G01010           1         <NA>
##    4 AT1G01010.1                                                        <NA>  AT1G01010.1
##    5        <NA>                                     AT1G01010.1        <NA>         <NA>
##    6        <NA>                                     AT1G01010.1        <NA>         <NA>
##   ..         ...             ...                             ...         ...          ...
##   68        <NA>                                     AT1G01030.1        <NA>         <NA>
##   69        <NA>                                     AT1G01030.1        <NA>         <NA>
##   70        <NA>                 AT1G01030.1,AT1G01030.1-Protein        <NA>         <NA>
##   71        <NA>                                     AT1G01030.1        <NA>         <NA>
##   72        <NA>                                     AT1G01030.1        <NA>         <NA>
##   -------
##   seqinfo: 7 sequences from an unspecified genome
## 
## ...
## <6 more elements>
unlist(sp) # Returns data as GRanges object
## GRanges object with 442 ranges and 10 metadata columns:
##           seqnames      ranges strand |   source           type     score     phase
##              <Rle>   <IRanges>  <Rle> | <factor>       <factor> <numeric> <integer>
##       1.2     Chr1   3631-5899      * |   TAIR10 mRNA                  NA      <NA>
##       2.3     Chr1   3631-5899      * |   TAIR10 mRNA                  NA      <NA>
##       3.4     Chr1   3760-5630      * |   TAIR10 protein               NA      <NA>
##       4.5     Chr1   3631-3913      * |   TAIR10 exon                  NA      <NA>
##       5.6     Chr1   3631-3759      * |   TAIR10 five_prime_UTR        NA      <NA>
##       ...      ...         ...    ... .      ...            ...       ...       ...
##   438.445     ChrM 11918-12241      * |   TAIR10        gene           NA      <NA>
##   439.446     ChrM 11918-12241      * |   TAIR10        mRNA           NA      <NA>
##   440.447     ChrM 11918-12241      * |   TAIR10        protein        NA      <NA>
##   441.448     ChrM 11918-12241      * |   TAIR10        exon           NA      <NA>
##   442.449     ChrM 11918-12241      * |   TAIR10        CDS            NA         0
##                            ID        Name                Note                          Parent
##                   <character> <character>     <CharacterList>                 <CharacterList>
##       1.2         AT1G01010.1 AT1G01010.1                                           AT1G01010
##       2.3         AT1G01010.1 AT1G01010.1                                           AT1G01010
##       3.4 AT1G01010.1-Protein AT1G01010.1                                                    
##       4.5                <NA>        <NA>                                         AT1G01010.1
##       5.6                <NA>        <NA>                                         AT1G01010.1
##       ...                 ...         ...                 ...                             ...
##   438.445           ATMG00030   ATMG00030 protein_coding_gene                                
##   439.446         ATMG00030.1 ATMG00030.1                                           ATMG00030
##   440.447 ATMG00030.1-Protein ATMG00030.1                                                    
##   441.448                <NA>        <NA>                                         ATMG00030.1
##   442.449                <NA>        <NA>                     ATMG00030.1,ATMG00030.1-Protein
##                 Index Derives_from
##           <character>  <character>
##       1.2           1         <NA>
##       2.3           1         <NA>
##       3.4        <NA>  AT1G01010.1
##       4.5        <NA>         <NA>
##       5.6        <NA>         <NA>
##       ...         ...          ...
##   438.445        <NA>         <NA>
##   439.446           1         <NA>
##   440.447        <NA>  ATMG00030.1
##   441.448        <NA>         <NA>
##   442.449        <NA>         <NA>
##   -------
##   seqinfo: 7 sequences from an unspecified genome
sp[1:4, "type"] # Subsetting of GRangesList objects is similar to GRanges objects.
## GRangesList object of length 4:
## $`1`
## GRanges object with 1 range and 1 metadata column:
##     seqnames    ranges strand |     type
##        <Rle> <IRanges>  <Rle> | <factor>
##   2     Chr1 3631-5899      * |     mRNA
##   -------
##   seqinfo: 7 sequences from an unspecified genome
## 
## $`2`
## GRanges object with 1 range and 1 metadata column:
##     seqnames    ranges strand |     type
##        <Rle> <IRanges>  <Rle> | <factor>
##   3     Chr1 3631-5899      * |     mRNA
##   -------
##   seqinfo: 7 sequences from an unspecified genome
## 
## $`3`
## GRanges object with 1 range and 1 metadata column:
##     seqnames    ranges strand |     type
##        <Rle> <IRanges>  <Rle> | <factor>
##   4     Chr1 3760-5630      * |  protein
##   -------
##   seqinfo: 7 sequences from an unspecified genome
## 
## $`4`
## GRanges object with 1 range and 1 metadata column:
##     seqnames    ranges strand |     type
##        <Rle> <IRanges>  <Rle> | <factor>
##   5     Chr1 3631-3913      * |     exon
##   -------
##   seqinfo: 7 sequences from an unspecified genome
lapply(sp[1:4], length) # Looping over GRangesList objects similar to lists
## $`1`
## [1] 1
## 
## $`2`
## [1] 1
## 
## $`3`
## [1] 1
## 
## $`4`
## [1] 1

Transcript Ranges

Storing annotation ranges in TranscriptDb databases makes many operations more robust and convenient.

library(GenomicFeatures)
download.file("http://cluster.hpcc.ucr.edu/~tgirke/Documents/R_BioCond/Samples/gff3.gff", "data/gff3.gff")
txdb <- makeTxDbFromGFF(file="data/gff3.gff", format="gff", dataSource="TAIR", organism="Arabidopsis thaliana")
## Warning in .extract_exons_from_GRanges(cds_IDX, gr, mcols0, tx_IDX, feature = "cds", : 163 CDS couldn't be linked to a transcript so were dropped (showing only the first 6):
##   seqid start  end strand   ID Name              Parent Parent_type
## 1  Chr1  3760 3913      + <NA> <NA> AT1G01010.1-Protein        <NA>
## 2  Chr1  3996 4276      + <NA> <NA> AT1G01010.1-Protein        <NA>
## 3  Chr1  4486 4605      + <NA> <NA> AT1G01010.1-Protein        <NA>
## 4  Chr1  4706 5095      + <NA> <NA> AT1G01010.1-Protein        <NA>
## 5  Chr1  5174 5326      + <NA> <NA> AT1G01010.1-Protein        <NA>
## 6  Chr1  5439 5630      + <NA> <NA> AT1G01010.1-Protein        <NA>
saveDb(txdb, file="./data/TAIR10.sqlite")
## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: TAIR
## # Organism: Arabidopsis thaliana
## # Taxonomy ID: 3702
## # miRBase build ID: NA
## # Genome: NA
## # Nb of transcripts: 28
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time: 2022-04-20 18:50:11 -0700 (Wed, 20 Apr 2022)
## # GenomicFeatures version at creation time: 1.46.1
## # RSQLite version at creation time: 2.2.9
## # DBSCHEMAVERSION: 1.2
txdb <- loadDb("./data/TAIR10.sqlite")
transcripts(txdb)
## GRanges object with 28 ranges and 2 metadata columns:
##        seqnames      ranges strand |     tx_id     tx_name
##           <Rle>   <IRanges>  <Rle> | <integer> <character>
##    [1]     Chr1   3631-5899      + |         1 AT1G01010.1
##    [2]     Chr1   5928-8737      - |         2 AT1G01020.1
##    [3]     Chr1   6790-8737      - |         3 AT1G01020.2
##    [4]     Chr1 11649-13714      - |         4 AT1G01030.1
##    [5]     Chr2   1025-2810      + |         5 AT2G01008.1
##    ...      ...         ...    ... .       ...         ...
##   [24]     ChrC    383-1444      - |        24 ATCG00020.1
##   [25]     ChrC   1717-4347      - |        25 ATCG00030.1
##   [26]     ChrM 11918-12241      + |        26 ATMG00030.1
##   [27]     ChrM     273-734      - |        27 ATMG00010.1
##   [28]     ChrM  8848-11415      - |        28 ATMG00020.1
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
transcriptsBy(txdb, by = "gene")
## GRangesList object of length 22:
## $AT1G01010
## GRanges object with 1 range and 2 metadata columns:
##       seqnames    ranges strand |     tx_id     tx_name
##          <Rle> <IRanges>  <Rle> | <integer> <character>
##   [1]     Chr1 3631-5899      + |         1 AT1G01010.1
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## $AT1G01020
## GRanges object with 2 ranges and 2 metadata columns:
##       seqnames    ranges strand |     tx_id     tx_name
##          <Rle> <IRanges>  <Rle> | <integer> <character>
##   [1]     Chr1 5928-8737      - |         2 AT1G01020.1
##   [2]     Chr1 6790-8737      - |         3 AT1G01020.2
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## $AT1G01030
## GRanges object with 1 range and 2 metadata columns:
##       seqnames      ranges strand |     tx_id     tx_name
##          <Rle>   <IRanges>  <Rle> | <integer> <character>
##   [1]     Chr1 11649-13714      - |         4 AT1G01030.1
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## ...
## <19 more elements>
exonsBy(txdb, by = "gene")
## GRangesList object of length 22:
## $AT1G01010
## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames    ranges strand |   exon_id   exon_name
##          <Rle> <IRanges>  <Rle> | <integer> <character>
##   [1]     Chr1 3631-3913      + |         1        <NA>
##   [2]     Chr1 3996-4276      + |         2        <NA>
##   [3]     Chr1 4486-4605      + |         3        <NA>
##   [4]     Chr1 4706-5095      + |         4        <NA>
##   [5]     Chr1 5174-5326      + |         5        <NA>
##   [6]     Chr1 5439-5899      + |         6        <NA>
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## $AT1G01020
## GRanges object with 12 ranges and 2 metadata columns:
##        seqnames    ranges strand |   exon_id   exon_name
##           <Rle> <IRanges>  <Rle> | <integer> <character>
##    [1]     Chr1 5928-6263      - |         7        <NA>
##    [2]     Chr1 6437-7069      - |         8        <NA>
##    [3]     Chr1 6790-7069      - |         9        <NA>
##    [4]     Chr1 7157-7232      - |        10        <NA>
##    [5]     Chr1 7157-7450      - |        11        <NA>
##    ...      ...       ...    ... .       ...         ...
##    [8]     Chr1 7762-7835      - |        14        <NA>
##    [9]     Chr1 7942-7987      - |        15        <NA>
##   [10]     Chr1 8236-8325      - |        16        <NA>
##   [11]     Chr1 8417-8464      - |        17        <NA>
##   [12]     Chr1 8571-8737      - |        18        <NA>
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## $AT1G01030
## GRanges object with 2 ranges and 2 metadata columns:
##       seqnames      ranges strand |   exon_id   exon_name
##          <Rle>   <IRanges>  <Rle> | <integer> <character>
##   [1]     Chr1 11649-13173      - |        19        <NA>
##   [2]     Chr1 13335-13714      - |        20        <NA>
##   -------
##   seqinfo: 7 sequences (2 circular) from an unspecified genome; no seqlengths
## 
## ...
## <19 more elements>

txdb from BioMart

Alternative sources for creating txdb databases are BioMart, Bioc annotation packages, UCSC, etc. The following shows how to create a txdb from BioMart.

library(GenomicFeatures); library("biomaRt")
txdb <- makeTxDbFromBiomart(biomart = "plants_mart", dataset = "athaliana_eg_gene", host="plants.ensembl.org")

The following steps are useful to find out what is availble in BioMart.

listMarts() # Lists BioMart databases
listMarts(host="plants.ensembl.org")
mymart <- useMart("plants_mart", host="plants.ensembl.org") # Select one, here plants_mart
listDatasets(mymart) # List datasets available in the selected BioMart database
mymart <- useMart("plants_mart", dataset="athaliana_eg_gene", host="plants.ensembl.org")
listAttributes(mymart) # List available features 
getBM(attributes=c("ensembl_gene_id", "description"), mart=mymart)[1:4,]

Efficient Sequence Parsing

getSeq

The following parses all annotation ranges provided by a GRanges object (e.g. gff) from a genome sequence stored in a local file.

gff <- gff[values(gff)$type != "chromosome"] # Remove chromosome ranges
rand <- DNAStringSet(sapply(unique(as.character(seqnames(gff))), function(x) paste(sample(c("A","T","G","C"), 200000, replace=T), collapse="")))
writeXStringSet(DNAStringSet(rand), "./data/test")
getSeq(FaFile("./data/test"), gff)
## DNAStringSet object of length 442:
##       width seq                                                                 names               
##   [1]  2269 GCCAGGTAGAATTCCAAAAATTGAAAGGCGTT...TGATCTGCCAGCGCCATTCGCTACGTCCGAAC Chr1
##   [2]  2269 GCCAGGTAGAATTCCAAAAATTGAAAGGCGTT...TGATCTGCCAGCGCCATTCGCTACGTCCGAAC Chr1
##   [3]  1871 CCTCCAAGTCTTAAAGCCTGCCTGCGACATTC...AAACCGCGAAATCTTATCGACCATTGTGTCTC Chr1
##   [4]   283 GCCAGGTAGAATTCCAAAAATTGAAAGGCGTT...GCGGTCACTTTGCTAGGGCTGGTAGCCATGCA Chr1
##   [5]   129 GCCAGGTAGAATTCCAAAAATTGAAAGGCGTT...TATCACACCTCTCTCTACTCGAATGACACGGC Chr1
##   ...   ... ...
## [438]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGA...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ChrM
## [439]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGA...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ChrM
## [440]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGA...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ChrM
## [441]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGA...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ChrM
## [442]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGA...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ChrM

extractTranscriptSeqs

Sequences composed of several ranges, such as transcripts (or CDSs) with several exons, can be parsed with extractTranscriptSeqs. Note: the following expects the genome sequence in a file with path data/test and a valid txdb defining the ranges for that genome.

library(GenomicFeatures); library(Biostrings); library(Rsamtools)
txdb <- loadDb("./data/TAIR10.sqlite")
indexFa("data/test") # Creates index for genome fasta
## [1] "data/test.fai"
txseq <- extractTranscriptSeqs(FaFile("data/test"), txdb, use.names=TRUE) 
txseq
## DNAStringSet object of length 28:
##      width seq                                                                  names               
##  [1]  1688 GCCAGGTAGAATTCCAAAAATTGAAAGGCGTTT...TGATCTGCCAGCGCCATTCGCTACGTCCGAAC AT1G01010.1
##  [2]  1623 CGCTTGTCCGTTCAAGAAAAGAATCTATTCCCC...CCACGCGTAGATTAGTTGGGGACGTAAGGACC AT1G01020.1
##  [3]  1085 CGCTTGTCCGTTCAAGAAAAGAATCTATTCCCC...GAGCAGGTTCGTTTTTAGAGACTCTATACCGT AT1G01020.2
##  [4]  1905 GTCCAGCCTACCCTATCGATACCTGTCACCATT...CAAAATCCGGCAGCATTTATCATGATCACGTA AT1G01030.1
##  [5]  1239 CGCGGGAGAAACATTCGTCGAGGAATGATGCTG...ATGGTTTGCTCGACGAAGCTTATTCTGTCCGC AT2G01008.1
##  ...   ... ...
## [24]  1062 GCCGTCCACAATGACGGTGTGGATCTGAGAACG...TGGGTATCAGCCTACCGGCTGTAGGGAAATAG ATCG00020.1
## [25]    72 CCCGAGACTAGCCGCACCACGTGAGCTTTGGTG...AAACATCCTGTGCACGTTTTTTTGCCGGCGAT ATCG00030.1
## [26]   324 CTGGGTTGGAATGGGGCAAATACTTTCAAGGAG...CAGGTGCAAATTGGCATCGCCGCACGGACGTT ATMG00030.1
## [27]   462 GGATCTGACGCCACCAGAACCACAACTCCTGGA...GTGCTATGATGGCGTATGAGCGGGGGGTCGTA ATMG00010.1
## [28]  2568 CCGCTGATAGAGGTCAAATTACCAACAATAGTG...TGATAAATGGGATGCTTTAAGGACTCGTCGGC ATMG00020.1

Homework 6

See here.

Session Info

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] GenomicFeatures_1.46.1      AnnotationDbi_1.56.2        rtracklayer_1.54.0         
##  [4] systemPipeR_2.0.4           ShortRead_1.52.0            GenomicAlignments_1.30.0   
##  [7] SummarizedExperiment_1.24.0 Biobase_2.54.0              MatrixGenerics_1.6.0       
## [10] matrixStats_0.61.0          Rsamtools_2.10.0            GenomicRanges_1.46.1       
## [13] BiocParallel_1.28.2         ggseqlogo_0.1               ggplot2_3.3.5              
## [16] seqLogo_1.60.0              Biostrings_2.62.0           GenomeInfoDb_1.30.0        
## [19] XVector_0.34.0              IRanges_2.28.0              S4Vectors_0.32.3           
## [22] BiocGenerics_0.40.0         BiocStyle_2.22.0           
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7           bit64_4.0.5            filelock_1.0.2         progress_1.2.2        
##  [5] RColorBrewer_1.1-2     httr_1.4.2             tools_4.1.3            bslib_0.3.1           
##  [9] utf8_1.2.2             R6_2.5.1               DBI_1.1.1              colorspace_2.0-2      
## [13] withr_2.4.3            prettyunits_1.1.1      tidyselect_1.1.1       curl_4.3.2            
## [17] bit_4.0.4              compiler_4.1.3         cli_3.1.0              xml2_1.3.3            
## [21] DelayedArray_0.20.0    labeling_0.4.2         bookdown_0.24          sass_0.4.0            
## [25] scales_1.1.1           rappdirs_0.3.3         stringr_1.4.0          digest_0.6.29         
## [29] rmarkdown_2.13         jpeg_0.1-9             pkgconfig_2.0.3        htmltools_0.5.2       
## [33] highr_0.9              dbplyr_2.1.1           fastmap_1.1.0          htmlwidgets_1.5.4     
## [37] rlang_1.0.2            RSQLite_2.2.9          farver_2.1.0           jquerylib_0.1.4       
## [41] BiocIO_1.4.0           generics_0.1.1         hwriter_1.3.2          jsonlite_1.8.0        
## [45] dplyr_1.0.7            RCurl_1.98-1.5         magrittr_2.0.2         GenomeInfoDbData_1.2.7
## [49] Matrix_1.4-0           Rcpp_1.0.8.2           munsell_0.5.0          fansi_0.5.0           
## [53] lifecycle_1.0.1        stringi_1.7.6          yaml_2.3.5             zlibbioc_1.40.0       
## [57] BiocFileCache_2.2.0    blob_1.2.2             parallel_4.1.3         crayon_1.4.2          
## [61] lattice_0.20-45        hms_1.1.1              KEGGREST_1.34.0        knitr_1.37            
## [65] pillar_1.6.4           rjson_0.2.20           codetools_0.2-18       biomaRt_2.50.1        
## [69] XML_3.99-0.8           glue_1.6.2             evaluate_0.15          blogdown_1.8.2        
## [73] latticeExtra_0.6-29    BiocManager_1.30.16    png_0.1-7              vctrs_0.3.8           
## [77] gtable_0.3.0           purrr_0.3.4            assertthat_0.2.1       cachem_1.0.6          
## [81] xfun_0.30              restfulr_0.0.13        viridisLite_0.4.0      tibble_3.1.6          
## [85] memoise_2.0.1          ellipsis_0.3.2

References

Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nat. Methods 12 (2): 115–21. https://doi.org/10.1038/nmeth.3252.

Lawrence, Michael, Wolfgang Huber, Hervé Pagès, Patrick Aboyoun, Marc Carlson, Robert Gentleman, Martin T Morgan, and Vincent J Carey. 2013. “Software for Computing and Annotating Genomic Ranges.” PLoS Comput. Biol. 9 (8): e1003118. https://doi.org/10.1371/journal.pcbi.1003118.

Last modified 2022-05-15: some edits (2f93ebfb)