Using merqury to assess completeness of mobile elements

[EmeraldSama94]

Hi Arang!

I'm trying to assess the completeness of mobile elements in Trypanosoma cruzi genomes, and I thought of using Merqury to do it. Instead of running the read set against the assembly, I'm trying to run the read set and the assembly (separately) against a multi-fasta file containing all the sequences annotated as mobile elements for this parasite. With this, I'm hoping to compare the completeness percentage of the read set with the same value for the assembly, so if they differ it would mean that they are not well represented in the assembly.

To do this, I generated 20-mers of the read library and the assembly using meryl, and then ran Merqury for each, considering the multi-FASTA file as the assembly. I expected that I would get comparable completeness values, but I was met with these results instead:

db_mobile [reads] all 429004 48535497 0.883897
db_mobile [assembly] all 15456 169646 9.11074

I have some questions about these results:

1. I thought that the third column, representing the solid k-mers in the assembly (db_mobile), should be the same regardless of the read set (in this case, the proper reads or the assembly), but it do not seem to be the case. Why is this happening?
2. Do you think that this use of merqury could work? I also thought of generating simulated reads of the assembly, with the same parameters of the original read library, and than running merqury with these two read sets, but I'm kind of stuck :/

Any help is much appreciated!

Best regards,
Samuel

[arangrhie]

Hi Samuel,

That's an interesting usecase of Merqury!
A couple things to consider here -

1. Merqury is collecting completeness stats against what's seen in the 'read meryl db'. Looks like you are querying your TE fasta file against the read and the assembly.
2. You can query the amount of TE kmers found in your assembly against what's 'expected' to be found. This could be either by using 1) the 20-mers collected from your TE fasta files or the 2) TE 20-mers found in your genome (reads).

Collect 20-mers in your TE fasta files

meryl count k=20 TE.fasta output TE.meryl

Query your assembly against that TE.meryl

1) Use the TE.meryl as the 'expected' kmer set
$MERQURY/merqury.sh TE.meryl asm.fasta TE_vs_asm

Query your assembly against the TE 20-mers found in your reads

2) Use the intersect between the TE.meryl and read.meryl as the 'expected' kmer set
meryl intersect reads.meryl TE.meryl output TE_in_reads.meryl
$MERQURY/merqury.sh TE_in_reads.meryl asm.fasta TE_reads_vs_asm

All of 1) and most of the results from 2) will make no sense other than the completeness metric. You might be able to draw the spetra-cn plots out of 2) but I'd expect that to require some manual tweaking after inspecting the .hist file. Once the range has been determined, you could tweak the -m -n parameters in Rscript $MERQURY/plot/plot_spectra_cn.R --help and see if that plots anything interesting.

Bonus, you can compare how many TE 20-mers are found in your reads. You can do the % calculation using the distinct number of kmers in #1 and #2; i.e. 100*#2 / #1.

meryl statistics TE.meryl | head # 1
meryl statistics TE_in_reads.meryl | head # 2

Let me know how this goes, sounds like an exciting experiment. :)

Best,
Arang

[EmeraldSama94]

Hello Arang!

Thank you very much for the inputs and suggestions! I did the second strategy you suggested, querying the assembly against the TE 20-mers found in the reads, as it would give me how much of the genomic information regarding TE is in the assembly, which is my ultimate goal. I used the meryl count-forward command for both the TE multi-fasta and the genomic reads, as the first is comprised of known sequences, making no sense to count k-mers in the reverse strand, and the read library is paired-end, so the reverse counting is already performed with the forward command. Then I ran Merqury, using the assembly as query and the TE_in_reads.meryl as database. I got a 41,59% of completeness, which is kinda cool and could make sense, as these sequences are interspersed in the genome and assemblies using only short reads have a hard time dealing correctly with them.

In order to check if this completeness value is supported by other means, I tried to look at how much of the reads corresponds to TE sequences, and to do this I divided the total (present) number of k-mers in the TE_in_reads.meryl database by the total (present) number of k-mers in the TE.meryl database, and I got a percentage of 4.72%. If I applied the same rationale to the assembly itself, with an asm.meryl and an TE_in_asm.meryl databases, I expected to see something like 1.96% of the assembly assigned as TE, based on the previous completeness value of 41.59% regarding the TE_in_reads database. Instead, I got a value of 6.12%, which I really do not understand. Clearly I'm assuming things wrong, but I do not know what... I know this is unrelated to the software itself, but any lead would be much appreciated. :-)

At last, the spectra-plots really did not make sense lol here is an exemple for your appreciation

Again, thank you very much for all your help!

Best,
Samuel

[arangrhie]

Hi Samuel,

Can you perform the same analysis with meryl count? I am not sure why count-forward is needed in your analysis.. this will distinguish i.e. AAAAGGG and CCCTTTT as different kmers, complicating the analysis. Once they integrate in the genome, TEs are TEs regardless of which strand it got into and the reverse complement should be treated as the same kind. There is always a chance that reads only from one strand got captured... Am I missing something?
I was expecting the % to be much higher than 4.72% - something over 90% since the reads came from the same genome that was assembled (assuming your assembly is close to complete).

If you'd like to trace back where the TE kmer is found in your assembly, try something like
meryl-lookup -bed -kmer TE_in_reads.meryl -sequence asm.fasta | bedtools merge -i - > TE_in_reads.bed
This will light up where your TE kmers are.

Best,
Arang

[EmeraldSama94]

Hello Arang!

I've been running some tests in the last few days, trying to understand the results I'm getting, and now I feel confident discussing them with you. First of all, I clearly misunderstood the count-forward function. I thought that it would only count k-mers in a single direction of a read/contig, and that the normal count would perform the counting in both forward and reverse senses. So, following this logic, if I was using a paired-end read library, I should use count-forward because the library itself should have the forward and reverse segments of a single sequence; if a read library is single-end, then the count method should be applied. But, if I correctly understood your explanation, the difference between count-forward and count is that the latter counts reverse complement k-mers as the same, while the first distinguishes between them, right? I'll stick to the count method nonetheless.

Now, regarding the results: first, I changed the assembly I was using as a test; before, I was using an in-house, in-progress assembly a colleague is working on, and I was not sure if the read set I was using really belonged to it. So, I decided to switch to an already published and considered to be high-quality assembly, so I could have more predictable results. I applied the same strategy you suggested, of querying the assembly against the TE 20-mers found in the reads. For this assembly, the TE completeness was of 93.71%, which made way more sense for a near-complete assembly! In addition, I could also obtain the amount of TE sequences in the read set, by dividing the total k-mer number in the TE_in_reads.meryl by the total k-mers in TE.meryl. With this, I saw that 4.70% of the genome (reads) is comprised of TEs, which also makes sense.

In the previous message, I tried to confirm the Merqury results by comparing the amount of TE sequences in the genome/reads (4.70%) with the amount of TE sequences in the assembly, which is 7.43% in this new dataset. I did not understand then, by now I can see the problem here. The assembly is a representation of an haploid genome, but Trypanosoma cruzi is an aneuploid organism. In terms of assembly, this means that there are less sequences being represented than it should be, and if this is true, the % of TEs regarding the assembly may be over represented if compared to the % of TEs in the genome/reads. So, what I'm trying to do now is to estimate each chromosome copy number by read depth and I'll use this information to artificially duplicate the supernumerary chromosomes, just to represent the correct amount of sequences for the organism. Then, I'll run meryl with this modified assembly and check if the % of TEs now matches the one of the read set. If so, I would have a strong signal that the completeness results for the assembly are correct.

I'm really sorry for the long post and for the broken English, and again I'm beyond then grateful for all the support. As soon as I have the new results, I'll be sure to inform you.

Best regards,
Samuel

[arangrhie]

To make sure I got it right; you did

meryl intersect reads.meryl TE.meryl output TE_in_reads.meryl
$MERQURY/merqury.sh TE_in_reads.meryl asm.fasta TE_reads_vs_asm

And compared the distinct num. of kmers in (1) TE_in_reads.meryl / (2) TE.meryl = 7.43%?
This is still very low, I was expecting something higher than this... Is the repeat database built from a different subspecies?

Regarding your plan for a modified assembly: Creating an additional haplotype by duplicating an existing one wouldn't boost up the completeness score unless it's properly unzipped (phased to contain the haplotype specific variant), it will only increase the number of times a distinct kmer is found in the assembly.
The completeness measure is taking the distinct kmers seen in the assembly, and calculates the portion from the given set (TE_in_reads.meryl in your case).