3
$\begingroup$

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!

$\endgroup$
1

2 Answers 2

1
$\begingroup$

It seems like the seq_start and seq_end are sensitive to the direction. Try swapping the two.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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