Zoonomia

The Zoonomia project compares genomes across many mammals across the animal kingdom. Paper is here with citation:

Zoonomia Consortium. A comparative genomics multitool for scientific discovery and conservation. Nature 587, 240–245 (2020). https://doi.org/10.1038/s41586-020-2876-6

The website is also a cool resource.

Approach

The key contributions of the Zoonomia project are to sequence many mammal genomes, and align them to find orthologous genes across species.

The authors used DISCOVAR used to assemble sequences de-novo, and used the Cactus pipeline (paper, code) to assemble large-scale aligned datasets. Tools for reading/modifying format are provided in the hal toolkit. In my opinion, this work is an impressive accomplishment of computational scaling and algorithmic innovation to push the limits of number of genomes processed, and makes smart use of evolutionary trees to make the many-to-many alignment problem more tractable.

The genome10k project seems to be a continuation and future direction for the Zoonomia project.

Browsing Output

(For more information on file formats for NGS, see note on Sequencing File Formats).

The data investigated here is called 241-mammalian-2020v2.hal (browser link). It is 806GB and contains the genomes for 241 mammals.

A naive implementation of alignment would require pairwise alignments, or roughly 30k alignments to do. The Cactus alignment algorithm provides significant optimization as it only requires about 500 alignments, by building an alignment tree - rather than aligning every single genome to all the others, inferred ancestor genomes can be aligned to other ancestor genomes, which allows computing pairwise alignments on-demand in a relatively quick manner. The HAL toolkit provides mechanisms for modifying this tree, as well as simple operations traversing the tree such as liftover to convert gene annotations from one genome to other genomes by traversing the alignment graph (BED, WIG).

The overall tree can be found on the the Zoonomia website. I made a copy with the Tree of Life project here for searching/browsing different branches in this post.

Analysis with HAL and BLAST

The following section attempts to perform some simple analysis on a subset of the Zoonomia dataset.

TP53

The Zoonomia paper hints that large mammals may contain extra copies of the tp53 gene, which could be part of the solution to Peto’s Paradox that cancer rates do not necessarily scale to number of cells in an organism.

The p53 protein is the “guardian of the genome” and it helps to slow reproduction and destroy damaged cells. The gene that codes p53 in humans is TP53 which is given here. P53 is one of the most studied genes, which is why we look at it here in the Zoonomia dataset.

Limiting analysis to Seals using halExtract

Limiting analysis to two similar seals can give us a sense for using HAL tools, and also test a hypothesis that larger sea mammals may have multiple copies of this gene.

Using halExtract we pull out a sub-tree with just two seals, Mirounga angustirostris (Northern Elephant Seal, genome assembly link) and Leptonychotes weddellii (Weddell Seal, arctic, genome assembly link):

$ halExtract --root fullTreeAnc216 241-mammalian-2020v2.hal seals.hal

$ halStats seals.hal confirms we extracted the data correctly.

Finding tp53 in Large Mammal genomes BLAST

Whole-genome analysis can be quite slow overall, so ideally we want to look at some subset of the genome of interest. The “long” way to do this is to do a manual BLAST search for a particular gene (in this case, tp53). We will use BLAST to find the gene within the Homo sapiens Zoonomia genome, then use HAL to compare across species. The main species of interest are the two seal genomes mentioned above, as well as the savanah elephant which is known to have multiple copies of this gene.

To find this gene (or a possibly orthologous one), we start by downloading 25 isoform genes for TP53 from here) and put them into the file tp53.fasta for querying.

To do the actual querying, we create a BLAST databases for each of the genomes of interest. Then we’ll go ahead and search the databases for the 25 isoform genes of tp53 downloaded above.

function make_blast_db() {
    # Pull out a genome into a FASTA file and index it into a BLAST DB.
    GENOME_NAME=$1
    echo "Extracting $GENOME_NAME FASTA File"
    hal2fasta 241-mammalian-2020v2.hal "$GENOME_NAME" > "${GENOME_NAME}.fasta"
    mkdir -p "blast/${GENOME_NAME}"
    makeblastdb -dbtype nucl \
        -in "${GENOME_NAME}.fasta" \
        -out "blast/${GENOME_NAME}/db"
}

function search_blast_db() {
    # Search an existing blast DB for sequences in a given FASTA file.
    GENOME_NAME=$1
    QUERY_FASTA_FILE=$2
    QUERY_NAME=$(basename $QUERY_FASTA_FILE | cut -d. -f1)
    echo "Searching $GENOME_NAME genome for $QUERY_FASTA_FILE"
    blastn \
        -db "blast/${GENOME_NAME}/db" \
        -query "$QUERY_FASTA_FILE" \
        -out "${GENOME_NAME}_${QUERY_NAME}.out"
}

# Index and search for tp53.
$ make_blast_db Mirounga_angustirostris
$ make_blast_db Leptonychotes_weddellii
$ search_blast_db Mirounga_angustirostris tp53.fasta
$ search_blast_db Leptonychotes_weddellii tp53.fasta

# Check matches in human and elephant genomes to confirm positive controls.
$ make_blast_db Homo_sapiens
$ search_blast_db Homo_sapiens tp53.fasta
$ make_blast_db Loxodonta_africana
$ search_blast_db Loxodonta_africana tp53.fasta

The results for searching for isoform g (NM_001276761.3) are:

• For positive control, as expected, the human genome (Homo sapiens) has perfect many perfect alignments within chomosome 17, which matches the expected location of 17p13.1. There are also many smaller sub-matches across other choromosomes. The matching position is on chr17 chromosome from positions 7669699 to 7668421 and positions 1230 to 2509 on the Plus(query)->Minus(`chr17) direction. There are also many hits on the same range of tp53 that is in the seal genome.
• The African Savana Elephant (Loxodonta africana) is said to have multiple copies of tp53, so we use it as a second positive control. I confirmed 16 copies of the gene with regions of ~800bp in length, as well as the same 100 length sequence found in the arctic seal and northern elephant seal.
• The arctic seal (Leptonychotes weddellii) has good alignment in the KB715705.1 800kbp scaffold, with a match length of 115bp identity. The position in the seal genome is KB715705.1 scaffold from indexes 333689 to 333803, which matches positions 807 to 921 in the tp53 isoform g.
• The northern elephant seal (Mirounga angustirostris) has good alignment in the MirAng_flattened_line_76382 15kbp scaffold, with a match length again of 115 with 91% identity. The position is MirAng_flattened_line_76382 scaffold from index 2706 to 2820, which matches positions 807 to 921 in the tp53 isoform g.

Finding tp53 in Seal genomes with halLiftover

The HAL format is optimized for finding similar sequences across genomes.

From the TP53 site we can create a BED file to annotate the Homo_sapiens for the gene in question. (Note that the BED file is expected to be tab-separated). For negative-strand matches returned in BLAST in the form (end, start) we return (start, end) and keep the strand as - per blast_to_bed.py.

# Put the location we found using BLAST into BED format. Use Python to print tabs.
python -c 'print(\
    "\t".join(["chr17", "7668420", "7669699", "tp53", "1000", "-"]))' \
    > tp53_homo_sapiens.bed

# Run Liftover from Homo_sapiens to Leptonychotes_weddellii.
halLiftover --bedType 6 \
241-mammalian-2020v2.hal \
Homo_sapiens tp53_homo_sapiens.bed \
Leptonychotes_weddellii Leptonychotes_weddellii_tp53_liftover.bed 

The result is that halLiftover found a much longer alignment along, 336374 - 335437, but allowed for more gaps. It’s likely that halLiftover was more permissive than with BLAST.

Retrying for positive control with the African elephant, Loxodonta_africana, gives large overlap spanning about 1kbp, but I think due to the way HAL handles the duplicated genomes it only finds one single chromosome rather than the 16 that BLAST found.

halLiftover --bedType 6 \
241-mammalian-2020v2.hal \
Homo_sapiens tp53_homo_sapiens.bed \
Loxodonta_africana Loxodonta_africana_tp53_liftover.bed

Looking at Seal Mutation Stats

Let’s look at the mutations across the species with halSummarizeMutations seals.hal:

To be honest, I’m not really sure what to make of these stats, but we can just keep them here for future reference.

Finding tp53 mutations in all of Zoonomia using halSnps

The analysis now concludes with a search for tp53 across all of Zoonomia mammals.

First, we need to confirm that the alignments we got from BLAST are matching the tp53 gene. To perform a quick sanity check, we export the sequence using hal2fasta and confirm that the resulting FASTA file’s reverse-complement (per bioinformatics.org) is matching the expected region for tp53 within the BLAST output.

hal2fasta --sequence chr17 --start 7668420 --length 1279 --upper \
--outFaPath Homo_sapiens_tp53_in_hal_file.fasta \
241-mammalian-2020v2.hal Homo_sapiens

Comparing the beginning of the match, we see:

BLAST:

Query  1230     CTCACTCCAGCCACCTGAAGTCCAAAAAGGGTCAGTCTACCTCCCGCCATAAAAAACTCA  1289
                ||| || |||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  7669699  CTC-CTACAGCCACCTGAAGTCCAAAAAGGGTCAGTCTACCTCCCGCCATAAAAAACTCA  7669641

RC FAST:        CTC CTACAGCCACCTGAAGTCCAAAAAGGGTCAGTCTACCTCCCGCCATAAAAAACTCA

The end of the match is the same:

BLAST:

Query  2490     ACAATAAAACTTTGCTGCCA  2509
                ||||||||||||||||||||
Sbjct  7668440  ACAATAAAACTTTGCTGCCA  7668421

RC FASTA:       ACAATAAAACTTTGCTGCCA

With more confidence that the extracted region is what we’ve intended it to be, we can run halSnps on our 4 previously studied genomes. Assuming halSnps handles strandedness as BLAST does, it should be fine to search for the reverse complement (seems like it does in the code).

halSnps will report the number of nucleotide matches (totalSnps) and cross-species mismatches (totalCleanOrthologousPairs) for orthologous genes to the target.

We run the command and check the output for positive control:

halSnps --refSequence chr17 --start 7668420 --length 1279 \
241-mammalian-2020v2.hal Homo_sapiens \
Homo_sapiens,Loxodonta_africana,Mirounga_angustirostris,Leptonychotes_weddellii

Predictably, Homo sapiens has a perfect match (so it has zero SNPs with the entire sequence length matching). The other species do contain more matches, which was not something discovered using BLAST. It seems like the orthologous genes may be “forced” in the alignment.

Now that we’ve confirmed reasonable output, we can run halSnps on all the zoonomia mammals. We test that at least 10% of the species (24) must have a mutation at the position for it to be labeled as a SNP, for smoothing. The 10% threshold was selected somewhat arbitrarily but should account for some noise in the sequencing.

function all_genome_names() {
    hal_file=$1
    echo $(halStats $hal_file | grep -v -e "fullTreeAnc\|hal v\|^$\|GenomeName" | cut -d "," -f1 | sort | uniq | paste -s -d, -)
}

halSnps --refSequence chr17 --start 7668420 --length 1279 --minSpeciesForSnp 24 \
241-mammalian-2020v2.hal Homo_sapiens \
"$(all_genome_names 241-mammalian-2020v2.hal)" | sed "s/ /,/g"

Certain primates share nearly identical copies of this gene. Future directions may plot similarity of this gene and others with expected lifespan. It could be interesting to cluster genes detected here, as there seems to be a modal peak with ~800bp in common with the human TP53 gene.

The full output for all 241 mammals is:

Summary

Overall, it seems like HAL and Zoonomia are interesting places to start with comparative genomics. Simply having all the genomes aligned and in one file is quite powerful and combined with AlphaFold (for certain proteins) might yield some results to better understand how mutations may affect function.

Some interesting directions might be to try to take the Fauna Bio approach to compare how hibernating mammals and reptiles can survive in sub-freezing temperatures, and try to understand using halLiftover if there are any genes that these reptiles and mammals share in common. For example, according to this paper on Reptile Freeze Tolerance, reptiles express genes related to iron binding, antioxidant defense, and serine protease inhibitors in their heart and liver after being exposed to freezing temperatures.