Hi,

I do not have a SMRT Portal at hand right now, but looking at https://github.com/PacificBiosciences/SMRT-Analysis/wiki/SMRT-Pipe-Reference-Guide-v2.0#-p_mapping-blasr-module suggests that there must be some protocol. The following link also suggests that it is possible to create SAM/BAM-formatted output (https://github.com/PacificBiosciences/cDNA_primer/wiki/A-hands-on-tutorial-of-three-aligners:-BLAT,-BLASR,-and-GMAP#blasr) although it is a command line-based tutorial. However, you might be able to find the corresponding protocols within SMRT Portal.

Shouldn't there be a matching protocol that allows you to select what you called "polished assembly" but only SMRT cells as input, you could follow the command-line steps as BLASR comes with SMRT Portal IIRC.

[EDIT] Actually, there seems to be the possibility to load _external_ references. For a video tutorial have a look at http://files.pacb.com/Training/HowtoLoadaReferenceSequence/story.html . And then you would probably use the P_Mapping protocol/module (s.a. https://speakerdeck.com/pacbio/secondary-analysis slide 8 of 18) in SMRT Portal.

Hope that this helps.

Best,

Cedric

Hi Cedric,

thanks for your suggestions, however these are not simple in a way I meant :-)

Of course mapping reads for polished contigs is possible even with many other software out of SMRT analysis. Command line is also ok. But mapping can lead (we can say "is going to lead") to different results than de novo assembly, because how many algorithms so many results. I was wondering if Quiver /algorithm for polishing which is part of RS_HGAP_Assembly.2 protocol/ can add information about insertions also to a SAM file, not only to contigs which Celera assembler made in previous step of the protocol.

If it's not possible, they've should totally upgrade this in the next version of SMRT analysis.

Hi Karel,

I think I get your point now a bit better. So let me try to rephrase it. You want to use the consensus step which is used in the assembly to report events such as insertions into a separate file (e.g., in SAM format). Is that right?

Best,

Cedric

I don't know if we understand each other. I just want to get correct SAM file with aligned reads for polished assembly (contigs in fasta). In SMRT portal you can only download SAM file with aligned reads for Celera assembly (contigs before polishinq with Quiver).

Just to be able to view it in assembly viewer, e.g. Tablet.

Let me refer then to the fermi paper by Heng Li (http://bioinformatics.oxfordjournals.org/content/28/14/1838) -> "In principle, every analysis based on whole-genome shotgun sequencing (WGS) data, such as SNP and insertion/deletion (INDEL) calling, can also be achieved with unitigs." and "In this sense, any unitig-based analysis has an equivalent read-based analysis, and vice versa.". If I correctly understand his main message in that paper, than he shows that (under certain assumptions) the variant-calling can be done either based on the reads or on the de-novo assembly-based unitigs (say contigs, for the sake of simplicity).

Hence, it should not make a *big* difference if you simply take the polished assembly and map your reads to that. After all, Quiver uses the quality values in addition to the alignment. That means, if Quiver called a certain base that would _not_ have been called by the alignment alone, it is due to the quality value information. This will then be reflected in the post-hoc alignment and could be inferred from the SAM/BAM output of BLASR (e.g., many reads that support an "A" at a certain position but the polished contig contains a "T" at that position).

A note on your comment about "download SAM file with aligned reads for Celera assembly": I assume that this is because HGAP is similar to an Overlap-Layout-Consensus type-of assembler. The Overlap step indeed needs to create the alignments and thus it is "easy" to return the alignments in SAM format. However, Quiver appears not be the "traditional-type-of-aligner" (which it also nowhere claims to be): "Quiver does not use the alignment provided by the mapper (BLASR, typically), except for determining how to group reads together at a macro level. It implicitly performs its own realignment, so it is highly sensitive to all variant types, including indels---for example, it resolves the example above with ease." from (https://github.com/PacificBiosciences/GenomicConsensus/blob/master/doc/QuiverFAQ.rst#what-is-quiver). As such, I would not expect it to report these alignments.

So, in summary, using the polished assembly as a reference (s. fermi above) and then simply aligning your reads to that via SMRT Portal should result in the kind of alignment you seem to be interested in. This alignment should the be viewable in your favorite assembly viewer, e.g., Tablet. Formulated another way around, assume you do a sequencing of an organism with an existing, high-quality genome (which I will equate here to what Quiver returns) and you would simply map your reads to that reference, you would then likely see differences in your reads compared to the reference. Whether these differences come from genuine differences or are errors needs of course to be determined. Coming back to your Quiver-polished assembly, you might also want to look at the assembly and use the GFF file returned by it (https://github.com/PacificBiosciences/GenomicConsensus/blob/master/doc/HowToQuiver.rst#step-3-run-quiver) to see where variants have been called. You would then not even need to perform the alignment, if I get this right.

Best,

Cedric

Hi Cedric,

thank you for your answer and references, I found them very interesting.

Thanks,

Karel