Hi everyone,
I'd like to extract intronic sequences for my study species using R. For that, I was planning to use the GenomicFeatures package. I got to the point when I have a GRanges object that I produced using intronicParts() which, I suppose, contains intron ranges. How should I proceed to get intron sequences in a fasta format? I did not find an answer in the manual: https://bioconductor.org/packages/release/bioc/manuals/GenomicFeatures/man/GenomicFeatures.pdf
Thank you in advance!
•
written
by
carroll •
Hi,
you could check out the getSeq function from BSgenome. It takes a GRanges object and a genome as an input and gives you the sequences. You can then save them as fasta using the write.fasta function of seqinr.
getSeq will return a Biostrings object, which you can save as a fasta using Biostrings::writeXStringSet, if you want to skip using seqinr.
I don't know what to do about the error.
As a workaround, you can also do something like this:
#For every intron
for(i in 1:nrow(Intronlist)){
#Get the strand the intron is located at
strand <- Intronlist$Strand[i]
#Get chromosome the intron is located at
chromosome <- Intronlist$Chromosome.scaffold.name[i]
#Get Intron border coordinates
upstream_end <- Intronlist$intron_upstream_end_cooordinate[i]
downstream_end <- Intronlist$intron_downstream_end_cooordinate[i]
#Get the respective genomic sequence (for introns of the minus strand, the upstrteam
test2 <- Views(DNAstringSet_of_Genome[[as.character(chromosome)]], start=upstream_end , end=downstream_end )
#For introns of the minus strand, the coordinates have to be switched, since the VIEWS function only handels end coordinates, which are higher than the start coordinates
if(strand == -1) test2 <- Views((DNAstringSet_of_Genome[[as.character(chromosome)]], start=downstream_end, end=upstream_end)
#Get the sequence from the Views object
test <- as.character(test2)
sequence_range <- as.character(test2)
test3 <- data.frame(test2@ranges)
#Views extracts the sequence on the plus strand, so for minus strand introns, replace the intron sequence with the reverse complement
sequence_range[strand=="-1"] <- gsub("A","t",sequence_range[strand=="-1"])
sequence_range[strand=="-1"] <- gsub("C","g",sequence_range[strand=="-1"])
sequence_range[strand=="-1"] <- gsub("T","a",sequence_range[strand=="-1"])
sequence_range[strand=="-1"] <- gsub("G","c",sequence_range[strand=="-1"])
sequence_range[strand=="-1"] <- toupper(sequence_range[strand=="-1"])
sequence_range[strand=="-1"] <- reverse(sequence_range[strand=="-1"])
#Store the intron sequence
Intronlist$intron_seq[i] <- sequence_range
}
Log in to answer this question.
