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.
Profile of conservation between bacterial sequences and human protein
sequence alignment
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.pyJust 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    149PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTGPLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTGPAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTGPLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTGPAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTGPAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTGPLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTGPLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTGPLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIGPAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVGPAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVGPAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIGPAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIGPAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVGPAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVGPLT_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 ADRVTVEFHYou'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:Blast sequence of interest.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.)Align all the sequences.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.