SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   General (http://seqanswers.com/forums/forumdisplay.php?f=16)
-   -   reads align but contigs do not?? (http://seqanswers.com/forums/showthread.php?t=74115)

sfh838t 02-08-2017 02:18 PM

reads align but contigs do not??
 
Here is something that puzzles me:
I have 4000 reads that align to my reference genome. I can visualize them, they cover maybe 80 % or a little more of the reference sequence.
I assemble these reads into contigs, get something like 160 contigs. ok.
however, when I try to align these contigs, only 27 of them still match my reference.
huh? how is that possible?
Any thoughts?

GenoMax 02-08-2017 02:55 PM

Have you checked to see how much duplication you have in your data? 4000 data tends to have more duplicates if you don't have the right library insert size. Use clumpify.sh from BBMap to get an estimate of total duplicates and optical duplicates. You may not be getting a good assembly if there are too many duplicates.

Brian Bushnell 02-08-2017 06:59 PM

Additionally, lengths and coverage are important here. Possibly, 99% of the bases map to the 27 contigs that match the reference, and the other contigs are short, degenerate contigs from errors or chimerism. Or maybe you're assembling various adapter sequences (another form of chimerism, in a way). It depends on how you're doing the alignment and assembly, and what kind of data you have. Can you describe the situation in more detail? Also, have you BLASTed the unaligned sequences, or compared them to known adapter sequences?

sfh838t 02-09-2017 06:10 AM

these are single end reads from an RNA sequencing project. adapter sequences have been removed.
out of some 37 mil reads I get 4000 reads reads aligning to a just under 8000bp long virus that we know is present in the sample. But these are not siRNA in 20nt length range, these reads range from about 30 to 100 nt long.
I am sure there is duplication, though I do not know how to find or eliminate that and would appreciate any hints as to where to read about this or find tools to work with this.
out of about 160 contigs, 40 are very short, the same length as the shortest reads in fact.
the reads matching the virus (bwa aligned, IGV visualized) seem to visually cover most of the virus sequence, the 27 contigs cover maybe 30 % of the virus sequence.
Thanks for the help !

GenoMax 02-09-2017 07:17 AM

Quote:

Originally Posted by sfh838t (Post 204048)
I am sure there is duplication, though I do not know how to find or eliminate that and would appreciate any hints as to where to read about this or find tools to work with this.

@Brian's post gives directions on how to use clumpify.sh from BBMap to identify, mark/eliminate duplicates. They can be optical or all (PCR+ optical).

sfh838t 02-09-2017 07:35 AM

I did see Brians post, but did not check on that so far. Looks like I will have to take the plunge and install yet another toolbox. Let's see what will go wrong this time.....:)

GenoMax 02-09-2017 07:40 AM

Nothing should go wrong. As long as you have Java 1.7 or greater available nothing else is needed.

sfh838t 02-09-2017 09:19 AM

Ok, I did install and use this clumpify tool.
and it did not make any difference. it found some 300 duplicates, I assemble slightly fewer reads, to similar n50 and total nt count and of the contigs produced I had 1 less contig match than before.....
I have heard the term redundant vs non-redundant. just to make sure I understand: removing these duplicates would make my set of reads non-redundant? Or am I completely lost?

GenoMax 02-09-2017 11:12 AM

What setting did you use for dupedist? For a 4000 flowcell data you need to include
Code:

dedupe=t optical=t dupedist=2500 dupesubs=0 spantiles=f
to your clumpify command line. (dupesubs=0 will only look for perfectly matching reads without any errors).

If you did not use even one of those options then I suggest that you try again.


All times are GMT -8. The time now is 02:29 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.