how to extract a different substring for each row of a dataframe in R?

Question

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!

Answer 1

Say your file looks like this:

cat test.fa
>chr01
ACCGGTT
>chr02
ATGCATGC

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

library(GenomicRanges)
library(BSgenome)
library(Biostrings)

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

as.character(SEQ)
[1] "C" "G"

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

library(Rsamtools)
fasta = FaFile("test.fa")
indexFa('test.fa')

as.character(getSeq(FaFile("test.fa"),adj))
chr01 chr02 
  "C"   "G"