# Profile of conservation between bacterial sequences and human protein

My goal is to make a conservation plot between bacterial sequences and a human protein.

So far, I have a FASTA file of the protein, and a FASTA file with the sequences of the proteins from the BLAST results.

In my initial attempts, I have tried to do this using msa software:

I have downloaded clustalx and run a profile alignment using the single protein sequence as Profile 1 and the sequences from the BLAST results as Profile 2, and selected the option "Align Sequences to Profile 1".

I'm having difficulty interpreting and dealing with the results. Since the majority of the sequences aren't aligned, the results are a total mess, with a huge number of dashes added.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

ie:

I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.

Is there a better way to achieve this or to salvage the results from the profile alignment?

Thanks.

• I curious why you want to use a profile alignment. Wouldn't it be easier to just add your sequence to the sequences from the blast results and run a standard multiple sequence alignment? – heathobrien Jul 12 '17 at 8:27
• Also, could you clarify what you're trying to achieve? You say you want the number of blast hits per amino acid, but the figure you show is plotting amino acid conservation. Do you care if the amino acids are the same at a given position or just if the high-scoring segment pairs overlap the position? – heathobrien Jul 12 '17 at 8:30
• We could help much more if you gave us the alignment in text form or even a few of your sequences so we can align them here. The image is informative but we can't really work with it. That said, please do explain what you're trying to do here. I am having trouble understanding what the number of blast hits per amino acid would convey. Why would you need that value? And are you sure you want all of those sequences? If the majority are not aligned, are you sure they're relevant? What is the actual question you are trying to answer using this alignment? – terdon Jul 12 '17 at 9:00
• A profile alignment generates a new alignment from two sets of aligned sequences (the "families"), where one particular case is between one family and one sequence. Therefore to run clustalx, you need to first align the proteins from BLAST results, before using them as your Profile 2. The other questions remain, about number of blast hits or final goal...maybe there is some confusion with the position-specific score matrix from PSI-BLAST, which also generates a score? – Leo Martins Jul 12 '17 at 13:11
• @bluescholar1212 have you looked at Jalview? It is designed as an MSA viewer, plus has many ways for viewing conservation which you can extend yourself. – ithinkiam Jul 13 '17 at 8:39

I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the same genome so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.

You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.

If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py

Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.

For example, given this MSA:

    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH


You'd get this plot:

I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.

What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.

As you say in the comments,

My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

1. Blast sequence of interest.
2. Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
3. Align all the sequences.
4. Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.
• this is an excellent answer. I think that this is the work flow that I will adopt. While researching this I also found this blog post link with a script that produces a pretty smooth conservation plot using the Savitzky-Golay technique. The script seems to follow a similar workflow to what you proposed. Now I will need to work on which blast parameters and metrics of conservation will be appropriate for this project. Thanks for your help. – bluescholar1212 Jul 17 '17 at 17:50
• If you don't mind abstracting the data slightly, the script can calculate a moving average/running mean which smooths the data a bit. – Joe Healey Jul 17 '17 at 19:59
• this may be worth posting a new question to answer, but in the link that I posted, the author mentions you could easily plug in a new function for the remote_blast function depending on usage. I have a list of parameters I would like to use for my blast, but ideally I would like to access the remote database. Do you know how I could change the function in that script to specify more specific parameters that I need? – bluescholar1212 Jul 17 '17 at 21:22
• I'm not sure what you mean exactly? With the python implementations of remote blast you should have most, if not all, the options of a standard blast. What sort of parameters are we talking about? It's easy enough to add additional arguments to the argparse framework if you want to specify blast thresholds etc. It's pretty much just a case of adding in another function or two to do the heavy lifting as far as generating the MSA for you goes, and then that can be passed directly to the other functions of the file. – Joe Healey Jul 18 '17 at 21:23