We changed our privacy policy. Read more.

Hot answers tagged

11

There are three possible problems that come to mind. Masking Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it. Score/e-value thresholds Another source of complication is that even if a match is found, that match will have very bad scores. Both ...


9

I'm trying to figure out why. This will be a longer read, Tl;dr at the end: Doing a match There is a good match of Q against LC074724 in 20170824: $ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \ -task blastn -max_hsps 1 | grep LC074724 Q LC074724.1 2581 Note that the bitscore of the LC074724 match is 2581. But there is ...


9

The E-value and the mapping qualities are two very different things. The E-value is "a parameter that describes the number of hits one can 'expect' to see by chance when searching a database of a particular size". More details can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect The mapping ...


9

I was able to reproduce your problem, even in a blank environment. $ conda create -n test $ conda activate test $ conda install bioconda::blast=2.7.1 Solving environment: failed PackagesNotFoundError: The following packages are not available from current channels: - bioconda::blast=2.7.1 - boost=1.64 - bioconda::blast=2.7.1 - gnutls - bioconda::...


8

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work: >K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15 So does removing everything after the 1st space: >K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 There is nothing in the FASTA definition that would ...


8

I believe you're looking for env_nr? It's listed as such, under Metagenomic proteins in the blastp webpage. It appears that the word within brackets should be supplied alongside the -db parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession. Edit: I followed up on this with a little more digging ...


8

-O is to specify the name of the output file. You want the -P option. So, from your example wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/


7

BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species. Since you are new to bioinformatics I recommend you search an ortholog database that already exists. Here are some databases: orthoDB OMA EggNog All of these should ...


7

The Blast E is the expected frequency of obtaining false positive. This will depend on your query size (number of nucleotide or amino acid residues) and the size of the database. With a short query and a large database you are more likely to have a sequence in your database that matches your query by simple chance. In summary, the smaller the query and the ...


6

With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of grep ^'>' fasta_file | wc -l ? This would show how many contigs you have in the file.


6

You can extract fasta sequence from a blastdb constructed from a fasta file using blastdbcmd which should be installed when you install blast/makeblastdb. blastdbcmd -entry all -db <database label> -out <outfile> If you had a database called my_database which contained the files: my_database.nhr my_database.nsq my_database.nin and you ...


6

You can specify send to stdout using out='-' in the Biopython wrapper. from Bio.Blast.Applications import NcbiblastnCommandline import pandas as pd cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-') output = cline()[0].strip() rows = [line.split() for line in output.splitlines()] cols = ['qseqid', 'sseqid', 'pident', '...


6

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,...


5

Try changing you evalue threshold to something like 0.001. If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of 1, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the evalue ...


5

Sure, so you start out with what is called a bitscore, which is a normalized to the score calculated from the alignment between the 2 seqs which depends on the following equation. It is independent of database size: (lambda * S - ln(k))/(ln)2 Then, the p-value of a local blast is just: 1/2^bitscore So if your bitscore is 15, you need 1/32768 ...


5

There are three processes wherein light is emitted. Bioluminescence: a chemical reaction releases light. The enzyme that does this is luciferase, while its substrate is luciferin, a small molecule. There are half a dozen organisms that have analogous (evolved independently) enzymes and totally different substrates, such as fireflies, Alivibrio bacteria, ...


4

You probably want to include query start (qstart) and query end (qend) in your blast output. Something like this: blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue' In R you can take the "qstart:qend" from each line for density plot. There are many ways in R to plot the densities of these start and ...


4

Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA


4

As I understand, the software tool Lambda is a viable, yet lesser known alternative to BLAST in the context of taxonomic classification of NGS data.


4

The answer is it depends. 40% pident with ~100% qcovs (assuming length of protein is not very small) tells you that the virulence factor (VF) and your prokka predicted protein are likely homologs, they belong to the same (or closely related) protein family. For some research studies, you might annotate your prokka proteins as 'putative' VFs. On the other ...


4

The Refseq team and also the NCBI resource coordinators team publish a new paper every few years, so check out the many papers (e.g. here or here), but to answer your 2nd question, non-redundancy here is (I think) defined very strictly as proteins that are identical in terms of sequence and length, so the clustering is trivial, without the need for a ...


4

You could take a look at kraken2 It should fit your purpose.


4

Update 2: I looked into this a little more, with the various data sources. This is related in part to the answer submitted by OP juanjo75es, in addition to discussion on chat. I don't entirely understand the logic, but the general thrust seems to be that SPAdes makes weird assemblies somehow. Some notes that I made: REFERENCE ASSEMBLIES FIV sequence U11820....


4

The Diamond protein aligner is a very popular tool for pairwise protein alignment which is much faster than blastp while having comparable sensitivity. Diamond also supports the use of using multiple compute nodes (similar to mpiBLAST). You can view the details of the new Diamond 2.0 release in this recent Nature Methods paper (source of the image as well). ...


3

I am the author of the question / answer here Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter? To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting. That is sadly not the case, though thats what I would say is ...


3

There might be a "cleaner" way to do this, but what I would do is first remove all of the spaces from the fasta sequence names and replace them with underscores using something like sed. sed 's/ /_/g' Run your BLAST after you remove the spaces and you will get the whole ID string in the result instead of only that first part (e.g. NP_567483.1). From here ...


3

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk: $ awk 'NR==FNR{a[$1]=0; next} { if($2 in a){ a[$2]=1 } }END{ for(i in a){ print i,a[i] } }' pangenome.txt blast....


3

Biopython has a parser for Blast's default text format (which I assume is the one you have). Just use the blast-text string to specify this format in SearchIO's parse method. from Bio import SearchIO results = SearchIO.parse('myfile.txt', 'blast-text') for entry in results: # extract the information you want query = entry.id for hit in entry....


3

The blast_formatter command allows you to produce output in any of the formats supported by BLAST. However, it requires that you run blastx with --outfmt 11 (ASN.1 format). I'm not sure there are any tools that will do the conversion from the full default BLAST report. This could potentially be scripted, but would be complex enough that you're probably ...


3

The evalue associated with a blast HSP (High-scoring Segment Pair) is a measure of how often you would expect to find such an HSP in a database of this size by pure chance. For instance, if you blast the sequence ACTG against the entirety of the NCBI's nr database, you would expect several hundred thousand hits by pure chance. That is, several hundred ...


Only top voted, non community-wiki answers of a minimum length are eligible