<p>That's a LOT of sequences! I'm doubtful you'd be able to get TripletRes or CCMPred to run at all with an MSA depth of 1.5 million. I suggest pruning your dataset to remove redundant sequences (sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA). I'm curious if your Neff (normalized number of effective sequences--as defined in the ResTriplet paper you link) is anywhere close to your number of sequences.</p>
<p>Also, I think you've misunderstood the use of DeepMSA. DeepMSA's main use is to find sequences to include in an MSA. I believe it uses much simpler programs to construct the MSA itself. So you'd use DeepMSA if you wanted to do 1 query sequence --&gt; an MSA of several hundred sequences. Since you already have the (very complete) set of sequences, it wouldn't make sense to use DeepMSA. You could instead just use something like <a href="http://www.clustal.org/omega/" rel="nofollow noreferrer">Clustal Omega</a> to generate the MSA. There's also a pretty easy-to-use <a href="https://www.ebi.ac.uk/Tools/msa/clustalo/" rel="nofollow noreferrer">web server</a> but it's limited to 4k sequences.</p>

CC BY-SA 4.0
2022-01-04T19:27:07.853Z
or have I misunderstood? Are the 1.5M not the same/related protein? If they aren't, I don't see why you'd want to limit your search to the 1.5M only as a larger search would make for a more accurate contact map.
jhschwartz
CC BY-SA 4.0
2022-01-05T02:19:29.427Z
Yes, it's a pretty big dataset - if you're curious, it's all known coronavirus N protein sequences.  I somewhat overestimated initially - it turns out there are only about 77k unique _protein_ sequences in there - taking into account synonymous codons and cases where only a small fragment was sequenced etc.
Alex I
CC BY-SA 4.0
2022-01-05T02:21:36.280Z
"I'm curious if your Neff" - You mean Nf, the same thing that is calculated by calcNf in the DeepMSA software package?  I'm still trying to figure that out, the program to calculate it runs out of memory :)
Alex I
CC BY-SA 4.0
2022-01-05T02:23:47.830Z
"sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA" - is it definitely true they don't add information? in this case the majority are just 1 or 2 aa changes, but there are a lot of them.  my hope was that this would be more informative than using a much smaller dataset of homologs from different species
Alex I
CC BY-SA 4.0
2022-01-05T16:02:37.440Z
You can use sequence weights to investigate how much information you get from a sequence. From the author of DeepMSA, see https://github.com/kad-ecoli/MSAParser/blob/master/calNf_ly.cpp
This should help you downsize your dataset to include the best sequences.
jhschwartz
CC BY-SA 4.0
2022-01-05T16:02:42.757Z
While it's not "definitely true" that a sequence that differs only 1 or 2 AA don't add information, it's likely that the amino acids which change in such a sequence could already change somewhere else in your MSA. Remember that TripletRes and CCMPred look for coevolution in your MSA, so the more variety you have across sequences, the more contacts you can derive from coevolution and the better structure you can model. I strongly recommend using DeepMSA as-is (it should be able to find the sequences in your dataset anyway) to make a good MSA for contact prediction.
jhschwartz
CC BY-SA 4.0
2022-01-05T16:03:04.747Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/132907/discussion-between-jhschwartz-and-alex-i).
jhschwartz
CC BY-SA 4.0
2022-01-04T19:25:03.650Z
2022-01-04T19:25:03.650Z
jhschwartz
<p>You might have some luck with <a href="https://github.com/marbl/parsnp" rel="nofollow noreferrer">Parsnp</a>. From the description:</p>
<blockquote>
<p>Parsnp was designed to align the core genome of hundreds to thousands of bacterial genomes within a few minutes to few hours. Input can be both draft assemblies and finished genomes, and output includes variant (SNP) calls, core genome phylogeny and multi-alignments. Parsnp leverages contextual information provided by multi-alignments surrounding SNP sites for filtration/cleaning, in addition to existing tools for recombination detection/filtration and phylogenic reconstruction.</p>
</blockquote>
<p>The tldr is that Parsnp is meant for aligning large sets of very closely related microbial sequences. For example, others in my lab have used parsnp to align over a million SARS-CoV-2 genomes. The only potential concern is that the Parsnp workflow also finishes with a non-optional phylogeny estimation step that takes quite long on such large datasets. Someone should really make it so that you can opt out of that step...</p>
<p>Spoiler/disclaimer, as the current maintainer of Parsnp, that someone is probably me :grimace:</p>
<p>The output of Parsnp is in <a href="http://darlinglab.org/mauve/user-guide/files.html" rel="nofollow noreferrer">XMFA</a>. You're probably looking for a fasta alignment, so you can use <a href="https://github.com/marbl/harvest-tools" rel="nofollow noreferrer">HarvestTools</a> to convert them:</p>
<pre><code>harvesttools -x sars_align/parsnp.xmfa -M sars_align/parsnp.maf
</code></pre>
<p>Both parsnp and harvest-tools are available on bioconda, and I'd highly recommend installing them that way.</p>

CC BY-SA 4.0
2022-01-07T04:19:39.203Z
2022-01-07T04:19:39.203Z
Throckmorton
<p>I'd like to try a contact prediction method (like TripletRes, CCMpred etc) using a custom set of sequences.  There are about 1.5M sequences, all very similar to the canonical sequence; their distribution reflects the occurence of mutations &quot;in the wild&quot; (they're not clustered or anything like that).  How should I make a MSA in this case?</p>
<p>For example in the original <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6851483/" rel="nofollow noreferrer"><strong>TripletRes paper</strong></a> it looks like the MSA comes from DeepMSA - this searches UniClust30 and UniRef90 using <a href="https://www.openstructure.org/docs/2.2/bindings/hhblits/" rel="nofollow noreferrer"><strong>hhblits</strong></a> and <a href="http://hmmer.org/" rel="nofollow noreferrer"><strong>jackhmmer</strong></a>.  I don't see any way to use DeepMSA to make a MSA from a fully custom set of sequences, and also it's not clear how to combine a custom set with the UniClust/UniRef search results.  I can make an alignment of just my data with jackhmmer, but I don't know if that's what the other tools expect.  For one thing, the alignment is pretty huge compared to the typical inputs (~million vs typical DeepMSA outputs of a few hundred); for another, it is a lot less variable (most sequences differ in just 1-2 aa, vs completely different species).</p>

<protein-structure><multiple-sequence-alignment>
How to make MSA for contact prediction?
# How to make MSA for contact prediction?

(Asked by: username_0 on 2022-01-01)

I'd like to try a contact prediction method (like TripletRes, CCMpred etc) using a custom set of sequences. There are about 1.5M sequences, all very similar to the canonical sequence; their distribution reflects the occurence of mutations "in the wild" (they're not clustered or anything like that). How should I make a MSA in this case?

For example in the original [**TripletRes paper**](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6851483/) it looks like the MSA comes from DeepMSA - this searches UniClust30 and UniRef90 using [**hhblits**](https://www.openstructure.org/docs/2.2/bindings/hhblits/) and [**jackhmmer**](http://hmmer.org/). I don't see any way to use DeepMSA to make a MSA from a fully custom set of sequences, and also it's not clear how to combine a custom set with the UniClust/UniRef search results. I can make an alignment of just my data with jackhmmer, but I don't know if that's what the other tools expect. For one thing, the alignment is pretty huge compared to the typical inputs (~million vs typical DeepMSA outputs of a few hundred); for another, it is a lot less variable (most sequences differ in just 1-2 aa, vs completely different species).

**Tags:** protein-structure, multiple-sequence-alignment

---

## 2 answers

### Answer by username_1 on 2022-01-04. Score: 1

That's a LOT of sequences! I'm doubtful you'd be able to get TripletRes or CCMPred to run at all with an MSA depth of 1.5 million. I suggest pruning your dataset to remove redundant sequences (sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA). I'm curious if your Neff (normalized number of effective sequences--as defined in the ResTriplet paper you link) is anywhere close to your number of sequences.

Also, I think you've misunderstood the use of DeepMSA. DeepMSA's main use is to find sequences to include in an MSA. I believe it uses much simpler programs to construct the MSA itself. So you'd use DeepMSA if you wanted to do 1 query sequence --> an MSA of several hundred sequences. Since you already have the (very complete) set of sequences, it wouldn't make sense to use DeepMSA. You could instead just use something like [Clustal Omega](http://www.clustal.org/omega/) to generate the MSA. There's also a pretty easy-to-use [web server](https://www.ebi.ac.uk/Tools/msa/clustalo/) but it's limited to 4k sequences.

> **Comments:**
> **username_1**: or have I misunderstood? Are the 1.5M not the same/related protein? If they aren't, I don't see why you'd want to limit your search to the 1.5M only as a larger search would make for a more accurate contact map.  
> **username_0**: Yes, it's a pretty big dataset - if you're curious, it's all known coronavirus N protein sequences.  I somewhat overestimated initially - it turns out there are only about 77k unique _protein_ sequences in there - taking into account synonymous codons and cases where only a small fragment was sequenced etc.  
> **username_0**: "I'm curious if your Neff" - You mean Nf, the same thing that is calculated by calcNf in the DeepMSA software package?  I'm still trying to figure that out, the program to calculate it runs out of memory :)  
> **username_0**: "sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA" - is it definitely true they don't add information? in this case the majority are just 1 or 2 aa changes, but there are a lot of them.  my hope was that this would be more informative than using a much smaller dataset of homologs from different species  
> **username_1**: You can use sequence weights to investigate how much information you get from a sequence. From the author of DeepMSA, see https://github.com/kad-ecoli/MSAParser/blob/master/calNf_ly.cpp This should help you downsize your dataset to include the best sequences.  
> **username_1**: While it's not "definitely true" that a sequence that differs only 1 or 2 AA don't add information, it's likely that the amino acids which change in such a sequence could already change somewhere else in your MSA. Remember that TripletRes and CCMPred look for coevolution in your MSA, so the more variety you have across sequences, the more contacts you can derive from coevolution and the better structure you can model. I strongly recommend using DeepMSA as-is (it should be able to find the sequences in your dataset anyway) to make a good MSA for contact prediction.  
> **username_1**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/132907/discussion-between-jhschwartz-and-alex-i).  

---

### Answer by username_2 on 2022-01-07. Score: 1

You might have some luck with Parsnp. From the description:

Parsnp was designed to align the core genome of hundreds to thousands of bacterial genomes within a few minutes to few hours. Input can be both draft assemblies and finished genomes, and output includes variant (SNP) calls, core genome phylogeny and multi-alignments. Parsnp leverages contextual information provided by multi-alignments surrounding SNP sites for filtration/cleaning, in addition to existing tools for recombination detection/filtration and phylogenic reconstruction.

The tldr is that Parsnp is meant for aligning large sets of very closely related microbial sequences. For example, others in my lab have used parsnp to align over a million SARS-CoV-2 genomes. The only potential concern is that the Parsnp workflow also finishes with a non-optional phylogeny estimation step that takes quite long on such large datasets. Someone should really make it so that you can opt out of that step...
Spoiler/disclaimer, as the current maintainer of Parsnp, that someone is probably me :grimace:
The output of Parsnp is in XMFA. You're probably looking for a fasta alignment, so you can use HarvestTools to convert them:
harvesttools -x sars_align/parsnp.xmfa -M sars_align/parsnp.maf

Both parsnp and harvest-tools are available on bioconda, and I'd highly recommend installing them that way.


---

