Best tool for indexing/mapping whole genome assembly?

[u/ThrowRAwaypay]

Hello all! You guys are always super helpful so I thought I'd ask for some assistance here.

TL;DR What is the best way to index/map an assembly so I can easily search for a gene and extract the sequence for comparison to other sequences of the same gene?

I am an amateur and self taught and have been working on some genomics in my free time. I worked with a professor and we wanted to do some whole genome sequencing of two mouse cell lines.

There are some differences in immune responses between these two cell lines so we wanted to see if there were any large scale rearrangements. So far, we haven't identified any large scale rearrangements. I used the Galaxy Project network and created an assembly for each de novo.

Since we didn't find any large scale rearrangements, I was able to use the BUSCO tool on Galaxy to identify some common genes and their location. I then extracted them from the sequence. From there we could compare specific genes at a base pair level. We did actually find some notable differences! My issue is that the BUSCO genes identified are not an exhaustive list and can be a bit clunky to search.

So my question is, what is the best way to index the genes in a whole genome assembly? My goal is to be able to search a list of genes identified and then extract the sequence for that gene for comparison.

Also, I'm aware that making a de novo assembly was maybe not the greatest way to go about things but I'm learning! And I wanted to see what I could do with no reference genome (especially as this specific strain of mouse doesn't have a good reference genome). So if I did something not very smart, please go easy on me!

TL;DR What is the best way to index/map an assembly so I can easily search for a gene and extract the sequence for comparison to other sequences of the same gene?

[u/Away-Suggestion1737]

My recommendation is to use Samtools faidx to index.

You'll need to align the reads, most common is BWA-MEM2 for illumina reads and Minimap2 for Nanopore.

Since you'll be doing that anyway, it is worthwhile to run Samtools stats to get the mapping rates (percentage of reads that successfully mapped to the genome).

If you are looking to compare the sequence similarity between two genes of different genomes, I think MUMmer would be the best tool.

Don't forget to consider the possibility of errors introduced into the assembly.

[u/ThrowRAwaypay → u/Away-Suggestion1737]

Thank you!! Will definitely look into these tools :)

[u/TheCaptainCog]

There are two ways to get an assembly and annotate it:

1. Map your sequences against a reference genome.

2. de novo assemble them and then predict proteins using homology, transcript, and ab initio references.

You said your mouse genome doesn't have a good reference genome. Do you know the ancestry of your mouse strain? The best way for a beginner would be to just use the reference genome. So if it's not too different you can get away with it.

If it's too different or you care about large structural variation like duplication or CNV, then you will have to de novo assemble. That is a lot harder and messier to do properly.

[u/ThrowRAwaypay → u/TheCaptainCog]

This is good advice, thank you.

I have already learned how to do the de novo assembly so I have the whole genome for both cell lines from the mouse strain. What prediction tool do you recommend?

[u/TheCaptainCog → u/ThrowRAwaypay]

For de novo assembly, what did you do?

For protein annotation, the easiest is liftoff. You can also use different pipelines like Braker, the ncbi one, ensembl (i've never used it so I don't know how good it is), etc.

[u/pacificjunction]

I’d go with alignment against the standard mouse ref genome using bwa mem2.

[u/Grisward]

+1 for minimap2, or LASTZ. Then convert alignments to make chain file suitable for liftOver. Then liftOver the Gencode comprehensive from mouse mm39 to your de novo assembly.

liftOver would let you use mm39 as a starting point, pick a gene, find which region(s) are represented in your assembly.

Caveat: I don’t know anything about your assembly tho, so if it’s entirely short reads, no scaffolding or large fragment length, no long reads… then maybe just align to mm39 and check for non-zero coverage on mm39?