0
$\begingroup$

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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.