---
# System prepended metadata

title: 'Lab 5: BLAST'

---

# Lab 5: BLAST

Blast+ options;
**BlastP**:
blastp- traditional blastp to compare a protein query to a protein database
blastp-short:	blastp optimized for queries shorter than 30 residues
blastp-fast:	a faster version of blastp that uses a larger word size (6 vs. 2-3)
**Blastn**:
blastn:	traditional blastn requiring an exact match of 11, to compare a nucleotide query to a nucleotide database
blastn-short:	blastn program optimized for sequences shorter than 50 bases
megablast:	traditional megablast used to find very similar sequences (e.g. intraspecies or closely related species)
dc-megablast:	discontinuous megablast used to find more distinct sequences (e.g. interspecies)

## 1.1 Command Line Blast
1. To get blast to runn a search from bi278, you need to tell it where to find your databases.To do this, we need to set an environmental variable called $BLASTDB.we run
echo $BLASTDB
i got nothing which meant that the environmental variable is empty
2. Tell Blast where these databases are:
```
export BLASTDB="$BLASTDB:/courses/bi278/Course_Materials/blastdb"
```
3. to check if this worked try:
```
echo $BLASTDB
```
I did this and it gave me the path of the file: :/courses/bi278/Course_Materials/blastdb
*BTW TO SEE HOW BLAST+ WORKS I CAN TYPE*:
blastn -help
blastp -help

I looked at the prokka output from last labs, i went into lab_03, into bbqs395, then got the faa(amino acid) and ffn (nucleotide) file

## 1.1 Run BLASTp
1. These files contain multiple sequences that each belong to a different gene, name of each gene declared in specfic way using >. we can use this symbol to seperate sequences in a multi=fatsa files. we can do this by using:
```
awk 'BEGIN {RS=">"} /hypothetical protein/ {print ">"$0}' ORS="" PATH/PROKKA_FILE
```
this looked like:
```
awk 'BEGIN {RS=">"} /hypothetical protein/ {print ">"$0}' ORS="" /home2/cibene26/lab_03/bbqs395/PROKKA_10032023.faa
```

(RS=">") is a record seperator
(/pattern/) match to particular pattern
When a record contains that pattern, then print the symbol > followed by the record with the pattern ($0)
When printing, it doesn’t use an additional output record separator (ORS="")

2. We had to move all of these hypothetical proteins to a new file we do this by doing `>` this is what it looked like( new file w all hyp proteuins is after the >):
```
awk 'BEGIN {RS=">"} /hypothetical protein/ {print ">"$0}' ORS="" /home2/cibene26/lab_03/bbqs395/PROKKA_10032023.faa > BWAHPPROKKA.faa
```


The entire **awk command** is contained within single quotes, and anything to do with the output format is outside these single quotes. After the command, you specify the input file you want awk to work on. Then you can either have everything printed to your screen, or have it sent to a file (as above using the > operator). Typically before sending anything to a file I will always check what the output looks like using the | operator to pipe my command to head or tail:
```
SOMECOMMAND | head
```

The absolutely simplest way to subset your data may be to just open up your file using a text editor, and copy and paste however many you want into a new file. But here is a way to do it in command line using awk by taking a percentage (for example, 0.2% = 0.002) of a file’s records:
```
awk 'BEGIN {RS=">"; srand()} (rand()<=.002) {print ">"$0}' ORS="" INPUTFILE`
```

for me this looked like:
```
awk 'BEGIN {RS=">"; srand()} (rand()<=.002) {print ">"$0}' ORS="" /home2/cibene26/lab_03/bbqs395/PROKKA_10032023.faa`
```

2. Take a subset of your hypotheical proteins 


```
awk 'BEGIN {RS=">"; srand()} (rand()<=.002) {print ">"$0}' ORS="" BWAHPPROKKA.faa >SUBHPPROKKA.faa`
```
*instead of using path you can move things to a new file this is how:
the use of the `>` will send info to whatever file i name after it

PROKKA_10032023.faa is all proteins
SUBHPPROKKA.faa (3-5) hypothetical proteins
BWAHPPROKKA.faa is all hps seperated



**We will be working with 2 protein databases: swissprot or refseq_select_prot and 2 nucleotide databases env_nt and nt_prok.**

3.Run a protein blast (blastp). This is the basic protein (amino acid) query blast to start with and adjust as you go with optional flags:

```
blastp -task blastp-fast -num_threads 2 -db PROTEIN_DB -evalue 1e-6 -query FAA_FILE
```

for me this looked like:

lab_05]$
```
blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -query SUBHPPROKKA.faa `
```


 *The evalue threshold (-evalue 1e-6) determines the stringency of your search. If you don’t get any matches, try removing this option. But when you look at your results, note that good e-values for BLAST searches are considered to be those that are 1e-6 or smaller (including 0). What this means is you can adjust your BLAST parameters to get matches, but not all matches are informative.*
 
 I added these specifications to the command above( for how we want our search outputs to look);
 -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 
 could also do this but i didnt try it:
 -outfmt 3 -num_descriptions 10 -num_alignments 10
 **Run blastp**
 the new command would look like:
```
blastp -task blastp-fast -num_threads 2 -db
refseq_select_prot -evalue 1e-6 -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 -query INPUTFILE > OUTPUTFILE
```
 or
```
`blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -outfmt 3 -num_descriptions 10 -num_alignments 10 -query INPUTFILE > OUTPUTFILE`
```

# STOPPED HERE# in my case this looked like:
```
blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 -query SUBHPPROKKA.faa > 2SUBHPPROKKA.faa
```

2SUBHPPROKKA.faa is a reorganized blastp ran through file

3 The `-max_target_seqs`, `-num_descriptions`, and `-num_alignments` options limit the results you see somewhere in the middle of the database search process. Will take some time When BLAST is running, it’s finding the best matches using a heuristic to search through all of these potential match subjects in the database.You can use & at the end of a command to run a search in the background

Here is what the output for a different search I ran looked like. All 10 results ($2) are for the same query sequence ($1).You can also use format 6 results to get the amino acid or nucleotide sequences of the subjects you are matching to.

Let’s say I’m interested in finding all the subject sequences of one of my queries BO395_01177. The subject sequence ids that match this query can be found in my output format 6 results. I just need to get the ones that match my specific query of interest ($1=="BO395_01177").
 
 
 Table 2. Available BLAST output format specifiers

specifier	meaning
qseqid	Query Seq-id
qgi	Query GI
qacc	Query accesion
sseqid	Subject Seq-id
sallseqid	All subject Seq-id(s), separated by a ‘;’
sgi	Subject GI
sallgi	All subject GIs
sacc	Subject accession
sallacc	All subject accessions
qstart	Start of alignment in query
qend	End of alignment in query
sstart	Start of alignment in subject
send	End of alignment in subject
qseq	Aligned part of query sequence
sseq	Aligned part of subject sequence
evalue	Expect value
bitscore	Bit score
score	Raw score
length	Alignment length
pident	Percentage of identical matches
nident	Number of identical matches
mismatch	Number of mismatches
positive	Number of positive-scoring matches
gapopen	Number of gap openings
gaps	Total number of gap
ppos	Percentage of positive-scoring matches
frames	Query and subject frames separated by a ‘/’
qframe	Query frame
sframe	Subject frame
btop	Blast traceback operations (BTOP)
staxids	unique Subject Taxonomy ID(s), separated by a ‘;’(in numerical order)
sscinames	unique Subject Scientific Name(s), separated by a ‘;’
scomnames	unique Subject Common Name(s), separated by a ‘;’
sblastnames	unique Subject Blast Name(s), separated by a ‘;’ (in alphabetical order)
sskingdoms	unique Subject Super Kingdom(s), separated by a ‘;’ (in alphabetical order)
stitle	Subject Title
salltitles	All Subject Title(s), separated by a ‘<>’
sstrand	Subject Strand
qcovs	Query Coverage Per Subject (for all HSPs)
qcovhsp	Query Coverage Per HSP

## 1.2 Run Blastn(ERROR)
To call up fatsa format sequences based on their newly indexed gene ids( we are using gene ids found in subset hypothetical; protein file, my SUBHPPROKKA.faa):
samtools faidx PATH/PROKKA_*faa GENEID1 GENEID2 …
in my case this looked like:
samtools faidx ./SUBHPPROKKA.faa bbqs295_02086 bbqs295_02612 bbqs295_03216 bbqs295_03241 > NUCPROKKA.faa

gave me error, this is what i ran bc the info saved to this file: 
blastn -task megablast -num_threads 2 -query SUBHPPROKKA.faa.fai -db nt_prok -evalue 1e-6
say:
Near line 1, there's a line that doesn't look like plausible data, but it's not marked as defline or comment.

# Run Blastn 10/22
i am using my subset hypothetical protein taking those gene idsa and doing samtools on it from the ffn file and sending thsie to a new file to run blast on
 this is how:
 1. make sure i got the nucletode sequences from the gene ids pulled from the ffn file( got ATGC for each hp gene id)
```
samtools faidx /home2/cibene26/lab_03/bbqs395/PROKKA_10032023.ffn bbqs295_02086 bbqs295_02612 bbqs295_03216 bbqs295_03241
```
 2.  link it back to the ffn to get the nucleotide file, send those nucletide sequences to a file to run blastn on
 ```
 samtools faidx /home2/cibene26/lab_03/bbqs395/PROKKA_10032023.ffn bbqs295_02086 bbqs295_02612 bbqs295_03216 bbqs295_03241 > BTPROKKA.ffn
```
BTPPROKKA is my nucletode sequences from my hp gen ids

# 1.3 Parse Blast Results Output

1.editing my results, using format 3, I ran I(could sned results fo a new file the wya i wnat themmorganized)
```
blastn -task megablast -num_threads 2 -query BTPROKKA.ffn -db nt_prok -evalue 1e-6 -outfmt 3 -num_descriptions 10 -num_alignments 10
```