[sagarutturkar]

HGAP Parameters

Hi rhall,

Thanks for your help. I am working on this single contig assembly and able to run base_modification correctly. I was trying to assemble another small genome with pacbio RS || data with 135x coverage. Through default HGAP method I got 3 contigs. What other HGAP parameters can be tweaked to get better (possibly 1 contig) assembly for these genome?

Also, what are the changes in the SMRTanalysis 2.0.1 version? Are there any changes that might improve the HGAP results from version 2.0.0?

Thanks
Sagar

[rhall]

The base modification pipeline will not break contigs. The number of contigs output should be the same as in the reference. If you use the reference that is one contig, by first downloading it from the HGAP job, then uploading it as a reference you should only get results for one contig.

[rhall]

Great!
The output from SMRT Portal HGAP goes through correction using the quiver algorithm, with >40x coverage you can expect a base QV of > 50. For the coverage vs QV plot see:

https://speakerdeck.com/pacbio/smrt-sequencing-overview

My normal workflow, given a single contig.
1. Check the degenerate file from the Celera Assembler stage (downloadable from the data section of SMRT Portal). It may contain a few to no sequences, blasting them is the quickest way to see if they are junk, or meaningful.. The reason I check this file is that Celera Assembler assumes even coverage, so any high copy number elements (in particular plasmids) will end up in this file not the final assembly.
2. Check the single contig for circularity using a dot plot:

http://www.mybiosoftware.com/sequence-analysis/6750

3. If the contig is circular, then circularize it. I found the easiest way to do this is to take the polished contig, arbitrarily introduce a break in the middle, then use the two contigs as input to minimus2:

http://sourceforge.net/apps/mediawiki/amos/index.php?title=Minimus2

which should then output a single circular contig by overlapping the start of one contig with the end of the other.
4. Run the circularized contig through RAST, BaSys, or using the original data and the contig as a reference run a Base Modification analysis in SMRT Portal.

http://basys.ca/

http://rast.nmpdr.org/

https://speakerdeck.com/pacbio/base-modification-overview

5. Publish.

[sagarutturkar]

Hi rhall,

Thanks for your reply. I am able to import the data into SMRTportal.

The RS || platform data looks excellent with average read_length of 6KB and longest reads of 25KB. I was able to assemble genome in single contig with HGAP method. This is amazing.

I also tried "RS_modification and motif analysis" protocol which generated consensus of 33 contigs. Does this protocol breaks the assembly at low confidence bases? Is it necessary to do the base-modification or I should go ahead with 1 contig assembly as high-confidence one?

Any suggestions about downstream analysis or any SMRTportal tools which will help me?

[rhall]

and.. The *.metadata.xml is a requirement if the data is to be used in SMRT Portal.

[rhall]

It looks like the files have been moved around and renamed, the simplest solution is probably to ask for the files again from the service provider in the original state.

https://github.com/PacificBiosciences/SMRT-Analysis/wiki/Cannot-import-new-SMRT-Cells

As a general rule, raw data from the instrument should not be moved from the generated directory structure and .h5 files should never be renamed.
Richard.

[sagarutturkar]

PacBio RS|| data

Hi,

We recently got some Data from PacBio RS|| system which looks like this.

10061/
10061/PACBIO_DATA/
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.ccs_stats.txt
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.ccs.fastq
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.1.1.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02_stats.txt
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.1.3.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.1.bas.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02_stats.txt
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.1.1.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.ccs.fastq
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.1.2.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.1.3.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.1.bas.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.1.2.bax.h5
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_D02.fastq
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.ccs_stats.txt
10061/PACBIO_DATA/XSBOR_20130726_RS42150_PL100036291-1_C02.fastq

However, When I copied same folder in inputs_dropbox folder and tried to import - It says no new SMRT cells to import. What might be the problem? After importing I changed the folder name. Also, I do not see the

*.metadata.xml
./Analysis_Results

folders as before. Any suggestions?

[rhall]

Some recent test data is available from:

https://github.com/PacificBiosciences/DevNet/wiki/E.-coli-20kb-Size-Selected-Library-with-P4-C2

This is Ecoli data, and has been used to demonstrate HGAP assembly. You can download the data, import it into SMRT Portal and used it as a test case for HGAP, Base Modification and resequencing pipelines.

[GenoMax]

GenoMax:
For 1.4, you had SMRT portal working, but not command line smrtpipe.py? Given that SMRT portal basically generates parameter files then runs smrtpipe.py this would seem like it should have been a simple fix.

One would think so but the command line usage is disconnected from the SMRTportal in ways we do not fully understand (and have run into for the last 2-3 software revisions).

We use SGE on one cluster. smrtpipe.py when submitted from a qsub script job generates sub-"qsub" jobs on its own .. which SGE does not like.

If you want to know more we can confer off-line as I do not want to hijack this thread.

[rhall]

GenoMax:
Unfortunately the reference guide for 1.4 is not as easy to link to (pdf file). That protocol in particular has not changed, and the documentation will be the same.
For 1.4, you had SMRT portal working, but not command line smrtpipe.py? Given that SMRT portal basically generates parameter files then runs smrtpipe.py this would seem like it should have been a simple fix.

[GenoMax]

rhall: You have linked the reference guide for smrtanalysis v.2.0. Is that applicable for v.1.4 (in case Sagar is still using that version).

We could never get the v.1.4 to work (smrtpipe.py) on the command line.

We are currently working with tech support to get just test data to work with v.2.0 (HGAP protocol) through SMRTportal. Have not thought about command line yet.

[rhall]

The error is simply that no reads have been found that link contigs. You can try changing the parameters:

https://github.com/PacificBiosciences/SMRT-Analysis/wiki/SMRT-Pipe-Reference-Guide-v2.0#parameters-5

I would also double check the reads, by mapping them to the contigs using a standard resequencing protocol.

[sagarutturkar]

I was able to run the AHA using command line and gave me scaffolds for the genome. I used fastX toolkit to convert filtered_subreads.fastq file to fasta format and used existing assembly as reference.


The command used to run smrtpipe was:

/opt/Pacbio/smartanalysis/analysis/bin/smrtpipe.py --params=params.xml xml:input.xml >& smrtpipe.err

However, when I revisit this for some other genome, I am getting errors again.

[ERROR] 2013-05-13 12:57:15,241 [pbpy.smrtpipe.SmrtPipeContext exit 557] Found error: Failed to find any links. Scaffolding failed.
[ERROR] 2013-05-13 12:57:15,241 [pbpy.smrtpipe.SmrtPipeMain run 648] SmrtExit Failed to find any links. Scaffolding failed.
Traceback (most recent call last):
  File "/data2/smart/smartanalysis/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/smrtpipe/SmrtPipeMain.py", line 612, in run
    self._runModules(sModules)
  File "/data2/smart/smartanalysis/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/smrtpipe/SmrtPipeMain.py", line 346, in _runModules
    instance.run()
  File "/data2/smart/smartanalysis/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/smrtpipe/modules/HybridAssembly.py", line 518, in run
    self._context.exit("Failed to find any links. Scaffolding failed.")
  File "/data2/smart/smartanalysis/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/smrtpipe/SmrtPipeContext.py", line 558, in exit
    raise SmrtExit(msg)
SmrtExit: SmrtExit Failed to find any links. Scaffolding failed.
[ERROR] 2013-05-13 12:57:15,243 [pbpy.smrtpipe.SmrtPipeMain exit 760] SmrtExit Failed to find any links. Scaffolding failed.

I have attached complete log files. Please help with this.

Thanks
Sagar

[sagarutturkar]

AHA commandline solved

For the command line job, the AHA workflow should be ran from a filtered subreads file. You will not get sensible results from a bas.h5 file, although it should work. pls2fasta can be used to convert from bas.h5 to fasta

I was able to run the AHA using command line and gave me scaffolds for the genome. I used fastX toolkit to convert filtered_subreads.fastq file to fasta format and used existing assembly as reference.

The XML files were created as described here:

https://github.com/PacificBiosciences/Bioinformatics-Training/wiki/AHA

The command used to run smrtpipe was:

/opt/Pacbio/smartanalysis/analysis/bin/smrtpipe.py --params=params.xml xml:input.xml >& smrtpipe.err

Thanks to rhall and GenoMax for all the help. I am still trying to figure out running the SMRTportal and re-downloading the bas.h5 data.

Thanks
Sagar

[rhall]

For the command line job, the AHA workflow should be ran from a filtered subreads file. You will not get sensible results from a bas.h5 file, although it should work. pls2fasta can be used to convert from bas.h5 to fasta
The SMRT portal job: It looks like the bas.h5 is corrupt, which would also explain the command line job failure. If this is the case, then pls2fasta will also fail, you can use h5dump to look at the contents of a h5 file.
I have seen cases of corrupt bas.h5 files when they are transferred over ftp incorrectly.