<p>It seems like the <code>seq_start</code> and <code>seq_end</code> are sensitive to the direction. Try swapping the two.</p>

CC BY-SA 4.0
2024-02-20T18:39:31.733Z
2024-02-20T19:12:46.153Z
2024-02-20T19:12:46.153Z
RaphiS
<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>

CC BY-SA 4.0
2024-02-20T22:23:11.523Z
2024-02-20T22:23:11.523Z
darked89
<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>

CC BY-SA 4.0
2021-03-28T17:04:17.203Z
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 ?
pippo1980
<python><sequence-annotation><blastn><entrez><efetch>
Finding annotated counterpart after BLASTn with efetch (Biopython)
# 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.


---

