# Lab #4 Multiple Sequence Alignment
## 1.Generate an MSA
When we align DNA bases or amino acid residues of two or more genes together, we are saying that not only are these genes homologous to each other, but the position along an alignment are also homologous to each other. To find homologous positions between two genes, an algorithm will apply a scoring system that rewards higher points to identical or similar bases and residues, and lower or negative points for differences including indels. Alignments are hypotheses and we an choose differnt algorithmic approaches or scroing schemes to represnet best how i think sequences evolve.We will compare two popular alignment software: muscle and mafft. both similar but maft has moe option if i have more complex datasets tp work with

1.Find the gene ids of all tssB sequenecs in your new genome that you annotated with prokka last week.I want to work with bbqs395
do this by:

grep tssB PATH/PROKKA_*.tsv

in my case this looked like ( had to get from course mateials hence the path):

grep tssB /courses/bi278/Course_Materials/lab_04/bbqs395/PROKKA_09162022*.tsv 

- the first column contains my gene id, Because of the way the tsv file is formatted, you probably have duplicate lines for gene and CDS. Just get the unique gene ids here.

my gene ids are:
BO395_00972 gene 579 tssB-5_1
BO395_02489 gene 507 tssB-5_2

2.Get fasta sequences for these genes, for protein coding genes, it always better to get amino acid sequences(exons)
I had to copy the ffa file from courses, i did this by:
cp /courses/bi278/Course_Materials/lab_04/bbqs395/PROKKA_09162022.faa .

**To get fatsa sequences, we need to index the appropiate fasta file that contains amino acid sequences for genes predicted in my genome we use a software called samtools for this:**
samtools faidx PATH/PROKKA_*faa
in my case, this looked like:
samtools faidx ./PROKKA_09162022.faa

3.**To call up fatsa format sequences based on their newly indexed gene ids( we are using gene ids found in step 1)**:
samtools faidx PATH/PROKKA_*faa GENEID1 GENEID2 ...
in my case this looked like:
samtools faidx ./PROKKA_09162022.faa BO395_00972 BO395_02489

resluts will look like a fatsa file eg.
>GENEID1
SEQUENCESEQUENCESEQUENCESEQUENCE
>GENEID2
SEQUENCESEQUENCESEQUENCESEQUENCE

4.we want to use these sequences to generate an alignment so we need to save them into a file.We can use the > operator to do this to write over whatever filename you put after it:
samtools faidx PATH/PROKKA_*faa GENEID1 GENEID2 ... > FILENAME
inmy case this looked like:
samtools faidx ./PROKKA_09162022.faa BO395_00972 BO395_02489 > PKG1_2

5.concatenate my tssB sequences with the ones downloaded from the T6SSS database:
cat /courses/bi278/Course_Materials/lab_04/t6ss_db.faa FILENAME > tssB_input.faa
in my case this looked like:
cat /courses/bi278/Course_Materials/lab_04/t6ss_db.faa PKG1_2 > tssB_input.faa

6.Align all sequences with muscle:
muscle -align tssB_input.faa -output tssB_muscle.afa
in my case this looked like:
muscle -align tssB_input.faa -output tssB_muscle.afa

7.Muscle has weird behavior of changing the order of sequences in its alignment, and we wnat to correct this.**To can check behavior itself by comparing the order of sequences by their headers:**
grep ">" tssB_input.faa | head
grep ">" tssB_muscle.afa | head

8.**We can use python script to reorder the sequences**:
python /courses/bi278/Course_Materials/lab_04/stable_py3.py tssB_input.faa tssB_muscle.afa | grep ">"
| head
^THIS DID NOT WORK SO I USED PYTHON3 INSTEAD OF PYTHON:
python3 /courses/bi278/Course_Materials/lab_04/stable_py3.py tssB_input.faa tssB_muscle.afa | grep ">" | head

9.**We have fixed the order now lets save it to a new file**:
python3 /courses/bi278/Course_Materials/lab_04/stable_py3.py tssB_input.faa tssB_muscle.afa >
tssB_muscle.faa
(note difference in afa to faa file, could change output filename if too similar) in my case:
python3 /courses/bi278/Course_Materials/lab_04/stable_py3.py tssB_input.faa tssB_muscle.afa > tssB_muscle1.faa

10.**Align with mafft(maaft doesnt not have reorder of sequencing problem like muscle)**:
mafft --maxiterate 3 tssB_input.faa > tssB_mafft.faa

# 2. Compare MSAs
we can now compare the two alignments using Multiple Sequence Alignmnet Viewer from NCBI

1.Click on Text in DataSource and copy and paste the text of an alignment file in the window. How I got the data to copy from a sequence file (copy them in seperate, open another window for other file, press upload, then click close):
cat tssB_mafft.faa
cat tssB_muscle1.faa

Q1. Both mucle and mafft, have large gaps at before sequence id Bmall 2 to accomodate for all sequencesto align.Both starts at 1.I would assume that the regions that are similar across my alignments share some type of ortholgy

Q2.With mfft, gaps do nyo often end the sequences as often as the alignment in muscle does.Mucsle has sugnificant gaps both at the start and at the end of most sequences, whereas mafft has more gaps at the end of sequences.

Q3.How dod you decide which alingmnet is better?
Whichever you decide is more convncing based on your philosophy

# 3.Generate a gene tree from you MSA
Alignments are usually not the final step of a sequence analysis. The next step for many alignments is to generate a phylogeny. Let’s make a quick phylogeny using approximate Maximum Likelihood with FastTree.

1.**Generate phylogenies from my alignments**(in my case, i had the muscle alignmnets under a different name which was: FastTree -lg < tssB_muscle1.faa > tssB_muscle1_ft.tre):
FastTree -lg < tssB_muscle.faa > tssB_muscle_ft.tre
FastTree -lg < tssB_mafft.faa > tssB_mafft_ft.tre

*FastTree uses a syntax that looks like this:
command < inputfile > outputfile
where it is using < and > to indicate to what to read in and what to write out.*

2.Now lets look at these trees by importing them one at a time.Go to Interactive Tree of Life. Go to upload, then copy and paste, then click upload. Get your tree file content using this command to get the content
cat tssB_mafft_ft.tre
cat tssB_muscle1_ft.tre

3.With these trees, I find that at least initially, looking at them in circular mode and ignoring branch lengths makes it easier to understand what I’m looking at. You can change these options in the [Control panel].Try different options here, and make sure you can find the sequences that you added to the larger number of sequences in the database file.

Just learned how to make and inspect alignments from homologous sequences, and generate and view trees