0
$\begingroup$

I have a very, very long list of genome coordinates of areas. I need to gather a sequence in order to compare it to another database of sequences that I have. How do I produce the genome sequences from these coordinates in R?

I am very new to this and need a detailed answer.

$\endgroup$
  • 1
    $\begingroup$ What genome? Namely, is it one for which there's already a BSgenome package? $\endgroup$ – Devon Ryan Dec 1 '17 at 22:06
  • $\begingroup$ The genome was given to me, it is not on there. My problem, more specifically I suppose, is that I was given a couple files. The genome was run through RepeatMasker to obtain coordinates of transposable elements. So, I have the coordinates from the map file that was produced from RepeatMasker and I would like to identify the sequences of these coordinates from my genome (which is too big to view). After I obtain the sequences, I will need to compare those sequences to a library of known elements that was also given to me in another file. So, I am lost especially since I am very new to this $\endgroup$ – Noob Dec 1 '17 at 22:36
  • 1
    $\begingroup$ OK, does this really have to be done in R? This particular task is a bit better done on the command line than in R. $\endgroup$ – Devon Ryan Dec 1 '17 at 22:39
  • $\begingroup$ Do you know what type of commands this would consist of? $\endgroup$ – Noob Dec 27 '17 at 3:30
2
$\begingroup$

This doesn't use R, but a simple and quite convenient method is to use bedtools getfasta with your fasta file. You'll need to make a BED file out of the RepeatMasker output:

awk 'BEGIN{OFS="\t"}{name=sprintf("%s_%s", $10, $11); strand="+"; if($9 == "C"){strand="-"}; print $5,$6-1,$7,name,".", strand}' repeatmasker.output > rmsk.bed

Then, given that and your genome.fa fasta:

bedtools getfasta -name -s -fi genome.fa -bed rmsk.bed -fo output.fa

This will produce a multiline fasta file names like AluSx_SINE/Alu (due to the -name option) where the sequences are reverse complemented if appropriate.

You can install bedtools either via conda (conda install -c bioconda bedtools) or github.

$\endgroup$

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.