What kind of BLAST do I need to do to accomplish this task?

Question

(I am a complete newbie to bioinformatics, so please bear with me!).

I recently used BLAST to compare a 500 bp nucleotide sequence in S. cerevisiae to a bunch of other 500 bp nucleotide sequences in various fungi (these other sequences are homologs of the S. cerevisiae gene). Ideally, when I click "BLAST", I want to get nothing back. What I'm trying to double check is the assumption that my S. cerevisiae sequence is NOT similar to my other input sequences.

Here's what I ran on the computer, using the BLAST website:

After I run this query, I only get one sequence back:

which tells me that the sequence I got as an output needs to be excluded from my dataset, as it's too similar to the S. cerevisiae sequence.

The problem is, I need to run this type of query and analysis about 6000 more times for 6000 different S. cerevisiae genes and their homologs. I don't know how to do this programmatically. The Biopython packages seem to only compare one sequence against the whole NCBI database, and I was really struggling to understand how to do this via the command line (was struggling with the package installations a bit, and it also seemed to only be for one sequence input).

Is what I'm doing here a multiple sequence alignment? When I looked into programmatically doing an MSA, the output seemed to only be the alignment, but I want to see an output like the website gives me. I'm not really sure what I should be looking into.

Any guidance would help!

Answer 1

Clarification

Just making sure - you have ~6000 datasets - one for each S.cerevisiae gene of interest - which are made up of homologues of each respective gene? And you want to filter out genes that are too similar?

In that case, it would be helpful to know how you're defining a sequence as 'similar enough' for it to be filtered. Given they're all homologues, it's likely there's going to be a high baseline of sequence similarity. In your example, you consider 97.40% identity across the whole query as different enough to be filtered out, so I'm getting the feeling you specifically want to filter out identical sequences?

Also, when you talk about genes here, do you want to filter out gene products (i.e. proteins) that are too similar to the actual coding sequences (CDS)? The CDS of two homologues might be X% different, but that difference might not be as high once translated into proteins due to how codons work. Proteins can go down to levels of similarity as low as ~30% identity and still have the same function, so depending on the context of your work (taxonomic, structural, functional) you might want to translate these CDS into protein, so further clarification would be helpful!

The actual bioinformatics question

To answer your question of programming, if the context of this comparison is always X S.cerevisiae gene vs. Y dataset of homologues, a MSA is definitely not the answer (you would do that if you wanted to see how all of the homologues compare to each other, rather than to specifically X gene). Also, using the web-client for NCBI BLAST will take a lot of time, so that's also not recommended for the kind of throughput you need.

Pairwise sequence alignments that you do on your local machine/server between X genes and all of its respective Y dataset sequences is probably the best way to go. I would suggest looking into needleall - it allows you to compare X sequences to Y sequences in a pairwise way, though it performs a global alignment rather than local like BLAST, so it might be slower.

I don't know what your computational resources are, so in case that might be too computationally heavy, the faster alternative would be to install the BLAST+ suite, which allows you to run BLAST queries like you did on the web browser, except on your own machine. Specifically, you could make your own custom individual BLAST databases out of each Y dataset using makeblastdb from the suite. Then you use blastn (or blastp if you decide to translate your sequences) on each dataset with its respective S.cerevisae sequence.

BLAST+ also allows you to customise the actual output as well, making it easier to parse your results programatically for when you want to remove sequences that are too similar. There are many examples online of how to do this, but have a go at it on your own first with the BLAST+ documentation I linked :)