<p>Say your file looks like this:</p>
<pre><code>cat test.fa
&gt;chr01
ACCGGTT
&gt;chr02
ATGCATGC
</code></pre>
<p>Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:</p>
<pre><code>library(GenomicRanges)
library(BSgenome)
library(Biostrings)

adj &lt;- data.frame(chr=c(&quot;chr01&quot;,&quot;chr02&quot;),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta &lt;- readDNAStringSet(&quot;test.fa&quot;)
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

as.character(SEQ)
[1] &quot;C&quot; &quot;G&quot;
</code></pre>
<p>or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:</p>
<pre><code>library(Rsamtools)
fasta = FaFile(&quot;test.fa&quot;)
indexFa('test.fa')

as.character(getSeq(FaFile(&quot;test.fa&quot;),adj))
chr01 chr02 
  &quot;C&quot;   &quot;G&quot;
</code></pre>

CC BY-SA 4.0
2020-06-30T18:50:57.290Z
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?
Joanne
CC BY-SA 4.0
2020-06-30T18:53:17.670Z
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
StupidWolf
CC BY-SA 4.0
2020-06-30T18:53:40.973Z
```adj``` should be GenomicRanges object...if it was clear from the solution above
StupidWolf
CC BY-SA 4.0
2020-06-30T18:59:49.420Z
thank you I will try
Joanne
CC BY-SA 4.0
2020-07-03T17:03:40.673Z
brilliant, worked like a charm! thank you
Joanne
CC BY-SA 4.0
2020-07-04T22:33:33.620Z
you're welcome :) glad it worked for you
StupidWolf
CC BY-SA 4.0
2020-06-30T08:53:24.433Z
2020-06-30T08:53:24.433Z
StupidWolf
<p>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 <code>seq.name</code> for the target chromosome and I've added an <code>adjacent_pos</code> variable with the location on the target chromosome that I want. Using <code>read.fasta</code> from the phylotools package gives me a genome in the format of a dataframe with <code>seq.name</code> and <code>seq.text</code> for each chromosome and unplaced scaffold.</p>
<p>Working on a test data set with a single chromosome, I successfully used something like this:</p>
<pre><code>adj &lt;- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base &lt;- str_sub(genome$seq.text, adj)
</code></pre>
<p>Now back to the whole genome, I can't wrap my head around how to direct <code>str_sub</code> to the correct <code>seq.name</code> for each <code>seq.position</code>. I'm sure there's an elegant solution, but I haven't been able to find it!</p>

CC BY-SA 4.0
2020-06-30T05:07:14.820Z
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.
thomas duge de bernonville
<r>
how to extract a different substring for each row of a dataframe in R?
# how to extract a different substring for each row of a dataframe in R?

(Asked by: username_0 on 2020-06-29)

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!

**Tags:** r
> **Comments:**
> **username_1**: 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.  

---

## 1 answers

### Answer by username_2 on 2020-06-30. Score: 1 | Accepted answer

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"

> **Comments:**
> **username_0**: 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?  
> **username_2**: 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  
> **username_2**: ```adj``` should be GenomicRanges object...if it was clear from the solution above  
> **username_0**: thank you I will try  
> **username_0**: brilliant, worked like a charm! thank you  
> **username_2**: you're welcome :) glad it worked for you  

---

