---
# System prepended metadata

title: Bi278 LAB 05 | 10/17/2023

---

**`` Command-line BLAST``**

++**Notes**++
___
Important task names for BLAST+
**Program: 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)

**Program: 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)

___


++**Lab START**++
```
echo $BLASTDB 
> 
```
++**Notes**++
This is an environmental variable. If the output is nothing, that means the environmental variable is currently empty. You need this before you can run BLAST during a new ssh session. 

```
export BLASTDB="$BLASTDB:/courses/bi278/Course_Materials/blastdb"
echo $BLASTDB
> :/courses/bi278/Course_Materials/blastdb
```


++**Notes**++
blastn -help
blastp -help
These requests can be made for a codec of resources on properly using BLAST
Alternatively, you may look here: https://www.ncbi.nlm.nih.gov/books/NBK279684/

```
cp -r /home2/asyero26/lab_04/bbqs395/*.faa ./


```

++**Notes**++
Amino acid/Protein sequence
> PROKKA_09162022.faa 

Nucleotide bases
> bbqs395_fastaseq

```
awk 'BEGIN {RS=">"} /hypothetical protein/ {print ">"$0}' ORS="" PROKKA_09162022.faa 
```
++**Notes**++
This allows you to delimit (seperate) records(sequences) in multi-fasta files
**Here's what each part does**
* Read the file you designate at the end of the command using whatever you said was the record separator (RS=">")
* Match a particular pattern (/pattern/) in each record, such as the phrase ‘hypothetical protein’
* 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="")


```
SOMECOMMAND | head

```
++**Notes**++
This allows you to preview an output before actually running the request. It's also useful for combining cmds where the left hand is the output you'd like to analyze 


**'awk' utilization**
```
awk 'BEGIN {RS=">"} /replication/ {print ">"$0}' ORS="" PROKKA_09162022.faa | grep "replication" |wc -l 
```
++**Notes**++
* performed a count and stratification for "replication"


```
awk 'BEGIN {RS=">"; srand()} (rand()<=.002) {print ">"$0}' ORS="" PROKKA_09162022.faa > subset#1

>BO395_00143 cysteine synthase B
MTLRLAEPIARHSIASDHPVQLATPLLEVHPGLLLKFESDNPGGSHKVRAARFIVRSAIA
SGDIVPGHTTVIEKTGGNFGFGLALACQAHNVPVELAVGLEFSAVKRRCLELFGARLIGL
DMLERGATPREVVEWHLEQANELERHYFYTDQFHNPASVAAHELETGPEIAMQLKAWPEV
KKLTFVSCAGTGASLTGIATALKHAGYAVDVVLVEPQGCDSRNGIFTSHKLEGMSVGVVP
PFLDWNLINEVRTVSHEEAMETRQTLARTLGFIAGNTAAACLTVAHTLARSASCEHKVLS
LMYDHGLWYLR
>BO395_00895 type III secretion system inner membrane export apparatus protein BpscS2 (T3SS-2)
MNYDALTHLTTQALTLCLLVSLPAVAIAACSGLAIAFLQAITSLQDASIGQSVKLVVVTL
VVVLSAPWGAAQIQNFAHVLLQVMFA
>BO395_01018 dehalogenase
MSATTALVPKAIIFDAYGTLFDVHSVVAAAEQMFPGQGDALSQLWRHKQIEYTHLRTLAD
PAGTHYQPFWNITLDALRFAAHTLGLSLGHQGQKRLMDEYACLSAFPDALPALRQLRDTL
PRHATGDSPRLALLSNGNPQMLDIAIKSAGMTGLFDHVLSADAVRAYQPAPAVYALGTQA
FGGVPAREIVFVSSNGWNAAGAAWFGYTTFWINRQHASPEELGVTFHGTGQGMHDLSPFL
QNFTFSSRSAGRSRSAATQARDSSQRQP
>BO395_01822 Cyclohexadienyl dehydrogenase
MIRVRGLRAARAFWVMDVAGFSFNKLVVAGVGLIGGSLARALRERDEIGGARQVIGIGRS
TDSTARALALGVIDQAVALTDEAALREALAGADIVLLAAPVAQTQPLLERIAPWLDAATI
ITDAGSTKSDAVAAAHAALGSRAGQFVPGHPIAGREASGVDAALPDLYVGRNVVLCPMAG
NTPEALERIASMWRAAGATVHEMSSVEHDRVFAAVSHLPHVLSFALVEQILASPDAALKF
SFAAGGFRDFTRIAASSPEMWRDVCVANRVALLEELDAYTAVLTRLRAAINTADGAELEA
VFARSRTARMAWQPGGQPVRLDDTAA
>BO395_03102 Membrane-bound lytic murein transglycosylase F
MKPGSDSLAQVARRQRIWPKAVLACVLCVAAQAQASAVVMMIGGMKNTAGQAASAQPVAA
PEPRSSVVTLYNPSGRVSTTFVPGMLMPVTTAPSGTVASRVMALTPLIAEVSHAANIDSA
LLMAVIDVESGGNPQAVSPKGATGLMQLMPATGKRHGAADLFDPRQNIMAGARYLHILMQ
QFGNIELALAAYNAGEGAVKKYGMQIPPYDETMAYVPKVLGRYQHYRTRTAAAATPGVPQ
```
++**Notes**++
This is how you subset data effectively. 
(rand()<=.002) = Asks for a subset of .002% of the inputfile

* After running this commend, I ended up with 5 proteins.


__++**BLAST**++__

```
blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -query subset#1
```
++**Notes**++
In the place of "refseq_select_prot" can be "refseq_select_prot" or any other database. 
swissprot = smaller (> 0.4 million sequences) manually curated database, best information if match is found [match unlikely]

refseq_select_prot = larger (> 63 million sequences). Plan to use this most of the 

```
-outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10

-outfmt 3 -num_descriptions 10 -num_alignments 10

```
++**Notes**++
There are various options that allow you to moderate the output of your sequences. 
The ``-max_target_seqs,`` ``-num_descriptions``, and ``-num_alignments`` options limit the results you see somewhere in the middle of the database search process. When you run BLAST for a project, DON’T USE THESE OPTIONS because you will get different results if these numbers are too small. Use either the defaults (250 or 500 depending on which option) and filter your results after the BLAST search, or at least use a larger number (e.g. 100+).
* Output fmt 6 is preferable, it outputs with the following columns "qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore"

```
blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -query subset#1 -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 > 09162022_blastp_refs_outf6

```


++**Notes**++
``1=="BO395_01177" ``= 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").

``blastdbcmd`` = find the entries from my search target database that match the subject sequence ids



++**Nucleotide Analysis**

```
blastn -task megablast -num_threads 2 -query bbqs395_fastaseq -db env_nt -evalue 1e-6 > bbqs395_blastn_env_nt.faa

blastn -task megablast -num_threads 2 -query bbqs395_fastaseq -db nt_prok -evalue 1e-6 > bbqs395_blastn_nt_prok.faa
```
++**Notes**++
Two options to replace "env_nt" with(nucleotide database):

**nt_prok** = Contains prokaryotic nucleotide sequences included in the regular nt (nucleotide) database

**env_nt** = Contains nucleotide sequences from environmental samples. Usually used for "invisible" organisms i.e COVID 

* Important Note: We're using blastn instead of blastp

***BLAST+ order for Profesor Noh***
qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs stitle sscinames staxid

* Important because BLAST+ has the unfortunate characteristic of truncating subject names

```
awk -F "\t" '{print $x,$y,$z}' OFS="\t" example_blastp_outfmt6.txt
```
++**Notes**++
You can think of "$x" as the column number. This allows you to output/print certain column numbers epending on what you're loking for

++** POST-ish LAB QUESTION **++

___
___
___
++**Q1**++
**blastn**
matrix = BLOSUM62 - Change depending on whether you're looking for matches between closely related sequences (higher) or distant (lower)
evalue = 10 (1000 for blastn-short) - Change when interested in different significance of matches
max_target_seqs = 500 - Change when interested in different # of aligned sequences to keep
num_descriptions = 500 - CHange when interested in more or less detailed descriptions 
num_alignments = 250 - Change when interested in diff # of db sequences to show alignments for (wider or smaller search)

**blastp**
matrix = BLOSUM62, PAM30 for blastp-short, 
evalue = 10
max_target_seqs = 500
num_descriptions = 500
num_alignments = 250

Thse filters would be applied after creating subsets, whle running the blast.



++**Q2**++

```
awk 'BEGIN {RS=">"} /BO395_03165 transporter protein/ {print ">"$0}' ORS="" PROKKA_09162022.faa > transporter_protein
```
```
blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -query transporter_protein -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 > transporter_protein_blastp_refs_outf6

blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 6e-6 -query transporter_protein -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 10 > transporter_protein_blastp_refs_outf6_v2
# changed eval to 6e-6

blastp -task blastp-fast -num_threads 2 -db refseq_select_prot -evalue 1e-6 -query transporter_protein -outfmt "6 std ppos qcovs stitle sscinames staxid" -max_target_seqs 3 > transporter_protein_blastp_refs_outf6_v3
# changed max_target_seqs to 4
 ```
 
 **Oservations:**
 V1 & V2 = These two pretty much had the same output. This likely indicates that there is very little chance that the outputs from V1 were due to chance/error 
 
 V3 = There were three matches outputted. They were the first three matches from V1 and V2.
 "Warning: [blastp] Examining 5 or more matches is recommended"
 
 General Observation: I noticed that they all take about the same time to finish.
 