Diffing Coronaviruses | Hacker News Reader

Diffing Coronaviruses

jameshfisher.com

119 pointsjamesfisher17 comments

diffing is not the right approach, it can't detect smaller blocks moving about.

since the genome is made out of only A,T,G,C for any change beyond the trivial diffing may generate the incorrect picture, it lines up the "wrong" bases.

Here is a much simpler solution that does a much better job and you don't have to muck around so much. First install some tools:

  conda install -c bioconda entrez-direct emboss

Now do a:

  # Get the data.
  efetch -db nuccore -id MG772933 -format fasta > MG772933.fa
  efetch -db nuccore -id MN988713 -format fasta > MN988713.fa
  
  # Align globally.
  stretcher -filter  MG772933.fa MN988713.fa  > aligned.txt

A much simpler explanation of what is going on (I picked a variable region here):

  ...
  MG7729 CACTCCTCACTATTCATAGAGGA-----GAC-CCT----ATGCCTAATAA
         :  : ::  :: : :::::: :      ::: :::    :: : :  : :
  MN9887 CTTTACTTGCTTTACATAGAAGTTATTTGACTCCTGGTGATTCTTCTTCA
  ...

it also reports:

  # Length:           29922
  # Identity:   26272/29922 (87.8%)
  # Gaps:         160/29922 ( 0.5%)

An example how domain knowledge makes data wrangling more effective and/or valid.

when all you have is a HMMER...

I did something similar recently using my custom terminal diff tool. Here is a more colourful terminal-based view of the differences between MN908947.3 and sars NC_004718.3:

https://twitter.com/RobertElderSoft/status/12226906294923796...

To produce that image I used this tool:

https://github.com/RobertElderSoftware/roberteldersoftwaredi... With these two genomes:

https://www.ncbi.nlm.nih.gov/nuccore/MN908947

https://www.ncbi.nlm.nih.gov/nuccore/NC_004718.3

with the leading numbers, spaces and newlines removed from each genome file.

My tool (as well as unix diff with --minimal) uses the myers diff algorithm which just compute a mathematically minimal diff. More advanced algorithms exist for computing phylogenetic trees that take into account the biological likelihood of certain sequence change (deletions vs. additions vs. translations etc.).

diff is a special case of a more or less obligate preprocessing step for sequence alignment algorithms.

You need to go further than diff to make meaningful phylogenetic comparisons, though!

https://en.wikipedia.org/wiki/Longest_common_subsequence_pro... has much more.

FYI—sequence alignment algorithms are really just a specialized subset of diffing algorithms that take into account the properties of biological sequences. (For example, certain edits are more or less likely to happen.)

https://en.wikipedia.org/wiki/Sequence_alignment

I would say that perhaps diff is a special case of alignments, but not in reverse.

Alignments differ from a "diff" because they are assigning a score to each operation (match, mismatch, gap). The "optimal" alignment only cares about maximizing the score, does not "care" about the objects that are aligned.

If one chooses the scores incorrectly they could end up with optimal (relative to the score) alignments that visually seem completely "wrong".

I can't see how one could tune a diff, to produce obviously "incorrect" output.

Maybe most correct is to say diff and seq alignments are specializations of some superset, given the specializations for both situations.

Diff algorithms absolutely do assign a "score" to edits, just in a very simple way. For example, by playing with the Myers diff algorithm, you can get an edit distance of one where you simply rewrite the entire file, or a large edit distance but with the smallest edits possible.

Definitely not a great approach.

It's worth pointing toward the actual algorithms commonly used for local and global alignments (global is the one to use in the coronavirus genome case).

Local alignments: Smith-Waterman [1]. Global alignments: Needleman–Wunsch [2]

[1] https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorit... [2] https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algor...

I've been using a kmer approach, just a few lines in Python, and it's interesting comparing coronavirus to many different genomes.

The neat thing about kmers is that it is alignment free, so can find similarities between genomes even if genes have been moved around.

Interestingly, it finds HIV is quite close to cov compared to a random sampling of refseq viruses. Maybe there is credence to that Indian paper after all.

I read once that we share 99% of our DNA with Chimpanzees but we also share 95% of DNA with fruit flies. So if the viruses are 89% similar does that mean the difference between SARs and nCOV is bigger than the difference between a human and fruit fly? Who knows what this analysis might be telling us.

That's not unreasonable to say, but virus sequences drift much faster and have a much shorter lifecycle. This means a 1% difference can be huge for us and insignificant for them.

the important thing to remember here is that there is identity and similarity. Depending on how we define similarity the percentages vary wildly.

Can we find millions of similar regions in the fruitfly as in the human? Yes. Are those regions on the same chromosome and on the same relative orientations - not even close.

If you start reading the genome from left to right, the human and fruitfly are not even remotely similar. If were to chop the DNA into 100bp pieces now you can find a lot of matches across the organisms, still many sequences will be very different. One basepair difference could lead to a completely different protein. Then the fruit fly genome is 300 million bases long, the human is 3 billion bases.

You can see how there is a lot more to similarity than initially meets the eye.

We also have more genes overall, so in absolute terms the difference between us and chimpanzees or fruit flies is pretty large compared to the differences between those two viruses. Also we basically have all the ingredients to procreate ourselves (we only need to team up with someone from the other sex), while viruses partly rely on the protein synthesis machinery of the host. So viruses can do away with conserved proteins they will find in most of their hosts, we can't.

This is a good start. However when trying to build phylogenies and understand evolutionary relationshios there are many next steps beyond sequence identity.

Interesting, in light of a recent discussion about viruses having very high gene variability: https://news.ycombinator.com/item?id=22276105

It seems NCBI's index is not updated. You should try MN996532, which is also from bats and the sequence was uploaded Jan 29. It shows 96% similarity.

Reply on news.ycombinator.com