Unconfigured Ad

**GenoMax** · 09-11-2013, 03:16 AM

"Rhall" from PacBio participates in this forum and will likely notice this query but you should open a ticket with PacBio tech support in case you have not done so.

**rhall** · 09-11-2013, 09:48 AM

What is the use case? Given the difficulty in barcoding long, none ccs reads, the usefulness of HGAP for barcoded reads is limited? HGAP in SMRT Analysis requires bas.h5 / bax.h5 input, the rich quality information is needed. It is possible to filter for specific reads - originally implemented for removing contamination:

Build software better, together

https://github.com/PacificBiosciences/Bioinformatics-Training/wiki/HGAP-Whitelisting-Tutorial

GitHub is where people build software. More than 225 million people use GitHub to discover, fork, and contribute to over 800 million projects.

It is also technically possible to use custom region tables for trimming the bas.h5, but unless you are very comfortable with h5 files this may be difficult.
Outside of SMRT Analysis, the development toolkit for HGAP allows much more freedom, and will work from fasta files, but it can be difficult to use.

Build software better, together

https://github.com/PacificBiosciences/HBAR-DTK

GitHub is where people build software. More than 225 million people use GitHub to discover, fork, and contribute to over 800 million projects.

**boetsie** · 09-13-2013, 07:37 AM

The goal was to assemble plasmid sequences. Multiple plasmid-pools were multiplexed on a single SMRT cell. After demultiplexing, we had several fastQ files, each belonging to a different plasmid-pool.

For those interested, the fastQ files could not be assembled with the HGAP software, since HGAP requires the bas.h5 files. Therefore, we had to use the original bas.h5 where all pooled plasmids are stored (and thus not demultiplexed). We had to provide the HGAP protocol a set of read-names (so-called whitelist) to indicate which reads to filter out. See this page about whitelists;

Code:

https://github.com/PacificBiosciences/Bioinformatics-Training/wiki/HGAP-Whitelisting-Tutorial

We had four samples multiplexed on this SMRTcell, so for each sample we made a whitelist containing the read-names of the other three samples. To get the read-names out of the fastQ files, we did;

Code:

awk -F '/' '/^@/ {print substr($1,2)"/"$2}' <FASTQ-file>| uniq > whitelist.txt

Next, we added the whitelist to the HGAP protocol (XML);

Code:

<param name="whiteList" label="Minimum Subread Length">
            <value>/PATH-TO/whitelist.txt</value>
            </param>

So, this will filter out the reads belonging to the different samples, but still the barcode of the to-be-assembled sample is present. Therefore, we added an extra trimming for the barcodes (within modulestage "assembly"):

Code:

<param name="layoutOpts" hidden="false">
      <value> --overlapTolerance 100 --trimHit 50 </value>
</param>

These were the steps we had taken (with the help of PacBio support) to get the results. Hope this helps for anyone in the future.

Regards,
Boetsie

**ebioman** · 12-11-2013, 05:35 AM

Hi I would like to tie in with part of the question.
Similarly like you, I discovered that HGAP2 would not take filtered sequences

One might wonder about the use of such a scenario. But e.g. I want to filter my reads and only assemble the chloroplastic genome. So, now I have to take ALL of my 30 smrtcells and then filter later ?? --> not very smart, but hmm...

My problem is in particular with the whitelist. Whenever I launch it, I get errors and it wont succeed. In the end I followed the tutorial step by step, but was still unlucky.

@boetsie: Could you please give me a little "head" of how your data looked like ? I am desperately looking for my mistake.

Error:

Code:

 smrtpipe.py --params=HGAP2.xml xml:input.xml
Traceback (most recent call last):
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/smrtpipe/engine/SmrtPipeWorkflow.py", line 609, in execute
    raise WorkflowError(task.error)
WorkflowError: task://jobId/moduleName/subreadRpt Failed

Whitelist:

Code:

m130905_150152_42182_c100546852550000001823086111241313_s1_p0/0
m130905_150152_42182_c100546852550000001823086111241313_s1_p0/1
m130905_150152_42182_c100546852550000001823086111241313_s1_p0/2
m130905_150152_42182_c100546852550000001823086111241313_s1_p0/3

Snippet of settings:

Code:

    <moduleStage name="filtering" editable="true">
        <module label="PreAssembler Filter v1" id="P_Filter" editableInJob="true">
            <description>Filter reads for use in the pre-assembly step of HGAP, the hierarchical genome assembly process.</description>
            <param name="minSubReadLength" label="Minimum Subread Length">
                <title>The minimum subread length. Shorter subreads will be filtered and excluded from further analysis.</title>
                <value>500</value>
                <input type="text" size="3"/>
                <rule type="number" min="0.0" message="Value must be a positive integer"/>
            </param>
            <param name="readScore" label="Minimum Polymerase Read Quality">
                <title>The minimum polymerase read quality determines the quality cutoff. Polymerase reads with lower quality will be filtered and excluded from further analysis.</title>
                <value>0.80</value>
                <input type="text" size="3"/>
                <rule type="number" min="0.0" message="Value must be between 0 and 1" max="1.0"/>
            </param>
            <param name="minLength" label="Minimum Polymerase Read Length">
                <title>The minimum polymerase read length. Shorter polymerase reads will be excluded from further analysis.</title>
                <value>500</value>
                <input type="text" size="3"/>
                <rule type="number" min="0.0" message="Value must be a positive integer"/>
            </param>
                <param name="whiteList" label="Minimum Subread Length">
    <value>/local/scratch/eschmid/HGAP_chloro/whitelist.txt</value>
</param>
        </module>
        <module label="PreAssemblerSFilter Reports v1" id="P_FilterReports" editableInJob="false"/>
    </moduleStage>

Remark: If I remove my whitelist, everything works well

**GenoMax** · 12-11-2013, 05:58 AM

This is an obvious question but I will still ask. Is "/local/scratch/eschmid/HGAP_chloro/whitelist.txt" readable by "smrtanalysis" user that runs the pipeline?

**ebioman** · 12-11-2013, 06:03 AM

Hi GenoMax
Thanks for the reply.
Yes, since I am running as a local user (eschmid) the pipeline - not in the smartportal but terminal mode. It is actually the current working directory in which I put all the necessary data.

Nevertheless I changed the access of the file to "everybody", to avoid it might complicate things within the pipeline - no change

**boetsie** · 12-11-2013, 06:24 AM

Hi ebioman,

My whitelist looks like:

HTML Code:

m130906_063502_42179_c100578172550000001823093601191416_s1_p0/513
m130906_063502_42179_c100578172550000001823093601191416_s1_p0/1219
m130906_063502_42179_c100578172550000001823093601191416_s1_p0/2063
m130906_063502_42179_c100578172550000001823093601191416_s1_p0/2266
m130906_063502_42179_c100578172550000001823093601191416_s1_p0/2572
m130906_063502_42179_c100578172550000001823093601191416_s1_p0/3088

I used the attached .xml (I had to rename it to .txt b/c seqanswers could not upload XML files). The command I used is;

HTML Code:

smrtpipe.py --params=HGAP_Assembly.txt xml:input.xml

Input.xml looks something like;

HTML Code:

<?xml version="1.0"?>
<pacbioAnalysisInputs>
  <dataReferences>
    <url ref="run:0000000-0002"><location>/home/boetsie/files/m130906.1.bax.h5</location></url>
    <url ref="run:0000000-0002"><location>/home/boetsie/files/m130906.2.bax.h5</location></url>
    <url ref="run:0000000-0002"><location>/home/boetsie/files/m130906.3.bax.h5</location></url>
  </dataReferences>
</pacbioAnalysisInputs>

Hope this helps.
Boetsie

Attached Files

HGAP_Assembly.txt (16.5 KB, 181 views)

**ebioman** · 12-11-2013, 06:36 AM

Thanks a lot boetsie
Indeed your whitelist looks very similar, as does the config and start-command.

I tried your config file which worked, but stopped at the same position and gave me again an error.
I think there might be really some problem with my whitelist - maybe it does not leave enough long reads or something similar. I will try to move from my minimal example to a bigger set to see whether it is a scaling problem.

**rhall** · 12-11-2013, 08:40 AM

It does not look like the filtering is failing, but rather that the filtering step does not result in any data.
task://jobId/moduleName/subreadRpt Failed
Is post whitelisting and filtering.
Does ./data/filtered_subreads.fasta file have any contents?
You can check specific errors in the ./log/P_Filter directory, subreadRpt.log should have a more verbose error.

**ebioman** · 12-12-2013, 12:58 AM

Hello rhall
Indeed I did not realize before that the filtered_subreads.fasta file was empty.

It is strange I broke down everything to a more minimal style

Code:

<?xml version="1.0"?>
<pacbioAnalysisInputs>
  <dataReferences>
    <url ref="run:0000000-0000"><location>/home/eschmid/Napoleom/raw/PacBio/0xJIqqkV_m130905_172422_42182_c100546852550000001823086111241314_s1_p0.2.bax.h5</location></url>
  </dataReferences>
</pacbioAnalysisInputs>

So, only one cell and the whitelist with only one entry from this cell:

Code:

m130905_172422_42182_c100546852550000001823086111241314_s1_p0/100233

Then I checked the filtered_subreads_summary.csv for the particular ZMV and it is present:

Code:

m130905_172422_42182_c100546852550000001823086111241314_s1_p0,100233,17452,20647,3195,0
m130905_172422_42182_c100546852550000001823086111241314_s1_p0,100258,1538,2287,749,0
m130905_172422_42182_c100546852550000001823086111241314_s1_p0,100258,2329,4161,1832,0

same for filtered_summary.csv

Code:

0xJIqqkV_m130905_172422_42182_c100546852550000001823086111241314_s1_p0,0xJIqqkV_m130905_172422_42182_c100546852550000001823086111241314_s1_p0/100233,20658,3195,0.7850,1,1,0,0
0xJIqqkV_m130905_172422_42182_c100546852550000001823086111241314_s1_p0,0xJIqqkV_m130905_172422_42182_c100546852550000001823086111241314_s1_p0/100234,863,0,0.0000,1,0,0,0

But then the log file for the filter Report (/log/P_FilterReports/subreadRPT.log) gives an error (see as well attachment):

Code:

echo 'Finished on' `date -u`;
# Success
exit 0
Running task://Anonymous/P_FilterReports/subreadRpt on Linux dgmserv01.vital-it.ch 2.6.18-348.16.1.el5 #1 SMP Wed Aug 21 04:00:25 EDT 2013 x86_64 x86_64 x86_64 GNU/Linux
Started on Thu Dec 12 08:34:03 UTC 2013
Validating existence of Input Files
Successfully found /local/scratch/eschmid/HGAP_chloro/data/filtered_subread_summary.csv
Successfully validated input files
Traceback (most recent call last):
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/bin/makeFilterSubreadReport.py", line 5, in <module>
    pkg_resources.run_script('pbpy==0.1', 'makeFilterSubreadReport.py')
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/redist/python2.7/lib/python2.7/site-packages/setuptools-0.6c11-py2.7.egg/pkg_resources.py", line 489, in run_script
    keys.append(dist.key)
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/redist/python2.7/lib/python2.7/site-packages/setuptools-0.6c11-py2.7.egg/pkg_resources.py", line 1207, in run_script
    def __init__(self,module):
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/EGG-INFO/scripts/makeFilterSubreadReport.py", line 205, in <module>
    sys.exit(main(sys.argv))
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/EGG-INFO/scripts/makeFilterSubreadReport.py", line 195, in main
    state = app.run()
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/EGG-INFO/scripts/makeFilterSubreadReport.py", line 126, in run
    dpi=self.dpi, useSdf=useSdf)
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/EGG-INFO/scripts/makeFilterSubreadReport.py", line 163, in _makeHistogram
    useSdf=useSdf, useMegabases=True)
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/pbpy-0.1-py2.7.egg/pbpy/plot/PlotHelpers.py", line 460, in makeHistogramWithCumulativePngThemeLpr
    hy, hx, hp = ax.hist(data, bins=60, ec=edgeColor, fc=barcolor)
  File "/software/UHTS/PacBio/SMRTPipe/2.0.1/analysis/lib/python2.7/matplotlib-1.0.1-py2.7-linux-x86_64.egg/matplotlib/axes.py", line 7582, in hist
subreadRpt.log

I really have difficulties to break down the problem here

**rhall** · 12-12-2013, 08:05 AM

If you look at the filtered_subreads.csv file, you need at least one entry with 1,1,1 in the final three columns. If a ZMW is included in the white list it does not necessarily produce sequence, if no sequence is produced the report fails. The last three columns of the filtered_subreads.csv indicate:
1. If the ZMW is a sequencing ZMW - A set number of ZMWs on the chip are non sequencing
2. The productivity of the ZMW (loading) - 0 = no polymerase, 1 = correctly loaded, 2 = multi loaded.
3. If the ZMW passed filter - This included the whitelist, RQ, subread and readlength filters.

**ebioman** · 12-12-2013, 08:25 AM

Thanks for clearing a little bit the fog

I went through the list and found a few reads which hat 0,0,0 but most of them were either
1,1,0 or 1,2,0 (similar for the ones in my example). What bothers me is, that I have no clue why he does not take my whitelist correctly - or, what I thought would be correct.

Since he succeeds if I remove the whitelist, my guess would be that the last ZMV pass filter doesnt fail due to any normal filter but my whitelist.

**rhall** · 12-12-2013, 09:27 AM

Have you tired whitelisting all the ZMWs? 0-150000

**ebioman** · 12-13-2013, 12:47 AM

I think I get slowly closer to the underlying problem rhall.
I checked for the particular smartcell how many ZMWs were present:

Code:

h5ls -r XWLUnE4W_m130905_172422_42182_c100546852550000001823086111241314_s1_p0.3.bax.h5

/PulseData/BaseCalls/ZMW/HoleNumber Dataset {54494/Inf}

Then I generated the whitelist:

Code:

for i in {0..54494} ; do echo m130905_172422_42182_c100546852550000001823086111241314_s1_p0/$i >> whitelist.txt; done

Which resulted in a format such as the one from @boetsie:

Code:

m130905_172422_42182_c100546852550000001823086111241314_s1_p0/0
m130905_172422_42182_c100546852550000001823086111241314_s1_p0/1
m130905_172422_42182_c100546852550000001823086111241314_s1_p0/2
m130905_172422_42182_c100546852550000001823086111241314_s1_p0/3

And again it fails with the same error:
WorkflowError: task://jobId/moduleName/subreadRpt Failed

I tried both, my protocol and the one from @boetsie which both succeed if I remove the whitelist.

The problem is reflected as well in the filtered_summary.csv. If I run it without the whitelist,
most reads get 1,1,1 whereas here again I got 1,1,0.

I have the impression that there is somehow a problem with the format of the whitelist or the position of the "module" in the protocol (which I would dismiss as well since @boetsies protocol produced the same error).

Attached Files

head.whitelist.txt (6.3 KB, 7 views)

Topics	Statistics	Last Post
Genome Sequencing Reveals Large, Repeat-Rich DNA in Trapdoor Spider Named for Stephen Colbert by SEQadmin2 Started by SEQadmin2, Yesterday, 09:57 AM	0 responses 20 views 0 reactions	Last Post by SEQadmin2 Yesterday, 09:57 AM
Single-Cell Sequencing of 5 Million Immune Cells Links Genetic Variants to Disease by SEQadmin2 Started by SEQadmin2, Yesterday, 09:51 AM	0 responses 11 views 0 reactions	Last Post by SEQadmin2 Yesterday, 09:51 AM
Exeter-Led Study Maps What an Average Human Genome Looks Like by SEQadmin2 Started by SEQadmin2, 10-22-2026, 11:13 AM	0 responses 67 views 0 reactions	Last Post by SEQadmin2 10-22-2026, 11:13 AM
Long-Read Sequencing Uncovers Cancer Mutations Conventional Testing Missed by SEQadmin2 Started by SEQadmin2, 10-22-2026, 10:48 AM	0 responses 47 views 0 reactions	Last Post by SEQadmin2 10-22-2026, 10:48 AM

Unconfigured Ad

HGAP with fastQ reads

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Latest Articles

ad_right_rmr

News