I read lecture notes about mutations, and wondered what kind of algorithms are there to detect mutations?

How do we know if the gene is mutated or whether it is a sequencing error?

I saw this thread which is related, but I am not sure how to use the CIGAR information, e.g., is it 100% accurate? What is the underline mechanism to detect mutation? Is there a way to predict what the mutation will cause?


I'll follow up to the great answer from Kamil S Jaron:

Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.

If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.

  • $\begingroup$ Nice summary, to add, PhastCons is the standard for conservation scores, although not the best as the Sipel lab has many newer tools alternatives and GERP++ is also better IMO. Also if you're looking at just a few individual loci, nothing beats a stare at a genome browser like UCSC $\endgroup$ – Chris_Rands Jun 11 '17 at 16:47

Part 1 : how to detect mutations

The keywords you are searching for are "variant calling". Basically you have to map sequencing reads to a reference genome (or gene) and then estimate for each position of the genome if the observed difference of mapped reads and the reference is more likely a sequencing error or a mutation (in genomic glossary - variant).

Popular tools for variant calling are GATK, FreeBayes or bcftools (previously part as Samtools package).

The question you linked asks for quick alternatives to simply looking for variants. Indeed, you can just visualise the mapped reads to the reference sequence and see if the variant is there or not. CIGAR is just a notation of read alignment used in sam files (files with mapped reads), you can find good explanation of CIGAR strings here.

The follow-up about "How to estimate an effect of mutation" is in @DevonRyan 's awesome answer.


Part 3

I thought I'd add a part 3 . Having found via VEP a truncation (missense or frameshift) or a missense, how do you find what are the consequences at the protein level?

The effect a mutation has on a protein can be diverse and how much imbalance in function can the cell's equilibrium absorb depends on the component. As a result you get recessive, dominant and embryonically lethal mutations. In this figure are shown some of the effects, divided into core and surface mutations: effects

Core variants

Core mutations are the most black-and-white to analyse, but require a protein structure available. A severe truncation and a highly deleterious missense will result in a protein that will likely be degraded. The difference in folding Gibbs free energy (∆∆G, where $\Delta \Delta G = \Delta G_{\rm{mutant}} – \Delta G_{\rm{WT}}$)) can be calculated for a missense, but how large this must be to result in a certain percentage being degraded depends on factors that cannot be calculated analytically. Alternatively, the Gibbs free energy of ligand binding is affected resulting in loss of catalysis.

To calculate the ∆∆G effect (stability), there are many tools:

  • Online SDM uses ML to estimate the ∆∆G quickly (do note that positive kcal/mol normally is a loss of stability, but here the sign is inverted so positive is a gain in stability
  • Rosetta toolkit. A highly flexible panel of movers (samplers) can calculate ∆∆G
  • FoldX
  • MD runs
  • Missense3D assesses a variant based on several structural criteria

In the case of ligand binding modelling as above but using the ligand bound structure and calculating the binding energy as the difference between bound and unbound (far far away). For altered catalysis, modelling as above but using the transition state and assessing the difference in the transition hill can be done. However, an active site variant is very likely to affect either ligand binding, transition state binding (catalysis) or product binding that it's surely deleterious.

This is all (mostly) straightforward. If the domain is a regulatory domain however, there is loss of inhibition, which results in a gain of function (dominant variant).

Surface variants

Then stuff gets messier when you get to surface variants and these can cause dominant effects. However, if you are close to a phosphorylation, ubiquitination or other PTM site, as identified in phosphosite plus, you likely have a case of deregulation. For linear motifs, ELM is useful, but requires a lot of literature scanning to confirm any find. Alternative a loss of an interface may break the protein function, however, what binds where can only be ascertained by available structures, cross-linking mass spec or in some cases well defined Western experiments, so this data is rare. For this type of variant there are two pipelines:

  • Miscast collates information on the protein
  • VENUS (beta) calculates ∆∆G and finds any relevant close-by PTM or gnomAD variants [disclaimer: this is mine].


As the variants can have a variety of effects, none of these score (unlike the genetics parts), but instead try to inform as to what may be going on.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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