[mira_talk] Re: not=4 does not appear to be working on my Mac OS; system freezes after long read name lengths

  • From: Chris Hoefler <hoeflerb@xxxxxxxxx>
  • To: "mira_talk@xxxxxxxxxxxxx" <mira_talk@xxxxxxxxxxxxx>
  • Date: Mon, 1 Sep 2014 15:50:24 -0500

> I’m not sure what was done to the raw data before we received it. The only
> thing we were told was: 'It looks like your output is around 76 megabases
> once the linker sequences are removed and the forward and reverse strands
> are edited into a high quality consensus.’


It sounds like they ran error-correction on it. Making some assumptions
about your library (10 kb bead-selected) and sequencing chemistry (P4-C2),
you should expect a yield of around 250 Mb per SMRT cell, so 3 cells will
give you around 750 Mb. After error-correction your yield will diminish a
bit, in my experience usually by 4-5X. It sounds like you took a pretty big
hit. I expect that is largely because you don't have enough coverage of the
genome. The error-correction only works efficiently if there is a
sufficient number of overlapping reads. Without that, it won't work well.

Does that refer to using the PacBio long reads as scaffolds to align short
> reads assembled as far as contains to?


Before error-correction methods were available, using the raw PacBio reads
to scaffold contigs was successful to some degree. Basically a modified
BLAST algorithm is used to align the reads to high-confidence contigs,
effectively joining them by some known distance (the length of the read -
the overlaps). The advantages of using PacBio reads to scaffold (as opposed
to a mate-pair library, for example), is the ability to know the precise
size of the gap in bases, and if you have several spanning reads you can
fill the gap with sequence from the reads. However, you still need a decent
coverage for this to work well, and in cases where you have that coverage
you are often better off just error-correcting and assembling de novo
instead. Also, if you have missassemblies in your contigs, scaffolding will
give you a lot of garbage, and you have to be careful with repeats.

Is just trying to get the short reads into even some degree of contigs
> using the long reads feasible?


Given the coverage you have, no, I wouldn't expect anything great. You may
have enough links to support scaffolds in a few places, but it is not going
to help the assembly much.

This actually would be something we’d be interested in. We’re immunologists
> by trade, and are primarily interested in the immune genes. What I’ve done
> so far is to use the 1,000 or so immune-related gene proteins from the
> urchin paper you referenced as tblastn queries (E < 0.001) against our
> IonTorrent database (I removed almost 2/3 of the reads using a filter of
> 90% of bases in a sequence had to have a quality quality cut-off value >=
> Phred=20). We got significant hits from all but a handful of genes


This is a good starting place. Put the known gene sequences that have hits
in your data set in a multi-fasta file and use as the bait. It is probably
best to just do a few at a time. You might try a few different k-mer sizes
depending on how similar the homologs are. The reads that are extracted can
then be used for a de novo assembly of those genes.

Is running mirabait doable on my quad core MacBook Pro with 16 GB memory,


Should be. Well, your Ion data sets total 12 Gb, so make sure you give as
much free memory available to Mira as you can.

I see this caveat in the manual:


Only important if you use -L. So don't use it.

Using MIRA? And could these be done on my Mac?


Yes. The number of reads you extract for each gene will be a small fraction
of your total reads, so it should be completely manageable on your Mac.

Again, I appreciate the depth of your responses.


No problem. Hope it helps.


On Mon, Sep 1, 2014 at 7:45 AM, John DeFilippo <defilippo.john@xxxxxxxxx>
wrote:

> Hi Chris,
>
> Again, I appreciate the depth of your responses.
>
> > you can't do much with PacBio data without error-correction
>
> I’m not sure what was done to the raw data before we received it. The only
> thing we were told was: 'It looks like your output is around 76 megabases
> once the linker sequences are removed and the forward and reverse strands
> are edited into a high quality consensus.’
>
> (For an 800 Mb genome, to me this only looks like 0.1x coverage - what am
> I missing?)
>
> > (scaffolding might be possible in some cases, but it is rarely the best
> option)
>
> Does that refer to using the PacBio long reads as scaffolds to align short
> reads assembled as far as contains to?
>
> > There is a lot that can be done short of a full genome assembly,
>
> Is just trying to get the short reads into even some degree of contigs
> using the long reads feasible?
>
> > A method that comes up periodically on this list is to fish out reads
> (using mirabait or similar) associated with specific features that you are
> interested in (ex. genes, metabolic pathways, ribosomal sequence, viral
> insertions, etc) and assemble them independently for comparison with a
> reference or whatever else you are trying to do.
>
> This actually would be something we’d be interested in. We’re
> immunologists by trade, and are primarily interested in the immune genes.
> What I’ve done so far is to use the 1,000 or so immune-related gene
> proteins from the urchin paper you referenced as tblastn queries (E <
> 0.001) against our IonTorrent database (I removed almost 2/3 of the reads
> using a filter of 90% of bases in a sequence had to have a quality quality
> cut-off value >= Phred=20). We got significant hits from all but a handful
> of genes.
>
> Is running mirabait doable on my quad core MacBook Pro with 16 GB memory,
> or would have to use our HPCC?
>
> I see this caveat in the manual: ‘Do not compute hash statistics from a
> file with sequences, but instead treat the baitfilename as file name of a
> valid mirabait hash statistics file and load it from disk… there are
> currently no fool-guards implemented. This means that the user must
> absolutely make sure to use the same mirabait value for 'k' both in the run
> which generated the hash statistics file and in the search using the
> pre-computed file or else results will be (horribly) wrong.’
>
> Does the fact that I have no idea what this means doom me?
>
> > and assemble them independently for comparison with a reference or
> whatever else you are trying to do.
>
> Using MIRA? And could these be done on my Mac?
>
> Thanks again for your help.
>
> John
> --
> You have received this mail because you are subscribed to the mira_talk
> mailing list. For information on how to subscribe or unsubscribe, please
> visit http://www.chevreux.org/mira_mailinglists.html
>



-- 
Chris Hoefler, PhD
Postdoctoral Research Associate
Straight Lab
Texas A&M University
2128 TAMU
College Station, TX 77843-2128

Other related posts: