# Bi278 Lab 6 v2 rpsblast and Orthofinder
### By Lee Ferenc 10/25/2023
## 1. Reversed position specific BLAST
Rpsblast is used to compare a protein query sequence against a database of protein profiles (NCBI has a collection of protein profiles on their conserved domain database (CDD)).
COGs (Clusters Of Orthologous Groups) were developed as a way of defining curated groups of proteins across genomes that should be orthologous to each other. This is to extropolate well-studied groups onto poorly-studied ones.
### Step 1: Establish BLAST onto ssh.
```
export BLASTDB="/courses/bi278/Course_Materials/blastdb"
echo $BLASTDB
```
### Step 2: Run rpsblast on your group of proteins.
You need to start with a protein fasta file.
```
rpsblast -query bbqs395_hypothetical.faa -db Cog -evalue 0.01 -max_target_seqs 25 -outfmt "6 qseqid sseqid evalue qcovs stitle" > bbqs395_hypothetical.rpsblast.table
```
#### My example:
```
rpsblast -query subset_hypo_proteins_prokka_09262023.faa -db Cog -evalue 0.01 -max_target_seqs 25 -outfmt "6 qseqid sseqid evalue qcovs stitle" > hypo_proteins_prokka.rpsblast.table
```
### Step 3: Keep only the best COG match for each query sequence that have query coverage of at least 70.
This means at least 70% of the query sequence is covered by the subject that matched.
```
awk -F'[\t,]' '!x[$1]++ && $4>=70 {print $1,$5}' OFS="\t" bbqs395_hypothetical.rpsblast.table > bbqs395_hypothetical.cog.table
```
The command:
* Takes the format 5 table and splits it by tabs and commas. We have to do this to get the COG #
* Match the same query ($1) as long as each query has a coverage of 70 or above
* Keeps best match and prints two columns (=fields)
#### My example:
`awk -F'[\t,]' '!x[$1]++ && $4>=70 {print $1,$5}' OFS="\t" hypo_proteins_prokka.rpsblast.table > hypo_proteins_prokka.cog.table
`
Okay so that didn't work so none of them were =>70. Oops. Trying w/ another subset of just toxin proteins.
```
rpsblast -query toxin.faa -db Cog -evalue 0.01 -max_target_seqs 25 -outfmt "6 qseqid sseqid evalue qcovs stitle" > toxin_prokka.rpsblast.table
awk -F'[\t,]' '!x[$1]++ && $4>=70 {print $1,$5}' OFS="\t" toxin_prokka.rpsblast.table > toxin_prokka.cog.table`
```
Output:
```
BHQS69_01563 COG2867
BHQS69_01564 COG2914
BHQS69_01744 COG2337
BHQS69_01799 COG2274
BHQS69_01800 COG0845
BHQS69_02317 COG5450
```
### Step 4: Look up the functional categories for the COG numbers you have.
Info in table: `cognames2003-2014.tab`
```
awk -F "\t" 'NR==FNR {a[$1]=$2;next} {if ($2 in a){print $1, $2, a[$2]} else {print $0}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/cognames2003-2014.tab bbqs395_hypothetical.cog.table > temp
```
Command:
* Takes `cognames2003-2014.tab` and makes a look-up table that can access categories by COG number
* Takes input file and goes thru row by row and finds what category each COG number belongs to
* Writes the row along with the COG category that it looked up and stores it in a new file (`temp`)
#### My example:
```
awk -F "\t" 'NR==FNR {a[$1]=$2;next} {if ($2 in a){print $1, $2, a[$2]} else {print $0}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/cognames2003-2014.tab toxin_prokka.cog.table > temp
```
Output of file head:
```
BHQS69_01563 COG2867 J
BHQS69_01564 COG2914 V
BHQS69_01744 COG2337 V
BHQS69_01799 COG2274 V
BHQS69_01800 COG0845 MV
BHQS69_02317 COG5450 K
```
### Step 5: Keep the first category for each COG number.
COGs can belong to multiple functional categories but the first one is its primary category.
```
awk -F "\t" '{if ( length($3)>1 ) { $3 = substr($3, 0, 1) } else { $3 = $3 }; print}' OFS="\t" temp > temp2
```
#### My example (same as above)
Output of file head:
```
BHQS69_01563 COG2867 J
BHQS69_01564 COG2914 V
BHQS69_01744 COG2337 V
BHQS69_01799 COG2274 V
BHQS69_01800 COG0845 M
BHQS69_02317 COG5450 K
```
### Step 6: Get full description of the functional category
Other table for this: `fun2003-2014.tab`
```
awk -F "\t" 'NR==FNR {a[$1]=$2;next} {if ($3 in a){print $0, a[$3]} else {print $0}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/fun2003-2014.tab temp2 > bbqs395_hypothetical.cog.categorized
```
#### My example:
```
awk -F "\t" 'NR==FNR {a[$1]=$2;next} {if ($3 in a){print $0, a[$3]} else {print $0}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/fun2003-2014.tab temp2 > prokka_toxins.cog.categorized
```
Output of file head:
```
BHQS69_01563 COG2867 J Translation, ribosomal structure and biogenesis
BHQS69_01564 COG2914 V Defense mechanisms
BHQS69_01744 COG2337 V Defense mechanisms
BHQS69_01799 COG2274 V Defense mechanisms
BHQS69_01800 COG0845 M Cell wall/membrane/envelope biogenesis
BHQS69_02317 COG5450 K Transcription
```
### Step 7: Clean up (remove temp files)
`rm temp temp2`
### (Optional) Step 8: Add what the individual COG numbers stand for to the final table
```
awk -F "\t" 'NR==FNR {a[$1]=$3;next} {if ($2 in a){print $1, $2, a[$2], $3, $4} else {print "no match"}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/cognames2003-2014.tab bbqs395_hypothetical.cog.categorized
```
#### My example:
```
awk -F "\t" 'NR==FNR {a[$1]=$3;next} {if ($2 in a){print $1, $2, a[$2], $3, $4} else {print "no match"}}' OFS="\t" /courses/bi278/Course_Materials/blastdb/cognames2003-2014.tab prokka_toxins.cog.categorized
```
Output:
```
BHQS69_01563 COG2867 Ribosome association toxin PasT (RatA) of the RatAB toxin-antitoxin module J Translation, ribosomal structure and biogenesis
BHQS69_01564 COG2914 Putative antitoxin component PasI (RatB) of the RatAB toxin-antitoxin module, ubiquitin-RnfH superfamily V Defense mechanisms
BHQS69_01744 COG2337 mRNA-degrading endonuclease, toxin component of the MazEF toxin-antitoxin module V Defense mechanisms
BHQS69_01799 COG2274 ABC-type bacteriocin/lantibiotic exporters, contain an N-terminal double-glycine peptidase domain V Defense mechanisms
BHQS69_01800 COG0845 Multidrug efflux pump subunit AcrA (membrane-fusion protein) M Cell wall/membrane/envelope biogenesis
BHQS69_02317 COG5450 Transcription regulator of the Arc/MetJ class K Transcription
```
## 2. Orthofinder
Orthofinder is a popular orthogroup detection software that finds homologous groups of genes (=orthogroup).
### Step 1: Collect predicted protein sequences (PROKKA*.faa equivalent) from the genomes you want to compare into a directory. Then run Orthofinder.
```
orthofinder -f ./ -X
```
`-f` Species where protein fasta files, (`./` means here) and `-X` species I don't want orthofiner to add species names to gene ids.
#### My example:
Due to the files taking too long I used the files from `/courses/bi278/Course_Materials/lab_06/Results_Oct23/` and didn't copy since there's a lot of folders
### Step 2: Check the overall statistics.
```
cat Comparative_Genomics_Statistics/Statistics_Overall.tsv
```
Looks something like:
```
Number of species 7
Number of genes 26880
Number of genes in orthogroups 25526
Number of unassigned genes 1354
Percentage of genes in orthogroups 95.0
Percentage of unassigned genes 5.0
Number of orthogroups 4357
Number of species-specific orthogroups 272
Number of genes in species-specific orthogroups 960
Percentage of genes in species-specific orthogroups 3.6
Mean orthogroup size 5.9
Median orthogroup size 7.0
```
Note: For all examples I'm using the same files so I don't see a reason to create examples when the path is the only different thing
### Step 3: Examine species tree file.
Orthofinder uses phylogenies (trees) to more accurately find orthogroups so we need to know what the species tree looks like.
```
cat Species_Tree/SpeciesTree_rooted_node_labels.txt
```
Output:
```
(Burkholderia_pseudomallei_K96243_132:0.139475,(((bbqs433_prokka:0.00032017,bbqs395_prokka:0.00012987)N4:0.00129303,bbqs859.prokka:0.0026247)N2:0.0949503,((bhqs11.prokka:0.00262492,bhqs155_prokka:0.000209001)N5:0.000404395,bhqs69_prokka:0.000737986)N3:0.10081)N1:0.139475)N0;
```
I input it into the tree of life (see lab 3) and it looked the same as the tutorial.
### Step 4: Examine file.
```
head Orthogroups/Orthogroups.tsv
```
Output: Incomprehensible cat for `.tsv` format but it's all the genes from each genome that belond to each detected orthogrou (first column). Listed largest to smallest.
You can just take the `.tsv` file and upload it to Google Slides so it's readable. Then you can see how many groups of size 3+. It was 453.
### Step 4: Examine File (But Different)
```
awk '$2=="0" && $6=="0" && $7=="0" && $8=="0" && $9>=3' ../Orthogroups/Orthogroups.GeneCount.tsv | wc -l
```
You also get 453. :)
### Step 5: Find fasta files for the genes by orthogroup id for each orthogroup:
Found here:
`Orthogroup_Sequences`
Remember, each fasta file is one gene.
### Step 6: Check out which genes are unique to individual genomes
Here:
```
Orthogroups/Orthogroups_UnassignedGenes.tsv
```
## 3. Make and use a custom BLAST database
### Step 1: Make your own directory to store BLAST databases in your bi278 home.
```
mkdir enfere24blastDB
```
### Step 2: Add this directory to your $BLASTDB
```
export BLASTDB="/courses/bi278/Course_Materials/blastdb:/home2/enfere24/enfere24blastDB"
```
### Step 3: Make protein databases
```
makeblastdb -in PATH/PROKKA.faa -dbtype prot -title NAME_prot -out enfere24blastDB/NAME_prot -parse_seqids
```
This will make a series of files: *.phr, *.pin, *.psq and more.
#### My Example:
```
makeblastdb -in lab_03/PROKKA_09262023.faa -dbtype prot -title BHQS69_prot -out enfere24blastDB/BHQS69_prot -parse_seqids
```
### Step 4: Make nucleotide databases
```
makeblastdb -in PATH/PROKKA.ffn -dbtype nucl -title NAME_nucl -out enfere24blastDB/NAME_nucl -parse_seqids
```
Files created: *.nhr, *.nin, *.nsq and more.
My example:
```
makeblastdb -in lab_03/PROKKA_09262023.ffn -dbtype nucl -title BHQS69_nucl -out enfere24blastDB/BHQS69_nucl -parse_seqids
```
(Don't worry for the project we'll make a different BLAST database this just testing. I did have to create a new DB cause I accidentally ran this. )
### Step 5: Combine BLAST databases
You can combine any BLAST databases of the same type (nucleotide or protein). This way you can search against multiple databases. Lab example:
```
blastdb_aliastool -dbtype prot -title 4burk_prot -out myblastdb/4burk_prot -dblist "bb395_prot bb433_prot bh155_prot bh69_prot"
```
Make sure to use the exact titles you give the databases you want to combine in `-dblist`.
### Step 6: Test.
```
blastp -query /courses/bi278/Course_Materials/lab_06/tssH_bpseudo.fasta -db YOUR_DB -outfmt 6
```
My example:
```
blastp -query /courses/bi278/Course_Materials/lab_06/tssH_bpseudo.fasta -db BHQS69_prot -outfmt 6
```
So rememer the database name is "BHQS69_prot" or whatever you named it in step 3
Output:
```
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01420 59.216 510 198 3 1 501 1 509 3.00e-166 509
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02560 48.356 517 256 4 1 515 53 560 1.22e-144 452
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02560 49.062 320 150 2 613 929 582 891 5.17e-78 271
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03640 44.372 462 249 2 1 454 15 476 5.89e-108 354
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03640 49.235 327 161 2 604 925 540 866 1.85e-86 295
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01886 41.935 496 250 7 10 495 3 470 2.20e-101 335
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01886 39.244 344 197 5 610 948 520 856 5.70e-73 256
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02364 42.995 414 206 10 32 430 21 419 3.34e-79 271
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02364 35.197 304 179 8 658 956 459 749 5.80e-46 174
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02351 72.222 18 5 0 984 1001 1311 1328 0.035 34.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02351 72.222 18 5 0 987 1004 1311 1328 0.035 34.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02351 72.222 18 5 0 990 1007 1311 1328 0.035 34.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02351 68.750 16 5 0 993 1008 1311 1326 1.5 29.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02351 61.111 18 7 0 981 998 1311 1328 1.8 29.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03064 35.484 62 37 2 736 795 137 197 0.11 32.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02258 28.846 52 33 1 303 350 125 176 0.12 32.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_00225 39.474 38 23 0 688 725 275 312 0.20 32.0
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03498 26.882 186 103 8 687 857 26 193 0.28 31.2
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02301 30.081 123 52 8 662 772 156 256 0.33 31.2
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02301 27.731 119 57 6 208 311 165 269 0.39 31.2
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01188 27.612 134 84 3 462 594 77 198 0.58 30.4
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_00631 30.612 49 31 1 676 724 149 194 2.1 28.5
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01307 40.476 42 24 1 693 733 31 72 3.4 27.7
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03134 38.776 49 29 1 572 619 20 68 4.2 26.9
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_03397 36.207 58 31 2 877 928 322 379 4.6 27.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_01126 45.455 22 12 0 949 970 27 48 4.9 25.0
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02698 35.526 76 37 2 25 91 29 101 5.4 26.6
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_02758 30.526 95 40 4 4 81 87 172 5.7 26.9
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE BHQS69_00111 25.974 77 52 1 225 301 101 172 6.7 26.9
```
As for NCBI questions, I genuinely don't have any right now. If I do, I'll bug my cousin and he'll probably have the answer or know who will. Only good family member. (He works for the NDB which is apart of the NIH and is specifically involved with coordinating with other NIH branches)
Sign in with Wallet
Connect another wallet