Pacbio Large spike at 4kb subread length

[kehey]

Hi all,

I'm running an HGAP on 450kb genome and it's the first time we get this kind of PacBio data. The assembly isn't very good (many contigs), and I was exploring the input data and noticed a number of things on which I would like your feedback:

• There is a large spike around the 4kb region in the subread length distibution. Are these potentially the control reads? And if so, what would this mean? I don't really find elaborate info on the purpose of these control reads, although I read in a post by rhall that these are not removed for instance in an HGAP for a reason.
• The yield from the pre-assembly is 43% (cutoff is 11kb, pre-assembled N50 is 4629)
• Because of the low pre-assembly yield, I went on to see what could be happening. When running a bridgemapper analysis, I noticed that around 60% of the reads are split, and that the split parts (prolog and epilog) are actually mapping onto the same genomic region as their primary mapping. Is this potentially adaptors missing between the forward and reverse (self-ligation)?
• Other stats such as productivity 1 are also low compared to what I find online as reasonable (16% in my case).

If possible, I would just like to have feedback on these observations.

[GenoMax]

@kehey: What chemistry version did this run use? How many SMRTcells was this data from? 16% productivity is not great but that metric has been very hard to influence intelligently in our hands. Small inserts more reliably give high productivity but that is not what one is looking for with PacBio (not for human samples anyway).

[kehey]

@GenoMax: 1 SMRT cell, P4 polymerase and C2 Chemistry.

[GenoMax]

Have you adjusted HGAP assembly parameter (don't have SMRTportal open now but it is under HGAP) appropriately for the small size of the genome you are working with? Default is set to 6 M.

I don't think PacBio has sold a control sized 4kb for a while (currently there is a 2kb control available). Have you asked your provider if they had added any control to your sample?

[kehey]

Yes, I have. Although to be precise, I'm not using the SMRTportal (not installed, using smrtpipe.py directly instead) but I have set that correctly in the xml file that I use to run it.

For reference: I have downloaded the E. coli dataset provided by PacBio as well and that ran perfectly using that xml (with adjusted genomeSize of course). So I'm quite confident in the HGAP run.

And that would not influence the observations from the raw data I described above. So I'm trying to figure out whether there is reason to believe that the data was 'suboptimal'. Especially the fact that these split reads have their prolog and epilog mapping at the exact same genomic position interval strikes me as odd.

We're going back to the provider for feedback on these things as well, but I wanted to have some more points of view on these things apart from the provider if possible. But could you tell me for instance what the control does and why the reads are kept in an HGAP analysis?

[rhall]

What is your Mb yield? It looks like most of the data is coming from the control, run a simple filtering protocol, which will remove the control reads to see your actual coverage for the genome of interest. The control reads are making up such a large portion of all the reads that they are effecting the calculation of seed length cutoff. Try turning off the auto seed calculation, and setting the seed length to 2000.

+ P4C2 had a 4kb control, so it is almost certainly the cause of the spike.

[kehey]

Thanks rhall, this is extremely usefull information! I will do the filtering run for the control reads (any comment on why these are not removed in an HGAP protocol?) and see what the amount of actual reads is. On filtering control reads, what is the difference between protocols/control/SControl.1.xml and RemoveControlReads.1.xml? I see that the latter is never used in any of the xmls while the former is used in e.g. RS_BridgeMapper.1.xml?

I have put a screenshot of the complete set of pre-assembly stats in attachment.

I had already tried lowering the cut-off manually to 4000, but I now realise that was just not low enough given the 4kb control is causing the issue.

On a different note, I already improved the assembly quite a bit by whitelisting my split reads (see first post where I describe the fact that I have 60% of reads that split in a bridgemapper run, of which 90% have their prolog/epilog mapping to the same genomic position, and on opposite strand). So these long reads that were actually mapping in a split fashion were also messing with the assembly. But as you can see from the attachment 2, the peak of the control reads is still in there, so they might still be causing harm.

Thanks!!

[rhall]

HGAP does have the option to remove the control reads, this is where RemoveControlReads.xml is used. Unfortunately it is not trivial to recreate the settings.xml to do this without first generating it via SMRT Portal. I have attached a valid settings.xml for use with smrtpipe.py, it will need some editing.
I didn't catch that the bridge mapped reads mapped on opposite strands. Close split mappings on the same strand can be explained by reference or read quality and the mappings being broken at low quality sequence. Opposite strand mappings are likely the result of sample prep issues. By any chance has this sample gone through amplification, WGA, MDA?

[kehey]

I've been assembling my own xmls by cutting and pasting the different parts I need from the common/protocols/<subfolder> based on what was needed described in the RS_<protocol> files in common/protocols, so thanks for the file! I can use it to compare and see whether I did something wrong.

We're planning to install the portal as well in order to generate the xmls, it just wasn't done when we first installed all the pacbio tools.

I already ran a run using SControl, both using the control 4kb and 4kb_c2 as references, and both times it returned only one read as being control read. I'll run the RemoveControlReads (slightly different settings than SControl?) version and see what that gives.

I'm not sure whether there was amplification involved. Would you be willing to share exactly how you think that could explain the observation? Then I know how to ask the question to the provider.

In any case, I'm reaching the point where I need to get back to the provider, but this information and tips have really helped me to gain confidence on my case (first time I ever had PacBio data in my hands, and only got them 6 weeks ago so I have little reference on data issues). Thanks a lot!

[rhall]

MDA in particular introduces systematic artifacts, specifically a switch-back palindrome as seen in http://www.biomedcentral.com/1472-6750/7/19, but because PacBio has no problems sequencing through palindromes we see this at a lot higher percentage than in the paper.

[kehey]

For future reference of people finding this topic: after consultancy with the provider, it turns out rhall had it right. The provider has discovered that these inverted repeat sequences are artifacts of a whole genome amplification step.

There are basically two ways to trace these problematic reads: based on a reference (either an external one or the one you build using HGAP) using BridgeMapper or by simply looking for reads that consist of long inverted repeats.

So be very cautious when using whole genome amplification. Best to try an reach enough starting material.