I have a dataframe with putative alignments of thousands of probe sequences. I would like to pull out the adjacent nucleotide from the genome for each alignment. My dataframe includes seq.name for the target chromosome and I've added an adjacent_pos variable with the location on the target chromosome that I want. Using read.fasta from the phylotools package gives me a genome in the format of a dataframe with seq.name and seq.text for each chromosome and unplaced scaffold.

Working on a test data set with a single chromosome, I successfully used something like this:

adj <- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base <- str_sub(genome$seq.text, adj)

Now back to the whole genome, I can't wrap my head around how to direct str_sub to the correct seq.name for each seq.position. I'm sure there's an elegant solution, but I haven't been able to find it!

  • $\begingroup$ please can you provide a small reproducible example for your alignments and genome objects? the merge() function from the data.table package will surely help you. $\endgroup$ Jun 30, 2020 at 5:07

1 Answer 1


Say your file looks like this:

cat test.fa

Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:


adj <- data.frame(chr=c("chr01","chr02"),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta <- readDNAStringSet("test.fa")
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

[1] "C" "G"

or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:

fasta = FaFile("test.fa")

chr01 chr02 
  "C"   "G"
  • $\begingroup$ but it's not just one location from each chromosome. I have about 160k alignments, each needing a different location pulled out from a 3G genome that has 21 chromosomes plus over 7500 unplaced scaffolds. Still perhaps the GenomicRanges object is the way to go? $\endgroup$
    – Joanne
    Jun 30, 2020 at 18:50
  • $\begingroup$ yes.. the code above works for as many locations as you have. You just need to construct the genomic ranges object like I did above, and I suggest using the Rsamtools solution if your file is huge $\endgroup$
    – StupidWolf
    Jun 30, 2020 at 18:53
  • $\begingroup$ adj should be GenomicRanges object...if it was clear from the solution above $\endgroup$
    – StupidWolf
    Jun 30, 2020 at 18:53
  • $\begingroup$ thank you I will try $\endgroup$
    – Joanne
    Jun 30, 2020 at 18:59
  • $\begingroup$ brilliant, worked like a charm! thank you $\endgroup$
    – Joanne
    Jul 3, 2020 at 17:03

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.