count number of GC for a given genomic ranges

I have Granges with start and end but I want to count number of gc in stretch of DNA. however, i'm not able to properly subset bsgenome to be able to do that. here's my approach.

> gr_pro
GRanges object with 3 ranges and 2 metadata columns:
      seqnames              ranges strand |      symbol                    key
         <Rle>           <IRanges>  <Rle> | <character>            <character>
  [1]     chr6 143510692-143512691      - |       FUCA2 FUCA2_ENSG0000000103..
  [2]     chr6   53615972-53617971      - |        GCLC GCLC_ENSG00000001084..
  [3]     chr6   41071944-41073943      + |        NFYA NFYA_ENSG00000001167..
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

##command used
lapply(gr_pro, function(x) Biostrings::countPattern("CG", BSgenome.Hsapiens.UCSC.hg38::Hsapiens[[x]]))

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'subject' in selecting a method for function 'countPattern': error in evaluating the argument 'i' in selecting a method for function '[[': GRanges objects don't support [[, as.list(), lapply(), or unlist() at
  the moment

pls see example data

> dput(head(gene_promoters_encode1kb_chr_subset[,c("symbol", "key")], 3))
new("GRanges", seqnames = new("Rle", values = structure(6L, .Label = c("chr1", 
"chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", 
"chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", 
"chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", 
"chrY"), class = "factor"), lengths = 3L, elementMetadata = NULL, 
    metadata = list()), ranges = new("IRanges", start = c(143510692L, 
53615972L, 41071944L), width = c(2000L, 2000L, 2000L), NAMES = NULL, 
    elementType = "ANY", elementMetadata = NULL, metadata = list()), 
    strand = new("Rle", values = structure(2:1, .Label = c("+", 
    "-", "*"), class = "factor"), lengths = 2:1, elementMetadata = NULL, 
        metadata = list()), seqinfo = new("Seqinfo", seqnames = c("chr1", 
    "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", 
    "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", 
    "chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", 
    "chrY"), seqlengths = c(NA_integer_, NA_integer_, NA_integer_, 
    NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
    NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
    NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
    NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
    NA_integer_), is_circular = c(NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA), genome = c(NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_)), elementMetadata = new("DFrame", rownames = NULL, 
        nrows = 3L, listData = list(symbol = c("FUCA2", "GCLC", 
        "NFYA"), key = c("FUCA2_ENSG00000001036.13", "GCLC_ENSG00000001084.12", 
        "NFYA_ENSG00000001167.14")), elementType = "ANY", elementMetadata = NULL, 
        metadata = list()), elementType = "ANY", metadata = list())

bsgenome biostrings epigenetics

• 2,039 views • link

updated 3.0 years ago by ATpoint

• written 3.0 years ago by Alewa

2 answers

There might be a better option, but here is something that combines the getSet() method from BSgenome to retrieve the DNA sequence from a BSgenome object based on a GRanges object, and then uses alphabetFrequency from Biostrings to get the number of ATCGs per GRanges entry, returning GC content as a slot in the original GRanges object:

library(Biostrings)
library(BSgenome.Ecoli.NCBI.20080805)
library(GenomicRanges)

gr <- GRanges(seqnames = rep(seqnames(BSgenome.Ecoli.NCBI.20080805)[1], 3),
              ranges = IRanges(start = c(300, 30000, 100000),
                               end = c(1000, 50000, 200000)))

getGC <- function(genome, ranges){

  require(BSgenome)
  require(Biostrings)
  require(GenomicRanges)

  s <- BSgenome::getSeq(genome, ranges)
  a <- Biostrings::alphabetFrequency(s)

  gc <- apply(a[,c("G","C")], 1, sum)
  al <- apply(a[,c("A", "T", "C", "G")], 1, sum)
  ranges$gc.content <- 100*gc/al
  ranges

}

getGC(BSgenome.Ecoli.NCBI.20080805, gr)

GRanges object with 3 ranges and 1 metadata column:
       seqnames        ranges strand | gc.content
          <Rle>     <IRanges>  <Rle> |  <numeric>
  [1] NC_008253      300-1000      * |    53.9230
  [2] NC_008253   30000-50000      * |    52.4374
  [3] NC_008253 100000-200000      * |    51.0375
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

No idea how this scales for GRanges objects with a large number of rows, probably poorly.

• 0 views • link

written 3.0 years ago by ATpoint

Just in case the OP was asking for occurrences of the "GC" sequence, instead of GC content, and using the same example as shown by ATpoint, you could do:

library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomicRanges)

# Example GR
gr <- GRanges(seqnames = rep(seqnames(BSgenome.Hsapiens.UCSC.hg38)[1], 3),
              ranges = IRanges(start = c(300, 30000, 100000),
                               end = c(1000, 50000, 200000)))

# Get genome sequences
seqs <- BSgenome::getSeq(BSgenome.Hsapiens.UCSC.hg38, gr)

# Count the pattern at each sequence
Biostrings::vcountPattern("GC",seqs)

# stringr alternative to count the pattern:
sapply(seqs,function(x){
  stringr::str_count(as.character(x),"GC")
})

• 0 views • link

written 3.0 years ago by Papyrus

Indeed, I might have misunderstood the question.

• 0 views • link

written 3.0 years ago by ATpoint