# working on the server
conda create -n phylo-tutorial-env #python=3.9
# feel free to name the environment whatever you like
# just try to make sure your name is descriptive, so you can remember what it was for
# the `python=3.9` is necessary to ensure a current version of Python is installed, instead of an outdated one
# You'll see a lot of scrolling text, and then need to confirm creation of the environment with "y"
# once setup is done, activate the environment: 
conda activate phylo-tutorial-env
# you can deactivate the environment later with: 
conda deactivate
# remember not to use that above command in your base conda environment!

# in the directory where you have your files for Applied Bioinformatics 
# for ex.: a directory named AppliedBioinfo/
mkdir PhyloWorkflow/
cd PhyloWorkflow/
mkdir Exercise/ Tutorial/
# note that `mkdir` can take multiple arguments
# so we've just created both directories with one line of code

cp -r /cephyr/NOBACKUP/groups/bbt045_2024/Phylogeny/Tutorial/* PATH/TO/YOUR/PhyloWorkflow/Tutorial/
cp /cephyr/NOBACKUP/groups/bbt045_2024/Phylogeny/*.py PATH/TO/YOUR/PhyloWorkflow/Tutorial/

# model usage: 
python clean_MSA_seqs.py input_fasta
# example application: 
python clean_MSA_seqs.py XP_001322682.1__MSAprep.fasta

# Adapt the command above to your file names!
# remember that this cell will assume you are running Python commands - you need to tell it you're running something on the command line

! mafft --localpair --maxiterate 1000 --amino XP_001322682.1__MSAprep_CLEAN.fasta > XP_001322682.1__MSA.fasta
# Adapt the command above to your file names!
# `--localpair --maxiterate 1000` tells MAFFT to use the L-INS-i algorithm
# `--amino` tells MAFFT that the input is a protein FASTA

/opt/alen/bin/alen XP_001322682.1__MSA.fasta
# A peek behind the curtain: 
# I struggled to add alen to the PATH, so it's best to use it by calling the full PATH to the executable within the container

! iqtree -s XP_001322682.1__MSA.fasta --prefix XP_001322682.1__MSA_IQ -m LG+I+R5 -seed 12345 -wbtl -B 1000 -T AUTO -ntmax 2
# Adapt the command above to your file names!
# -s is the option to specify the name of the alignment file that is always required by IQ-TREE to work.
# -m is the option to specify the model name to use during the analysis. 
# The special MFP key word stands for ModelFinder Plus, which tells IQ-TREE to perform ModelFinder 
# and the remaining analysis using the selected model.
# Here, the model to use has been pre-selected: LG+R5
# To make this reproducible, need to use -seed option to provide a random number generator seed.
# -wbtl Like -wbt but bootstrap trees written with branch lengths. DEFAULT: OFF
# -T AUTO: allows IQ-TREE to auto-select the ideal number of threads
# -ntmax: set the maximum number of threads that IQ-TREE can use

%%script false --no-raise-error
! iqtree -s XP_001322682.1__MSA.fasta --prefix XP_001322682.1__MSA_IQ -m LG+I+R5 -seed 12345 -wbtl -B 1000 -T AUTO -ntmax 2
# -B specifies the number of bootstrap replicates where 1000 is the minimum number recommended

# run the python script to visualize your tree in this cell!
! python visualize_PhyloTree_base.py XP_001322682.1__MSA_IQ.treefile

#TrimAl example usage:
trimal -in XP_001322682.1__MSA.fasta -out XP_001322682.1__MSA_trim1.fasta -gappyout
# CIAlign example usage: 
CIAlign --infile XP_001322682.1__MSA.fasta --outfile_stem XP_001322682.1__MSA_trim2 --clean

%%bash
# write your trimming code here

%%bash
# you may want to generate new trees and visualize them here

%%bash
# visualize the bootstrapped trees

# ref: https://stackoverflow.com/questions/73702044/how-to-get-a-consensus-of-multiple-sequence-alignments-using-biopython
# ref: https://biopython.org/docs/1.75/api/Bio.Align.AlignInfo.html

# Part 1: Import necessary modules, determine input & output files

# import relevant modules
import os # allows access to the operating system
# modules necessary for the function to work
from __future__ import annotations
from pathlib import Path
from itertools import chain
# necessary Biopython modules
import Bio
from Bio import AlignIO
from Bio.Align import MultipleSeqAlignment
from Bio.Align.AlignInfo import SummaryInfo


# load input & output files
input_msa = "XP_001322682.1__MSA.fasta"
# extract the file basename
base = os.path.basename(input_msa)
outfile_base = os.path.splitext(base)[0]
# designate the output file nam
output_file = outfile_base + '__ConSeq.fasta'


# create a SeqRecord object with Biopython
SeqRecord = Bio.SeqRecord.SeqRecord


# Part 2: Create the consensus sequence

# define a function to return the consensus sequence
# ref: https://stackoverflow.com/a/73705206/18382033
def get_consensus_seq(filename: Path | str) -> SeqRecord:
    # create an alignment object from the file
    common_alignment = MultipleSeqAlignment(chain(*AlignIO.parse(filename, "fasta")))
    # create a summary object for the MSA
    summary = SummaryInfo(common_alignment)
    # create the consensus sequence
    consensus = summary.dumb_consensus(threshold=0.35)
    # note that dumb_consensus() raises a deprecation error, but you can ignore this
    # return the consensus sequence
    return consensus

# use the function to retrieve the consensus sequence
consensus_seq = get_consensus_seq(input_msa)
# take a look at the sequence
print(consensus_seq)


# Part 3: Write out result in FASTA format

# write out the consensus sequence to a file
with open(output_file, "w") as outfile: 
    # open the output file for writing
    # write the FASTA header line first - you can change this if you'd like
    outfile.write(">XP_001322682.1__MSA__Consensus" + "\n")
    # then write out the consensus sequence to the file
    outfile.write(str(consensus_seq))

# Part 1: Import necessary modules, determine input & output files

# import relevant modules
import os # allows access to the operating system


# set up input & output files
input_seq_file = "XP_001322682.1__MSAprep_CLEAN.fasta"
consensus_file = "XP_001322682.1__MSA__ConSeq.fasta"
# extract the file basename
base = os.path.basename(consensus_file)
outfile_base = os.path.splitext(base)[0]
# get only the original query protein name
outfile_split_base = outfile_base.split("__")[0]
# designate the output file name
output_file = outfile_split_base + '__functQuery.fasta'


# Part 2: Import the query data

# create a list of query proteins
query_list = ["XP_001322682.1", "KAJ4462103.1"]

# create an empty dictionary for the query sequences
fasta_dict = {}

# parse the large fasta file to retrieve the query sequences
with open(input_seq_file, "r") as infile: 
    # open the large sequence file for reading
    for line in infile: 
        # iterate over the file line by line
        if line.startswith(">"): 
            # identify the header lines
            # and save them to a new variable
            header = line
            # check if any of the list elements are in the header line
            # ref: https://stackoverflow.com/questions/6531482/how-to-check-if-a-string-contains-an-element-from-a-list-in-python
            if any(gene_id in header for gene_id in query_list): 
                # use any() & list comprehension to find the headers for the query protein IDs
                # save the sequence which is in the next line of the file
                fasta_seq = next(infile)
                # and add both the header & sequence to the dictionary
                fasta_dict[header] = fasta_seq


# add the consensus sequence to the fasta dictionary
with open(consensus_file, "r") as infile: 
    # open the consensus FASTA sequence file for reading
    for line in infile:
        # iterate over the file line by line
        if line.startswith(">"):
            # identify the sequence header line
            # and save it to a variable
            header = line
        else: 
            # for the sequence line
            # save the sequence to a variable
            fasta_seq = line
    # add the consensus sequence & header to the dictionary
    fasta_dict[header] = fasta_seq


# Part 3: Write out resulting FASTA to file

with open(output_file, "w") as outfile: 
    # open the output file for writing
    for key in fasta_dict.keys():
        # iterate over the dictionary via its keys
        # and write the key:value pair to the outfile
        outfile.write(key + fasta_dict[key])
        # note that since I didn't strip the endline characters, they don't need to be added here

%%bash

# create results directory for EggNOG files
mkdir EggNOG_Results

# dowload results files
wget -P EggNOG_Results FILE/PATH/HERE

%%bash

# create a directory
mkdir IPRScan_Results

# download the files
wget -P EggNOG_Results FILE/PATH/HERE

Phylogeny Tutorial¶

Introduction¶

Setting up the environment¶

(BONUS) Creating the `conda` environment & container¶

(Mandatory) Running the container¶

Setting up your directory system¶

Obtaining & Exploring Data¶

Data Location¶

Preliminary data exploration¶

Finding homologous sequences¶

Cleaning the data¶

Multiple Sequence Alignment¶

Generating the MSA¶

Viewing the MSA¶

Creating the tree¶

Visualizing Trees¶

Editing the MSA¶

AliView¶

Editing using `conda`¶

Comparing trees¶

Bootstapping¶

Applications of phylogenomics¶

Structural Bioinformatics¶

Functional profiling of proteins¶

Summing up¶

Citations¶

Phylogeny Tutorial¶

Introduction¶

Setting up the environment¶

(BONUS) Creating the conda environment & container¶

(Mandatory) Running the container¶

Setting up your directory system¶

Obtaining & Exploring Data¶

Data Location¶

Preliminary data exploration¶

Finding homologous sequences¶

Cleaning the data¶

Multiple Sequence Alignment¶

Generating the MSA¶

Viewing the MSA¶

Creating the tree¶

Visualizing Trees¶

Editing the MSA¶

AliView¶

Editing using conda¶

Comparing trees¶

Bootstapping¶

Applications of phylogenomics¶

Structural Bioinformatics¶

Functional profiling of proteins¶

Summing up¶

Citations¶

(BONUS) Creating the `conda` environment & container¶

Editing using `conda`¶