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