---
# System prepended metadata

title: Lab 6

---

# Lab 6

`mkdir labwk6` 
## Exercise 1: Reversed position specific BLAST

1. `export BLASTDB="/courses/bi278/Course_Materials/blastdb"` to establish where the system should look for the BLAST database 
2. `echo $BLASTDB` got the output /courses/bi278/Course_Materials/blastdb so I know it worked
3. Next I needed to run rpsblast on my group of proteins, starting with a fasta file `rpsblast -query hypPROKKA -db Cog -evalue 0.01 -max_target_seqs 25 -outfmt "6 qseqid sseqid evalue qcovs stitle" > hyp.rpsblast.table`
4. NOTE: Keep only the best COG match for each query sequence that have query coverage of at least 70, meaning at least 70% of the query sequence is covered by the subject that matched. Use the commad `awk -F'[\t,]' '!x[$1]++ && $4>=70 {print $1,$5}' OFS="\t" hyp.rpsblast.table > hyp.cog.table`
5. Next we can look up what each of these COGs is and categorize them to the 17 broad functional categroies the individuals COGs belong to. First we want to look up the functional categories for the COG numbers we have. This is stored in `cognames2003-2014.tab`
Use the command: `awk -F "\t" 'NR==FNR {a[$1]=$2;next} {if ($2 in a){print $1, $2, a[$2]} else {print $0}}' OFS="t\t" /courses/bi278/Course_Materials/blastdb/cognames2003-2014.tab hyp.cog.table > temp` 
Then `head temp` to make sure this worked and to see what was created 
6. I now want to keep the first category for each COG number, because the first one is the primary category `awk -F "\t" '{if ( length($3)>1 ) { $3 = substr($3, 0, 1) } else { $3 = $3 }; print}' OFS="\t" temp > temp2`
I then opened temp2 to look at the new file
7. Now I wanted to have a full description of the functional category which can be looked up in the table `fun2003-2014.tab` 
Use the command `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 > hyp.cog.categorized` 
8. `rm temp temp2` 

## Exercise 2: Make and use a custom BLAST database

9. `cd ~` to make sure the directory for the BLAST database is in my bi278 home`mkdir blastdb` 
10. In order to search against the new database I need to add this to the $BLASTDB using the command `export BLASTDB="/courses/bi278/Course_Materials/blastdb:/home2/slpuro26/blastdb"` 
11. To make protein databases, use the command `makeblastdb -in /courses/bi278/Course_Materials/labwk4/prokka.faa -dbtype prot -title genomename_prot -out blastdb/genomename_prot - parse_seqids` 
BUT modify on an amino acid multi-fasta file. The command now looks like`makeblastdb -in ~/labwk5/hypPROKKA -dbtype prot -title hypPROKKA_protdb -out blastdb/hypPROKKA_protdb -parse_seqids`
12. To make a nucleotide database, use the command `makeblastdb -in ~/labwk5/PROKKA_09162022.ffn -dbtype nucl -title PROKKA_nucldb -out blastdb/PROKKA_nucldb -parse_seqids`
13. Lastly, combine any BLAST databases of the same type to search against multiple datbases at once. To do this, list the exact titles given to the databases I want to combine in `-dblist`. 
command example `blastdb_aliastool -dbtype prot -title 4burk_prot -out myblastdb/4burk_port -dblist "bb395_prot bb433_prot bh155_prot bh69_prot "` I could do this using the files in the lab_06 directory 

14. to check that it worked and searched against the new databases I created `blastp -query tssH_bpseudo.fasta -db hypPROKKA_protdb -outfmt 6 ` this worked! my results: 
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_03207	29.927	137	81	4	687	811	26	159	0.32	29.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02418	34.146	41	25	1	200	240	437	475	0.39	29.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02877	38.776	49	29	1	572	619	20	68	1.3	26.6
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02020	34.694	49	32	0	501	549	565	613	1.8	27.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02020	56.667	30	13	0	954	983	610	639	5.9	25.4
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02020	28.947	76	52	1	507	580	619	694	7.1	25.0
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02020	45.238	42	23	0	955	996	629	670	8.0	25.0
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_01993	46.512	43	23	0	77	119	104	146	2.5	26.6
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_01456	45.000	40	20	1	423	462	36	73	3.3	26.2
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_00274	44.444	27	15	0	566	592	43	69	3.5	24.3
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_03341	28.125	32	23	0	511	542	369	400	5.3	25.4
tr|A0A8A4DWX9|A0A8A4DWX9_BURPE	BO395_02667	23.256	86	50	2	431	512	44	117	8.3	24.3

## Exercise 3: Orthofinder and its output

15. Collect all the predicted protein sequences into a directory `mkdir forortho` and then copied the four .faa files from the lab_06 file into this new forortho file `cp bhqs69_prokka.faa bbqs395_prokka.faa bbqs433_prokka.faa bhqs155_prokka.faa ~/labwk6/forortho`
16. then run orthofinder within that folder`#orthofinder -f ./ -X`
17. `mv OrthoFinder/Results_Oct25/ .`
18. `rmdir OrthoFinder`
19. `cd Results_Oct25`

These commands are not actually being run

20. Now `cd /courses/bi278/Course_Materials/lab_06/Results_Oct22` 
21. `cat Comparative_Genomics_Statistics/Statistics_Overall.tsv`
22. To look at the species tree file: `cat Species_Tree/SpeciesTree_rooted_node_labels_txt` or can visualize using iToL as we have in a previous lab

24. To look at the file that describe the orthogroups of genes shared by more than one genome `head Phylogenetic_Hierarchical_Orthogroups/N0.tsv` 
This lists out all of the genes that belong to each of the detected orthogroups, losted from smallest to largest. This file makes it easy to sort out the members or orthogroups that are only found in a paticular subset of genomes
25. For example, to find genes that are shared between genome 1 and genomes 4 and 5. Use the command `awk -F "\t" '$4!="" && $5=="" && $6=="" && $7!="" && $8!="" {print}' OFS="\t" Phylogenetic_Hierarchical_Orthogroups/N0.tsv
26. To see which genes are unqiue to individual genomes `head Orthogroups/Orthogroups_UnassignedGenes.tsv`
my results for this were 
Orthogroup	Burkholderia_pseudomallei_K96243_132	bbqs395_prokka	bbqs433_prokka	bhqs155_prokka	bhqs69_prokka
OG0004292	BPSL0020				
OG0004293	BPSL0024				
OG0004294	BPSL0025				
OG0004295	BPSL0035				
OG0004296	BPSL0042				
OG0004297	BPSL0043				
OG0004298	BPSL0056				
OG0004299	BPSL0057				
OG0004300	BPSL0059