Answers.Body: <p>It seems like the <code>seq_start</code> and <code>seq_end</code> are sensitive to the direction. Try swapping the two.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-20T18:39:31.733Z
Answers.LastActivityDate: 2024-02-20T19:12:46.153Z
Answers.LastEditDate: 2024-02-20T19:12:46.153Z
Answers.UserDisplayName: RaphiS
Answers.Body: <p>Lets assume you can get the <code>blastn</code> output in a <code>BED</code> format:</p>
  <pre><code>genomeA_contig001    hit_start    hit_end    query_transcript    score   strand
  genomeB_contig021    hit_start    hit_end    query_transcript    score   strand
  </code></pre>
  <p>which will contain all the unknown hits to all genomes.</p>
  <p>All you need is a set of <code>GTF</code> files specific for <code>genomeA</code>, <code>genomeB</code> and so on where you add the suitable prefix (<code>genomeA_</code> for <code>GTF</code> describing genes/exons in that genome and so on). For transcript hits you may start with just the annotated exons. Next run:</p>
  <p><code>bedtools intersect -wao -a all.transctipts.hits.bed -b all_annotated_exons.gtf my_result.tsv</code></p>
  <p>to get all the overlaps with number of overlapping bases.
  Filter and sort the <code>my_result.tsv</code> possibly importing it into a data frame using <code>pandas</code> or <code>polars</code>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-20T22:23:11.523Z
Answers.LastActivityDate: 2024-02-20T22:23:11.523Z
Answers.UserDisplayName: darked89
Body: <p>I am creating a pipeline for the identification of unknown transcripts. After a local BLASTn search of the transcripts, I have a large list of the respective hits with different genomes. I have the accession code of the respective genome, start/end coordinates, and the strand where the hit is found. Now, I would like to find if the genome where the hit was found is annotated within the coordinates returned by BLAST.</p>
  <p>I have found the possibility of using efetch to search for each genome via Biopython, but this only returns the whole annotation, even though I limit the search to the specific coordinates with:</p>
  <pre><code>Entrez.efetch(db=&quot;nuccore&quot;, id=&quot;CP054847.1&quot;, strand=&quot;2&quot;, seq_start=&quot;326053&quot;, seq_end=&quot;326786&quot;, retmode=&quot;xml&quot;) 
  </code></pre>
  <p>Of course, I could parse the XML and look for possible annotations within the given range, but not only would this be quite computationally intense since I expect around 10.000+ hits to be searched, but it also seems to return the same results for the plus/minus strand, even though the annotations should be strand-specific.</p>
  <p>I have also intended to use the standalone tool of efetch, but it would be easier to integrate the search to the pipeline using Biopython.</p>
  <p>Do you have any ideas why <code>Entrez.efetch</code> is returning the whole <code>xml</code> or is there any fast way of filtering the output? Do you have maybe another alternative/idea for this step? I would appreciate any suggestion!</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-03-28T17:04:17.203Z
Comments.Text: can you find any difference between https://www.ncbi.nlm.nih.gov/nuccore/CP054847.1/ and https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP054847.1 database entries ?
Comments.UserDisplayName: pippo1980
Tags: <python><sequence-annotation><blastn><entrez><efetch>
Title: Finding annotated counterpart after BLASTn with efetch (Biopython)
ThreadText: # Finding annotated counterpart after BLASTn with efetch (Biopython)
  
  (Asked by: username_0 on 2021-03-24)
  
  I am creating a pipeline for the identification of unknown transcripts. After a local BLASTn search of the transcripts, I have a large list of the respective hits with different genomes. I have the accession code of the respective genome, start/end coordinates, and the strand where the hit is found. Now, I would like to find if the genome where the hit was found is annotated within the coordinates returned by BLAST.
  
  I have found the possibility of using efetch to search for each genome via Biopython, but this only returns the whole annotation, even though I limit the search to the specific coordinates with:
  
  ```
  Entrez.efetch(db="nuccore", id="CP054847.1", strand="2", seq_start="326053", seq_end="326786", retmode="xml") 
  
  ```
  
  Of course, I could parse the XML and look for possible annotations within the given range, but not only would this be quite computationally intense since I expect around 10.000+ hits to be searched, but it also seems to return the same results for the plus/minus strand, even though the annotations should be strand-specific.
  
  I have also intended to use the standalone tool of efetch, but it would be easier to integrate the search to the pipeline using Biopython.
  
  Do you have any ideas why `Entrez.efetch` is returning the whole `xml` or is there any fast way of filtering the output? Do you have maybe another alternative/idea for this step? I would appreciate any suggestion!
  
  **Tags:** python, sequence-annotation, blastn, entrez, efetch
  
  > **Comments:**
  > **username_1**: can you find any difference between https://www.ncbi.nlm.nih.gov/nuccore/CP054847.1/ and https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP054847.1 database entries ?  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2024-02-20. Score: 1
  
  It seems like the `seq_start` and `seq_end` are sensitive to the direction. Try swapping the two.
  
  
  ---
  
  ### Answer by username_3 on 2024-02-20. Score: 1
  
  Lets assume you can get the `blastn` output in a `BED` format:
  
  ```
  genomeA_contig001    hit_start    hit_end    query_transcript    score   strand
  genomeB_contig021    hit_start    hit_end    query_transcript    score   strand
  
  ```
  
  which will contain all the unknown hits to all genomes.
  
  All you need is a set of `GTF` files specific for `genomeA`, `genomeB` and so on where you add the suitable prefix (`genomeA_` for `GTF` describing genes/exons in that genome and so on). For transcript hits you may start with just the annotated exons. Next run:
  
  `bedtools intersect -wao -a all.transctipts.hits.bed -b all_annotated_exons.gtf my_result.tsv`
  
  to get all the overlaps with number of overlapping bases.
  Filter and sort the `my_result.tsv` possibly importing it into a data frame using `pandas` or `polars`.
  
  
  ---
  
  