<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Multisequence Alignment (MSA)
Proteins are made up of sequences of amino acids chained together. Their amino acid sequence determines their structure and function. Finding proteins with similar sequences, or homologous proteins, is very useful in identifying the structures and functions of newly discovered proteins as well as identifying their ancestry. Below is an example of what a protein amino acid multisequence alignment may look like, taken from [2].
Colab
This tutorial and the rest in this sequence can be done in Google colab. If you'd like to open this notebook in colab, you can use the following link.
HH-suite
This tutorial will show you the basics of how to use hh-suite. hh-suite is an open source package for searching protein sequence alignments for homologous proteins. It is the current state of the art for building highly accurate multisequence alignments (MSA) from a single sequence or from MSAs.
References
Step1:  Be sure to input 'y' into the cell below to allow conda to install hh-suite to this colab.
Step2:  Using hhsearch
hhblits and hhsearch are the main functions in hhsuite which identify homologous proteins. They do this by calculating a profile hidden Markov model (HMM) from a given alignment and searching over a reference HMM proteome database using the Viterbi algorithm. Then the most similar HMMs are realigned and output to the user. To learn more, check out the original paper in the references above.
Run a function from hhsuite with no parameters to read its documentation.
Step3:  Let's do an example. Say we have a protein which we want to compare to a MSA in order to identify any homologous regions. For this we can use hhsearch. 
Let's start by downloading a database to work with. hh-suite provides a set of HMM databases that will work with the software, which you can find here
Step5:  Now let's take some protein sequence and search through the dbCAN database to see if we can find any potential homologous regions. First we will specify the sequence and save it as a FASTA file or a3m file in order to be readable by hhsearch. I pulled this sequence from the example query.a3m in the hhsuite data directory.
Step6:  Then we can call hhsearch, specifying the query sequence with the -i flag, the database to search through with -d, and the output with -o.
Step8:  The 'Prob' column describes the estimated probability of the query sequence being at least partially homologous to the template. Probabilities of 95% or more are nearly certain, and probabilities of 30% or more call for closer consideration. The E value tells you how many random matches with a better score would be expected if the searched database was unrelated to the query sequence. These results show that none of the sequences align well with our randomly chosen protein, which is to be expected because our query sequence was chosen at random.
Now let's check the results if we use a sequence that we know will align with something in the dbCAN database. I pulled this protein from the dockerin.faa file in dbCAN.
Step9:  As you can see, there are 2 sequences which are a match for our query sequence. 
Using hhblits
hhblits works in much the same way as hhsearch, but it is much faster and slightly less sensitive. This would be more suited to searching very large databases, or producing a MSA with multiple sequences instead of just one. Let's make use of that by using our query sequence to create an MSA. We could then use that MSA, with its family of proteins, to search a larger database for potential matches. This will be much more effective than searching a large database with a single sequence.
We will use the same dbCAN database. I will pull a glycoside hydrolase protein from UnipProt, so it will likely be related to some proteins in dbCAN, which has carbohydrate-active enzymes.
The option -oa3m will tell hhblits to output an MSA as an a3m file. The -n option specifies the number of iterations. This is recommended to keep between 1 and 4, we will try 2.
<ASSISTANT_TASK:>
Python Code:
!wget -c https://repo.anaconda.com/miniconda/Miniconda3-4.5.4-Linux-x86_64.sh
!chmod +x Miniconda3-4.5.4-Linux-x86_64.sh
!bash ./Miniconda3-4.5.4-Linux-x86_64.sh -b -f -p /usr/local
Explanation: Multisequence Alignment (MSA)
Proteins are made up of sequences of amino acids chained together. Their amino acid sequence determines their structure and function. Finding proteins with similar sequences, or homologous proteins, is very useful in identifying the structures and functions of newly discovered proteins as well as identifying their ancestry. Below is an example of what a protein amino acid multisequence alignment may look like, taken from [2].
Colab
This tutorial and the rest in this sequence can be done in Google colab. If you'd like to open this notebook in colab, you can use the following link.
HH-suite
This tutorial will show you the basics of how to use hh-suite. hh-suite is an open source package for searching protein sequence alignments for homologous proteins. It is the current state of the art for building highly accurate multisequence alignments (MSA) from a single sequence or from MSAs.
References:
[1] Steinegger M, Meier M, Mirdita M, Vöhringer H, Haunsberger S J, and Söding J (2019) HH-suite3 for fast remote homology detection and deep protein annotation, BMC Bioinformatics, 473. doi: 10.1186/s12859-019-3019-7
[2] Kunzmann, P., Mayer, B.E. & Hamacher, K. Substitution matrix based color schemes for sequence alignment visualization. BMC Bioinformatics 21, 209 (2020). https://doi.org/10.1186/s12859-020-3526-6
Setup
Let's quickly install anaconda and use it to properly install hhsuite.
End of explanation
!conda install -c conda-forge -c bioconda hhsuite  
Explanation: Be sure to input 'y' into the cell below to allow conda to install hh-suite to this colab.
End of explanation
!hhsearch
Explanation: Using hhsearch
hhblits and hhsearch are the main functions in hhsuite which identify homologous proteins. They do this by calculating a profile hidden Markov model (HMM) from a given alignment and searching over a reference HMM proteome database using the Viterbi algorithm. Then the most similar HMMs are realigned and output to the user. To learn more, check out the original paper in the references above.
Run a function from hhsuite with no parameters to read its documentation.
End of explanation
%%bash
cd /content/
mkdir hh
cd hh 
mkdir databases; cd databases
wget http://wwwuser.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs/dbCAN-fam-V9.tar.gz
tar xzvf dbCAN-fam-V9.tar.gz
Explanation: Let's do an example. Say we have a protein which we want to compare to a MSA in order to identify any homologous regions. For this we can use hhsearch. 
Let's start by downloading a database to work with. hh-suite provides a set of HMM databases that will work with the software, which you can find here: http://wwwuser.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs
dbCAN is a good one for this tutorial because it is a small download.
End of explanation
with open('protein.fasta', 'w') as f:
    f.write(
>Uncharacterized bovine protein (Fragment)
--PAGGQCtgiWHLLTRPLRP--QGRLPGLRVKYVFLVWLGVFAGSWMAYTHYSSYAELCRGHICQVVICDQFRKGIISGSICQDLCHLHQVEWRTCLSSVPGQQVYSGLWQGKEVTIKCGIEESLNSKAGSDGAPRRELVLFDKPSRGTSIKEFREMTLSFLKANLGDLPSLPALVGRVLLMADFNKDNRVSLAEAKSVWALLQRNEFLLLLSLQEKEHASRLLGYCGDLYVTEGVPLSSWPGATLPPLLRPLLPPALHGALQQWLGPAWPWRAKIAMGLLEFVEDLFHGAYGNFYMCETTLANVGYTAKYDFRMADLQQVAPEAAVRRFLRGRRCEHSADCTYGRDCRAPCDTLMRQCKGDLVQPNLAKVCELLRDYLLPGAPAALRPELGKQLRTCTTLSGLASQVEAHHSLVLSHLKSLLWKEISDSRYT
)
Explanation: Now let's take some protein sequence and search through the dbCAN database to see if we can find any potential homologous regions. First we will specify the sequence and save it as a FASTA file or a3m file in order to be readable by hhsearch. I pulled this sequence from the example query.a3m in the hhsuite data directory.
End of explanation
!hhsearch -i /content/protein.fasta -d /content/hh/databases/dbCAN-fam-V9 -o /content/protein.hhr
Explanation: Then we can call hhsearch, specifying the query sequence with the -i flag, the database to search through with -d, and the output with -o.
End of explanation
with open('protein2.fasta', 'w') as f:
    f.write(>dockerin,22,NCBI-Bacteria,gi|125972715|ref|YP_001036625.1|,162-245,0.033
SCADLNGDGKITSSDYNLLKRYILHLIDKFPIGNDETDEGINDGFNDETDEDINDSFIEANSKFAFDIFKQISKDEQGKNVFIS
)
    
!hhsearch -i /content/protein2.fasta -d /content/hh/databases/dbCAN-fam-V9 -o /content/protein2.hhr
Explanation: The 'Prob' column describes the estimated probability of the query sequence being at least partially homologous to the template. Probabilities of 95% or more are nearly certain, and probabilities of 30% or more call for closer consideration. The E value tells you how many random matches with a better score would be expected if the searched database was unrelated to the query sequence. These results show that none of the sequences align well with our randomly chosen protein, which is to be expected because our query sequence was chosen at random.
Now let's check the results if we use a sequence that we know will align with something in the dbCAN database. I pulled this protein from the dockerin.faa file in dbCAN.
End of explanation
!wget -O protein3.fasta https://www.uniprot.org/uniprot/G8M3C3.fasta
!hhblits -i /content/protein3.fasta -d /content/hh/databases/dbCAN-fam-V9 -oa3m query.a3m -n 2 
Explanation: As you can see, there are 2 sequences which are a match for our query sequence. 
Using hhblits
hhblits works in much the same way as hhsearch, but it is much faster and slightly less sensitive. This would be more suited to searching very large databases, or producing a MSA with multiple sequences instead of just one. Let's make use of that by using our query sequence to create an MSA. We could then use that MSA, with its family of proteins, to search a larger database for potential matches. This will be much more effective than searching a large database with a single sequence.
We will use the same dbCAN database. I will pull a glycoside hydrolase protein from UnipProt, so it will likely be related to some proteins in dbCAN, which has carbohydrate-active enzymes.
The option -oa3m will tell hhblits to output an MSA as an a3m file. The -n option specifies the number of iterations. This is recommended to keep between 1 and 4, we will try 2.
End of explanation
<END_TASK>