<p>You can get specific sequences efficiently using the following packages:</p>

<pre><code>library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data &lt;- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges &lt;- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %&gt;% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %&gt;% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %&gt;% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %&gt;% pull(strand) # the strand your sequence is located on
)
&gt; getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT
</code></pre>

<p>To write the sequences to a file, use <code>writeXStringSet()</code></p>

<p>For the second part of what you need, <code>GenomicFeatures</code> has a number of functions that obtain the annotated sequences of the whole genome by its type:</p>

<pre><code>transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)
</code></pre>

<p>And this is how you use one of them to find whether your sequences overlap a particular type of DNA:</p>

<pre><code>Find_overlap &lt;- function(Ref_Data,
                         Overlap_Message){

  Ref_Data &lt;- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data &lt;- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %&gt;%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %&gt;%
    as.data.table()

  Methylation_data[Ref_Data, 
                   Overlap := Overlap_Message, 
                   on = .(start &gt;= start, 
                          end &lt;=  end)] # wherever the range overlaps, assign an identifier to it
}


Methylated_in_5UTR &lt;- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %&gt;% 
  Find_overlap(., "5' UTR") %&gt;% 
  print
</code></pre>

<pre><code>seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2020-05-17T06:45:22.653Z
2020-05-17T06:45:22.653Z
Anurag N. Sharma
<p>I have a problem in R. I have the following dataSet (the first three rows shown) (the 5th number is the methylation level in its region):</p>

<pre><code>chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +
</code></pre>

<p>and then I'm trying to get the sequences from every row in the dataSet (427.000 rows in total):</p>

<pre><code>library(BSgenome.Hsapiens.UCSC.hg19)
genome &lt;- BSgenome.Hsapiens.UCSC.hg19
chr&lt;-as.matrix(as.character(mydataSet[,1]))
start&lt;-as.matrix(as.integer(as.character(mydataSet[,2]))-50)
end&lt;-as.matrix(as.integer(as.character(mydataSet[,3]))+50)
Sequences_for_50_plus_bases&lt;-getSeq(genome,chr,start=start,end=end)
</code></pre>

<p>I'd like to find a new fuction where I can have this dataSet, and find for its rows (which represent sequences in the human genome) if they refer to a promoter, or an enhancer or any other DNA regions' information, because I'd like to compare the differences between the regions' methylation levels.</p>

<p>Is there a way to do this?</p>

<r><bioconductor><phylogenetics><methylation><hg19>
Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19)
# Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19)

(Asked by: username_0 on 2020-05-16)

I have a problem in R. I have the following dataSet (the first three rows shown) (the 5th number is the methylation level in its region):
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +

and then I'm trying to get the sequences from every row in the dataSet (427.000 rows in total):
library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet[,1]))
start<-as.matrix(as.integer(as.character(mydataSet[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet[,3]))+50)
Sequences_for_50_plus_bases<-getSeq(genome,chr,start=start,end=end)

I'd like to find a new fuction where I can have this dataSet, and find for its rows (which represent sequences in the human genome) if they refer to a promoter, or an enhancer or any other DNA regions' information, because I'd like to compare the differences between the regions' methylation levels.
Is there a way to do this?

**Tags:** r, bioconductor, phylogenetics, methylation, hg19

---

## 1 answers

### Answer by username_1 on 2020-05-17. Score: 1 | Accepted answer

You can get specific sequences efficiently using the following packages:
library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data <- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges <- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %>% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %>% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %>% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %>% pull(strand) # the strand your sequence is located on
)
> getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT

To write the sequences to a file, use writeXStringSet()
For the second part of what you need, GenomicFeatures has a number of functions that obtain the annotated sequences of the whole genome by its type:
transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)

And this is how you use one of them to find whether your sequences overlap a particular type of DNA:
Find_overlap <- function(Ref_Data,
                         Overlap_Message){

  Ref_Data <- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data <- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %>%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %>%
    as.data.table()

  Methylation_data[Ref_Data, 
                   Overlap := Overlap_Message, 
                   on = .(start >= start, 
                          end <=  end)] # wherever the range overlaps, assign an identifier to it
}


Methylated_in_5UTR <- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %>% 
  Find_overlap(., "5' UTR") %>% 
  print

seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
```


---

