0
$\begingroup$

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!

$\endgroup$
1
  • $\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 '20 at 5:07
0
$\begingroup$

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"
$\endgroup$
6
  • $\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 '20 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 '20 at 18:53
  • $\begingroup$ adj should be GenomicRanges object...if it was clear from the solution above $\endgroup$
    – StupidWolf
    Jun 30 '20 at 18:53
  • $\begingroup$ thank you I will try $\endgroup$
    – Joanne
    Jun 30 '20 at 18:59
  • $\begingroup$ brilliant, worked like a charm! thank you $\endgroup$
    – Joanne
    Jul 3 '20 at 17:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.