SMRT Analysis: Read filtering

From mn/ibv/bioinfwiki
Jump to: navigation, search


PacBio stores reads and the associated quality values in HDF5-encoded files, which is a binary file format. This means you cannot open these files using a text editor program, unlike for instance fastq or SAM files. If working only with the smrtpipe program, it often may not be necessary to manually inspect the content of HDF5 files. But if using external programs for quality assessment, assembly or other steps, it will become necessary to access these files. This also applies to some specialized usages of the smrtpipe.

Accessing PacBio HDF5 files

PacBio uses two main types of HDF5 files, the bas.h5/bax.h5 files and the cmp.h5 files. Reads are stored in the former format, while alignments (resulting from intermediate steps in the smrtpipe pipeline) are contained in cmp.h5 files.

Prior to the PacBio RS II technology, reads were encoded in one single bas.h5 file. The current version uses three bax.h5 files instead; the (single) bas.h5 still present only holds indices to the reads now contained in the bax.h5 files. The three bax.h5 file names will terminate with “.n.bax.h5”, n giving the part number (i.e. n = 1, 2 or 3). Using PacBio RS II, the bax.h5 files will serve as the input to the smrtpipe program; here, the bas.h5 file can often be ignored.

In the following, only the term bas.h5 will be used, meaning both bas.h5 and bax.h5

Extracting reads from the bas.h5/bax.h5 files

The script, included in the smrtanalysis package, can be used to extract the reads contained in a bas.h5 file:

module load smrtanalysis/2.3.0 --readType subreads --outType fasta /path/to/bas.h5

The resulting file will have the same name as the bas.h5 file, but ending with a “.fasta” extension (in case of using bax.h5 files, the resulting file will NOT contain the bax.h5 part number).

The script can export reads both in the fasta and fastq format. In the latter case, quality values are encoded as phred33 quality scores. Also, three distinct types of PacBio reads can be exported: --readTYPE {ccs, subreads, unrolled}.

  • the “ccs” option exports circular consensus reads, i.e. sequences that have been read multiple times by one sequencing detector.
  • the “subreads” option exports all the subreads present in the raw reads, i.e. the raw reads split by the adaptor sequences.
  • the “unrolled” option exports all raw reads.

For more information about these read types, see the PacBio documentation.

More options for the script are obtained by using the “-h” option, including for instance quality filtering of the exported reads.

Working with cmp.h5 files

The cmp.h5 files contain alignment information produced by certain steps in the smrtpipe program, and serve themselves as input files for subsequent smrtpipe steps. As for the bas.h5 files, a script is available for manipulating cmp.h5 files, included in the smrtanalysis package. Use the “-h” option for further details about this script.

General tools for reading HDF5 files:

The HDF group has made available a number of tools used to manipulate HDF5 files. The HDFView program ( is a Java program providing a graphical user interface to the content of HDF5 files. Additionally, a number of command line tools are available (, most of which are installed on Abel.

For instance, the following command displays the individual sections present in a bas.h5 file:

h5ls -r /path/to/bas.h5

The individual lines of this output can then be used in the h5dump command, so as to read the contends of the referred section. For instance, the following command will display the number of ZWM holes utilized: 

h5dump -y -d  /PulseData/BaseCalls/ZMW/HoleNumber  /payth/to/bas.h5 | head

In the output, the number "nnnnn" in "DATASPACE SIMPLE { ( nnnnn ) / ( H5S_UNLIMITED ) }" indicates the number of holes, i.e. the potential number of reads in the file. Notice, however, that not all holes may produce reads that pass the quality control.

Filtering reads PacBio reads

PacBio quality values

When exporting PacBio reads to a fastq file, the resulting quality values are phred33-encoded quality values. However, these values are calculated from a number of different quality scores encoded within the bas.h5 file. Thus, the precise meaning of the PacBio values may differ from for instance Illumina quality values. For example, the main errors in PacBio sequencing are indels, not incorrect base calls. The inverse of a quality value for a given base gives the probability of this base in reality being something else, including a gap. But it is unclear where to look for the probability that certain base is missing in the PacBio sequence (i.e. a quality value for a gap in the PacBio sequence). When looking at unfiltered PacBio reads, I seeing mean quality scores of around 80%, to be expected for PacBio reads. Thus, the exported phred33 quality values are probably precise enough to allow a meaningful assessment of sequence qualities using tools developed for other sequencing techniques.

FastQC and PacBio reads

The more recent versions of the FastQC program allow the assessment of PacBio read qualities. The FastQC results are presented as for any other sequencing technology. For instance, the display of per-base quality values uses the standard ordering into good (green), mediocre (yellow) and bad (red) categories. I.e. FastQC does not take into account that PacBio reads are by definition of poorer quality than for instance Illumina reads. This of course means that a “fail” status on the per base sequence quality is not meaningful for these reads.

Using the smrtpipe P_Filter module

It is easy to use the filtering module included in the smrtpipe program on its own. Most (maybe all) protocols would include the usage of this module anyway, so an isolated usage of this module makes only sense if subsequent processing is done with tools not included in the smrtpipe. A protocol using the P_Filter module basically stops after the filtering step.

The following protocol produces filtered reads from the bas.h5 files included in the input.xml file:

<?xml version="1.0" encoding="utf-8"?>


  <protocol id="RS_HGAP_Assembly.3" version="2.3.0" editable="false">


  <module id="P_Fetch">


  <module id="P_Filter">

   <param name="minSubReadLength">



   <param name="readScore">



   <param name="minLength">




  <module id="P_FilterReports">




The P_Fetch module is always required in smrtpipe protocols; the P_FilterReports module generates a HTML-based report (it can be excluded if not necessary). In this example, the P_Filter module accepts reads with a mean read score of better than 0.80. This corresponds to a phred33 quality score of around 6, as can be verified if running FastQC on the produced output.

For details of the P_Filter module, such as filtering parameters, input and output file names, see the smrtpipe user guide (

The functionality of the P_Filter module is partially overlapping with that of the script. Using identical readscore parameters, these two programs give identical results. But the P_Filter module provides more additional options than the script, such as for instance a defined signal-to-noise ratio.

Filtering reads based on sequence identity (whitelisting)

In addition to excluding reads below a given quality score, often one wants to exclude certain reads based on their sequence. For instance, it may be necessary to filter out bacterial reads from a eukaryotic sample. In this case, the problematic reads may be identified by aligning all reads to a bacterial genome, keeping only the non-aligning reads.

As a first step, the alignment is done using the BLASR program. This program is capable of aligning long PacBio reads to genomic sequences. It is included in the smrtanalysis package, and is also used internally in the smrtpipe, for instance during read correction.

Next, the ids of all aligned reads are recorded (constituting in effect a “blacklist”).

Finally, going through all reads in the bas.h5 input file, all reads not found to be blacklisted are written to a text file. This text file constitutes the whitelist, and must be included in the protocol.xml file (P_Filter module):

<param name="whiteList" label="Read IDs to whitelist">




Using the whitelisting procedure includes using HDF5 commandline tools to extract read ids, and running a custom script for generating the whitelist. See for details.