# BI278 Group Project
## Olivia Schirle & Kaz Hirata
## Download fna and faa files from NCBI
First, we downloaded the genomes from NCBI for each of our selected organisms. Using the following access numbers, we searched for each of the genomes in the NCBI genome database - https://www.ncbi.nlm.nih.gov/genome. At the top of the page, under "Reference genome", where it says "Dowload sequences in FASTA format for genome, protein", we selected "genome" to download the genome file (.fna) and "protein to download the protein file (.faa). We moved the files from our Desktop into our working directory (/courses/bi278/oeschi23). We then unzipped the files. Lastly, we renamed the files to more intuitive names that included the genus and species.
**Organism (access number)**
* Mycobacterium leprae TN (NC_002677)
* Mycobacterium avium 104 (NC_008595)
* Mycobacterium intracellulare subsp. intracellulare MTCC 9506 (GCF_000298095.1)
* Yersinia pestis KIM10+ (GCF_000006645.1)
* Yersinia pseudotuberculosis IP 32953 (GCF_000047365.1)
* Yersinia similis (GCF_000582515.1)
* Leptospira interrogans (GCF_022559725.1)
* Leptospira biflexa Patoc 1 (GCF_000017685.1)
**Code**
```
ssh oeschi23@bi278
mkdir group_project
cp /personal/oeschi23/BI278/BI278_Group_Project_Genome_Sequences/* ~/group_project/
gunzip GCF* # Unzip the genomefiles
# Rename the files
mv GCF_000006645.1_ASM664v1_genomic.fna ./Y.pestis.fna
mv GCF_000195855.1_ASM19585v1_genomic.fna ./M.leprae.fna
mv GCF_000014985.1_ASM1498v1_genomic.fna ./M.avium.fna
mv GCF_000298095.1_ASM29809v1_genomic.fna ./M.indicus_pranii.fna
mv GCF_000047365.1_ASM4736v1_genomic.fna ./Y.pseudotuberculosis.fna
mv GCF_000582515.1_ASM58251v1_genomic.fna ./Y.similis.fna
mv GCF_013137915.1_ASM1313791v1_genomic.fna ./F.nucleatum.fna
mv GCF_004006635.1_ASM400663v1_genomic.fna ./F.necrophorum.fna
mv GCF_022559725.1_ASM2255972v1_genomic.fna ./L.interrogans.fna
mv GCF_000017685.1_ASM1768v1_genomic.fna ./L.biflexa.fna
```
## Collect Basic Genome Stats
We used the genome files to compute the genome size and GC% for each of the organisms.
```
nano
```
We used our shell script from one of our labs to prints the genome size and number of G and C nucleotides.
Shell script:
```
#!/bin/bash
grep -v ">" $1 | tr -d -c GCATgcat | wc -c
grep -v ">" $1 | tr -d -c GCgc | wc -c
```
```
# Run the shell script on each of the genomes
sh genome_stats.sh Y.pestis.fna
```
| Organism | Pathogenicity | Genome Size (bp) | GC % |
| ------------------------------ | ------------------ | ---------------- | -------- |
| *Yersinia pestis* | High Pathogenicity | 4701745 | 0.47695 |
| *Yersinia pseudotuberculosis* | Low Pathogenicity | 4840898 | 0.47550 |
| *Yersinia similis* | Non pathogen | 4964409 | 0.46979 |
| *Mycobacterium leprae* | High Pathogenicity | 3268203 | 0.57797 |
| *Mycobacterium avium* | Low Pathogenicity | 5475491 | 0.68988 |
| *Mycobacterium indicus pranii* | Non pathogen | 5589007 | 0.68028 |
| *Fusobacterium nucleatum* | High pathogenicity | 2473944 | 0.26963 |
| *Fusobacterium necrophorum* | Low pathogenicity | 2678415 | 0.34030 |
| *Leptospira interrogans* | High pathogenicity | 5521455 | 0.35418 |
| *Leptospira biflexa* | Non pathogenic | 3951448 | 0.38895 |
### Compare Genome Stats
## BLAST
Goal: Compare similarity of pathogenic and non-pathogenic genomes
First, we created a new directry for BLAST database using `mkdir blastdb_pathogens`
Next, we used .faa files of pathogens we downloaded to create the protein BLAST databse.
The files are
```
M.avium.faa Y.pseudotuberculosis.faa F.necrophorum.faa
L.interrogans.faa Y.pestis.faa F.nucleatum.faa M.leprae.faa
```
We ran
`makeblastdb -in species_name.faa -dbtype prot -title species_name_prot -out blastdb_pathogens/species_name_prot -parse_seqids`
Run this command for all of the files above. This gives you protein BLAST database for all the pathogens.
Next, repeat the steps with .faa files of non-pathogens.
The files are
```
M.indicus_pranii.faa Y.similis.faa L.biflexa.faa
```
## Number of toxin-antitoxin modules
| Organism | Pathogenicity | Total number of toxin-antitoxin systems | Type I | type II | Number of hypothetical proteins |
| ------------------------------ | ------------------ | --------------------------------------- | ------ | ------- | ------------------------------- |
| *Yersinia pestis* | High Pathogenicity | 19 | 1 | 18 | 338 |
| *Yersinia pseudotuberculosis* | Low Pathogenicity | 15 | 1 | 12 | 367 |
| *Yersinia similis* | Non pathogen | 32 | 5 | 23 | 436 |
| *Mycobacterium leprae* | High Pathogenicity | 1 | 0 | 1 | 689 |
| *Mycobacterium avium* | Low Pathogenicity | 10 | 0 | 9 | 680 |
| *Mycobacterium indicus pranii* | Non pathogen | 8 | 0 | 8 | 591 |
| *Fusobacterium nucleatum* | High pathogenicity | 30 | 0 | 12 | 475 |
| *Fusobacterium necrophorum* | Low pathogenicity | 13 | 0 | 8 | 330 |
| *Leptospira interrogans* | High pathogenicity | 22 | 0 | 20 |1291 |
| *Leptospira biflexa* | Non pathogenic | 19 | 0 | 18 |837 |
We used `grep "toxin-antitoxin" species_name.faa | wc -l` to count the number genes containing toxin-antitoxin modules.
We then used the same comand, but counted the number of hypothetical proteins by running `grep "hypothetical protein" group_project/species_name.faa | wc -l`
Finally, we counted the number of type I and type II by running `grep "type I" group_project/species_name.faa | wc -l` and `grep "type II" group_project/species_name.faa | wc -l`
## OrthoFinder
Goal: Identify missing genes
First, we made directories for each of our phyla containing the FASTA protein sequences for the species in that phyla.
```
mkdir protein_seqs
mkdir protein_seqs/Yersinia
mv protein_seqs/Y.pestis.faa protein_seq/Yersinia
```
Then we ran orthofinder for each of our phyla.
```
cd protein_seqs/Yersinia/
orthofinder -f ./ -X
```
Find the genes unique to each species.
**If using species-specific orthogroups,**
Note: Not all species have species-specific orthogroups based on the N0.tsv results. However, the Statistics_PerSpecies.tsv (path: OrthoFinder/Results_Nov13/Comparative_Genomics_Statistics/Statistics_PerSpecies.tsv) file suggests that there are more species-specific orthogroups than what N0.tsv shows.
Access the OrthoFinder results of COGs containing only a single species:
```
awk -F "\t" '$4=="" && $5=="" && $6!="" {print}' OFS="\t" Yersinia/OrthoFinder/Results_Nov13/Phylogenetic_Hierarchical_Orthogroups/N0.tsv > Yersinia/Y.similis_unique_orthogroups.txt
```
Manipulate the orthogroup data to get a list of the gene ids and remove the trailing commas and hidden character (\r)
```
awk '{for(i=4;i<NF;i++) print $i}' Y.similis_unique_orthogroups.txt | sed 's/,//g' | sed 's/\r//g' > Y.similis_unique_genes.txt
```
To run the process of searching the fasta file and appending the matches to a new file for all of the genes, we wrote a shell script.
```
sh gene_search.sh Y.similis_unique_orthogroups_cleaned.txt Y.similis.faa
```
Shell script (gene_search.sh):
```
#!/bin/bash
# input text file containing a list of gene ids
input=$1
# protein sequence fasta file (.faa) for the species of interest
fasta=$2
# read the input file line by line
for line in `cat "$input"`; do
# search for the gene within the FASTA file
samtools faidx $fasta "$line" >> "$input.faa"
done
```
**If using unassigned genes,**
The OrthoFinder results file containing the unassigned genes:
```
OrthoFinder/Results_Nov13/Orthogroups/Orthogroups_UnassignedGenes.tsv
```
The first column is the Orthogroup id and each of the next columns corresponds to one of the species and contains the gene ids of that species (one per row).
Access the unassigned genes for a particular species, remove the first row, which contains the header, and remove the hidden character (\r):
```
awk -F "\t" '$4!="" {print $4}' OFS="\t" OrthoFinder/Results_Nov13/Orthogroups/Orthogroups_UnassignedGenes.tsv | sed '1d' | sed 's/\r//g' > Y.similis.UnassignedGenes.txt
```
Run a shell script to search for each of those gene ids in the protein sequence FASTA file and append the matching gene sequences to a new file.
```
sh gene_search.sh Y.similis.UnassignedGenes.txt Y.similis.faa
```
*Old shell scripts:*
Shell script version 0:
```
#!/bin/bash
# input file containing a list of gene sequence ids
input=$1
# read the input file line by line
while read -r line
do
# search for the gene within the FASTA file
echo | awk 'BEGIN {RS=">"} /$line/ {print ">" $0}' ORS=" " Y.similis.faa
echo -e "$line\n"
done < "$input"
```
Shell script version 1 (file_parse_temp_2.sh):
```
#!/bin/bash
# input file containing a list of gene sequence ids
input=$1
# read the input file line by line
for line in `cat "$input"`; do
# search for the gene within the FASTA file
samtools faidx Y.similis.faa "$line" >> "$input.output"
done
```
Shell script version 2 (collect_gene_seqs.sh):
Working shell script based on OrthoFinder N0.tsv data (includes manipulating the data)
```
#!/bin/bash
# input text file containing subset of OrthoFinder results file N0.tsv
input=$1
fasta=$2
# parse the orthogroup data to get a list of the gene ids and remove the trailing commas
awk '{for(i=4;i<NF;i++) print $i}' $input | sed 's/,//g' | sed 's/\r//g' > "$input_temp$
# read the input file line by line
for line in `cat "$input_temp_list.txt"`; do
# search for the gene within the FASTA file
samtools faidx $fasta "$line" >> "$input.faa"
done
```
Next steps:
Run rpblast to assign COGs
Look up functional categories for COG number
## rpsBLAST
I ran series of commands
First, `rpsblast -query Y.similis_unique_genes.txt.faa -db /courses/bi278/Course_Materials/blastdb/Cog -evalue 0.01 -max_target_seqs 25 -outfmt "6 qseqid sseqid evalue qcovs stitle" > Y.similis.rpsblast.table` to run rpsblast using the .faa file Olivia obtained using the above shell script and the professor's Cog RPS database.
Next, `awk -F'[\t,]' '!x[$1]++ && $4>=70 {print $1,$5}' OFS="\t" Y.similis.rpsblast.table > Y.similis.cog.table` to keep only the best COG match for each query sequence that have query coverage of at least 70.
Then, `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 Y.similis.cog.table > temp` to get functional categoris for the COG number, which was stored in the class material cognames2003-2014.tab.
After that, `awk -F "\t" '{if ( length($3)>1 ) { $3 = substr($3, 0, 1) } else { $3 = $3 }; print}' OFS="\t" temp > temp2` to keep the first category for each COG number, which is the primary category.
Finally, `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 > Y.similis.cog.categorized` to have a full description of the functional category.
This gave us a list looking like this:
**WP_025384545.1 COG2805 N Cell motility
WP_025381189.1 COG2804 N Cell motility
WP_025381278.1 COG3497 X Mobilome: prophages, transposons
WP_025381364.1 COG1344 N Cell motility
WP_025381362.1 COG1344 N Cell motility
...**
We repeat these commands for each of the species
COG Categories:
A - RNA processing and modification
B - Chromatin structure and dynamics
C - Energy production and conversion
D - Cell cycle control, Cell division, Chromosome partitioning
E - Amino acid transport and metabolism
F - Nucleotide transport and metabolism
G - Carbohydrate transport and metabolism
H - Coenzyme transport and metabolism
I - Lipid transport and metabolism
J - Translation, Ribosomal structure and biogenesis
K - Transcription
L - Replication, Recombination and repair
M - Cell wall/Membrane/Envelope biogenesis
N - Cell motility
O - Posttranslational modification, protein turnover, chaperones
P - Inorganic ion transport and metabolism
Q - Secondary metabolites biosynthesis, Transport and catabolism
R - General function prediction only
S - Function unknown
T - Signal transduction mechanisms
U - Intracellular trafficking, secretion, and vesicular transport
V - Defense mechanisms
W - Extracellular structures
X - Mobilome: prophages, transposons
Y - Nuclear Structure
Z - Cytoskeleton
Shell script to count the number of each functional category present in a given set (named cog_count.sh):
```
#!/bin/bash
# input text file containing the primary functional category for each gene
input=$1
a=0
b=0
c=0
d=0
e=0
f=0
g=0
h=0
i=0
j=0
k=0
l=0
m=0
n=0
o=0
p=0
q=0
r=0
s=0
t=0
u=0
v=0
w=0
x=0
y=0
z=0
# get the list of the functional category alphabetic index
awk '{print $3}' $input > "cog_count_temp.txt"
# read the input file line by line
for line in `cat "cog_count_temp.txt"`; do
# count the number of instances of each functional category
case $line in
"A")
a=$((a+1))
;;
"B")
b=$((b+1))
;;
"C")
c=$((c+1))
;;
"D")
d=$((d+1))
;;
"E")
e=$((e+1))
;;
"F")
f=$((f+1))
;;
"G")
g=$((g+1))
;;
"H")
h=$((h+1))
;;
"I")
i=$((i+1))
;;
"J")
j=$((j+1))
;;
"K")
k=$((k+1))
;;
"L")
l=$((l+1))
;;
"M")
m=$((m+1))
;;
"N")
n=$((n+1))
;;
"O")
o=$((o+1))
;;
"P")
p=$((p+1))
;;
"Q")
q=$((q+1))
;;
"R")
r=$((r+1))
;;
"S")
s=$((s+1))
;;
"T")
t=$((t+1))
;;
"U")
u=$((u+1))
;;
"V")
v=$((v+1))
;;
"W")
w=$((w+1))
;;
"X")
x=$((x+1))
;;
"Y")
y=$((y+1))
;;
"Z")
z=$((z+1))
;;
esac
done
printf "A $a\nB $b\nC $c\nD $d\nE $e\nF $f\nG $g\nH $h\nI $i\nJ $j\nK $k\nL $l\nM $m\nN $n\nO $o\n$
rm cog_count_temp.txt
```
Using the shell script:
```
sh cog_count.sh Y.pestis.Unassigned.cog.categoraized
```
The output is a text file with two columns - the first column contains a single capital letter (A in the first row to Z in the last row) representing the COG categories and the second column is an integer representing the number of genes in the input file with that COG category.
## OrthoFinder with all Taxa Together
Google Drive
Gist
Clipboard
Sign in with Wallet
