Hi all,
I have performed a de novo assembly and have a set of contigs for a specific virus. There is currently not a great reference genome but the beginning and end of the genome are very highly conserved. So, I am wondering if there is an easy way to identify the ends of the genome in my contigs and trim off extra reads that got assembled off of the ends?
1 answer
Before trimming anything, worth looking at coverage across that last 100 bp. Assembler tips run out of read support and stitch in noise, so if the tail sits at 5x while the body is at 500x, that's your answer right there. Flat coverage all the way to the end points more at adapter or primer that survived read trimming, and you'd want to fix that upstream rather than clipping it off the contigs.
For the trim itself, anchor on the conserved sequence instead of cutting a fixed length. Grab ~25 bp from the known genome start, find it in each contig allowing a couple of mismatches, and cut everything upstream of the hit - seqkit locate -p SEQ -m 2 gives you the coordinate, seqkit subseq -r does the cut. Or map each contig to whatever reference you do have with minimap2 -c and take the query start/end of the alignment block, which handles both ends in one pass.
Log in to answer this question.
Is it a circular viral genome?
if not, why would you want to trim off the ends?
Not a circular viral genome. My thought is that there may be extra reads assembled on to what I can pretty confidently say is the beginning and end of the genome, as an artifact of my assembler or extra junk that got sequenced. For example, say the start of the genome is highly conserved and begins with ATCGATCG... and my contig instead starts with TTAAGGCCATCGATCG... I'm then hoping to to cut off the TTAAGGCC, as it might make my sample look more divergent than it actually is.
How do you know the edges are errors? If it is linear, you can look at read alignment to make sure that the edge being trimmed is due to miss-assembly, or whether it is supported by a good number of reads.
I have checked using multiple sequence alignment and the samples assembled using de novo assembly tend to have ~100 bp added onto the ends of the genome that do not align to anything else.
It sounds like you are talking about MSA with another assembly?
I am suggesting you align the raw reads back to your de-novo assembly and check whether the reads support this additional 100bp. If the number of reads supporting your assembly drops before this additional 100bp (say at your suspected beginning) then it seems fair to trim.
If you just have one other assembly for comparison (and no molecular data etc), I would not use this along as evidence for the beginning/end of the genome.
What kind of input sample (DNA/RNA) did you actually start from? is it a metagenomic sample (== contains likely more than 1 species) or is it an isolate sample (== contains only 1 species) ?