2
$\begingroup$

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 it looks like the MSA comes from DeepMSA - this searches UniClust30 and UniRef90 using hhblits and jackhmmer. 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).

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

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 to generate the MSA. There's also a pretty easy-to-use web server but it's limited to 4k sequences.

Improve this answer
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – jhschwartz
    Jan 4, 2022 at 19:27
  • $\begingroup$ 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. $\endgroup$
    – Alex I
    Jan 5, 2022 at 2:19
  • $\begingroup$ "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 :) $\endgroup$
    – Alex I
    Jan 5, 2022 at 2:21
  • $\begingroup$ "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 $\endgroup$
    – Alex I
    Jan 5, 2022 at 2:23
  • $\begingroup$ You can use sequence weights to investigate how much information you get from a sequence. From the author of DeepMSA, see github.com/kad-ecoli/MSAParser/blob/master/calNf_ly.cpp This should help you downsize your dataset to include the best sequences. $\endgroup$
    – jhschwartz
    Jan 5, 2022 at 16:02
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.