SEQanswers

Go Back   SEQanswers > Applications Forums > De novo discovery



Similar Threads
Thread Thread Starter Forum Replies Last Post
Why MAQ consensus seq better than SAMtools consensus ?? av_d Genomic Resequencing 5 10-18-2015 03:44 AM
To Build a New Index or Not To Build? arcolombo698 RNA Sequencing 0 11-30-2013 01:04 PM

Reply
 
Thread Tools
Old 03-01-2016, 10:32 AM   #1
sfh838t
Member
 
Location: Mountain Grove, MO, USA

Join Date: Apr 2014
Posts: 23
Default localized consensus build instead of scaffolding

I am trying to make something out of a denovo, heterozygous plant WGS.
it won't go past a contig stage (scaffolding improved average contig length by 1 nt!) for several reasons.
However, the ref seq (species of one of the parents) is covered to 85% with contigs up a few hundred contigs deep in places. I would like to create localized consensuses by collapsing all those reads into one stretch of consensus, without "inserting" the refseq (like GATK alternate ref maker would do). this should be possible by checking bedtools genome coverage for stretches with coverage, isolating all reads for that region from SAM file, then creating a consensus for that stretch, essentially making a scaffold. This would make it easier to make some use of the assembly attempt, for instance for blasting against to find certain genes or promotor regions etc. however for an entire plant genome even I will run out of patience trying to do this by hand!!

Anybody good with python or other languages that could write a script? like I said, it should be possible by finding uninterrupted coverage stretches from bedtools genome coverage, grouping reads with their alignment info from SAM files and then using something like GATK alternate reference maker on each stretch.
Maybe all that's needed would be an executable Linux wrapper script?
Any suggestions? maybe there is a tool already out there that I missed stumbling upon?
sfh838t is offline   Reply With Quote
Old 03-16-2016, 06:40 AM   #2
Markiyan
Member
 
Location: Cambridge

Join Date: Sep 2010
Posts: 85
Default Was it a hiseq 2x100 bp lane run?

So what sequencing technology had been used and total volume of data after adapter/quality trimming?

Assuming it is some sort of the illumina run and no much improvement after scaffolding, I assume it is truseq or nextera shotgun library. Ideally it should have been 2x250 on hiseq or a few miseq runs (2x300), but in quite a few cases people just do 2x100 or 2x125 cycle hiseq, and then wonder what to do with millions of contigs with N50<1kb...

If you are really keen of a few particular genes -

I would start with making legacy blast database from the preprocessed raw reads, aiming for the coverage of 25-50x, (divide and manage it with *.nal files, if formatdb volume count is over 99), stick it on NVMe SSD :-) and try some blasts of 1-2 kb probe dna or aa sequences (use megablast -n T option and up the wordsize to 32 - 64 bp...)

If longer probe sequence - then divide onto 1-2kb chunks and conquer.

If the range is still to short - try nextera matepair, Pacbio, Molecullo, 10K, BAC libraries.

If you wish than you can pull out the matching reads and try assembling them with any assembler (even PHRAP works ok on up to 100K illumina reads).
Markiyan is offline   Reply With Quote
Reply

Tags
consensus sequence, partial chromsome, scaffold

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:44 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO