I need to identify the sequence of a gene in the complete genome sequence . I thought it was simple, instead it is not a straightforward task !

My method was the following: I downloaded the FASTA format file of the complete genome sequence of the archaea organism that I am studying, that is CP003685, and then I downloaded the CDS file from NCBI by clicking on ">Send to>Coding Sequencing>FASTA Nucleotide>Create File". In this way I have a txt file (always FASTA format file) in which I find a sequence for each gene (if I am not wrong).

So for example, I should find, in the complete genome sequence, the sequences of the CDS file (the last that I downloaded), right?

For example, if I consider:

>lcl|CP003685.1_cds_AFN03983.1_1007 [locus_tag=PFC_05180] [protein=DNA protection protein DPS] [protein_id=AFN03983.1] [location=complement(940023..940580)] [gbkey=CDS]

I should find in the complete genome sequence, this sequence, right ? But this is not. Why?

I suspect that I miss very important passages.

Thank you in advance.

  • $\begingroup$ See How to identify genes from a genome assembly of C. Elegans?. Exonerate and wise2 should be able to do this. $\endgroup$
    – terdon
    Jun 9 at 10:58
  • $\begingroup$ I simply copy the sequence of the gene (like the one I have written in the post) and in the file with the complete genome sequence (CP003685) I click ctrl+F and paste the reverse complement of the sequence of the gene and search it. But it does not find it. Sorry if I did not really specified this very crude way to do it :') . For the introns I have no idea... $\endgroup$
    – Manuela
    Jun 9 at 11:48
  • $\begingroup$ It is a non-professional way but I do not understand why it does not work since the two sequences should be the same... the two files I am considering (as said below) are only two different ways to represent the genome DNA sequence, right ? $\endgroup$
    – Manuela
    Jun 9 at 11:56
  • 3
    $\begingroup$ Well yes, kinda, but you have line breaks. You will never find the entire thing unless you convert it to a one-line format first. Try searching for a smaller section, say ~10 nucleotides and if the first one doesn't work try some others. You need to find a sequence that appears in a single line otherwise you won't find a match. $\endgroup$
    – terdon
    Jun 9 at 12:11
  • $\begingroup$ Thank you very much @terdon. I did not know this and now it seems to be finally clear putting all pieces together ! $\endgroup$
    – Manuela
    Jun 9 at 12:31

I would just blast it. When blasting locally, you need to first make a database from your genome, so assuming you got the command-line version of blast installed you can do something like

makeblastdb -in my_study_genome.fa -dbtype nucl
blastn -max_target_seqs 10 -db my_study_genome.fa -evalue 1e-10 -outfmt 6 -query my_downloaded_gene_of_interest.fasta -out gene_in_my_study_genome.blast

blast search for anything that is "similar enough" and also considers both strands, so for example, if there is an SNP in the middle of the gene. Hence even if the searched gene and the reference genome are from different sequencing runs / individuals, blast will still find it. Blast will also handle linebreaks for you (in fasta files the nucleotide sequence is quite often wrapped by 60 or 100 characters a line). Long story short, there is not much that can go wrong with blast (while with full text search you can easily shoot yourself in foot).

-- edit 2 --

As @terdon suggested, potentially if you would locating a transcript sequence (excluding introns), you would do better with software like exonerate or genewise, because they model exon/intron boundaries. However, even in those cases blast will give you an approximate position of the gene.

  • $\begingroup$ @Manuela this really should be done using a tool that models genes (try exonerate or genewise). Blast is great but will not be able to model introns correctly if your species has splicing. $\endgroup$
    – terdon
    Jun 9 at 12:38

In your sequence header you will see:


It's a reverse strand gene. Reverse complement the sequence, then search the genome sequence.

  • $\begingroup$ Thank you very much for the answer but it is still not working . I do not understand why. I used this software to make the reverse complement, copied and searched in the complete genome sequence but I still do not find it. $\endgroup$
    – Manuela
    Jun 9 at 10:22

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.