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.

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.