Efficient digest of high-throughput sequencing data in a reproducible report

Abstract

Background: High-throughput sequencing (HTS) technologies are spearheading the accelerated development of biomedical research. Processing and summarizing the large amount of data generated by HTS presents a non-trivial challenge to bioinformatics. A commonly adopted standard is to store sequencing reads aligned to a reference genome in SAM (Sequence Alignment/Map) or BAM (Binary Alignment/Map) files. Quality control of SAM/BAM files is a critical checkpoint before downstream analysis. The goal of the current project is to facilitate and standardize this process.

Results: We developed bamchop, a robust program to efficiently summarize key statistical metrics of HTS data stored in BAM files, and to visually present the results in a formatted report. The report documents information about various aspects of HTS data, such as sequencing quality, mapping to a reference genome, sequencing coverage, and base frequency. Bamchop uses the R language and Bioconductor packages to calculate statistical matrices and the Sweave utility and associated LaTeX markup for documentation. Bamchop's efficiency and robustness were tested on BAM files generated by local sequencing facilities and the 1000 Genomes Project. Source code, instruction and example reports of bamchop are freely available from https://github.com/CBMi-BiG/bamchop.

Conclusions: Bamchop enables biomedical researchers to quickly and rigorously evaluate HTS data by providing a convenient synopsis and user-friendly reports.

Figure 1

System architecture of bamchop program.

Figure 2

Estimation of summary statistics by randomly selected sequencing reads. The x-axis indicates the number of reads selected from a BAM file while the y-axis represents the values of eight summary statistics estimated using the selected reads. Each "violin" in the plots represents the distribution of estimated statistics based on 100 resamplings while the horizontal lines correspond to the global averages.

Figure 3

This graphic index represents the sequencing depth along chromosomes. It can be used to quickly identify large regions with extraordinarily high or low depth.

Figure 4

Sequencing depth of an exome sequencing sample. (A) The number and percentage of genomic locations with sequencing depth exceeding given values. (B) Mean sequencing depth of different genomic features. As expected for this exome sample, the exons have much higher mean depth than the other regions.

Figure 5

Sequencing quality. (A) Global distribution of single-base quality scores. An overall high quality of HTS data is suggested as most of the bases were scored around 30 (p = 0.001). (B) Heat map of sequencing quality as a function of read position. The blue solid line and black dashed line indicate the position-specific medians and means of scores, respectively. This figure suggests that the overall sequencing quality of this sample was substantially reduced after 80 bases.

Figure 6

Read mapping statistics. (A) Map quality thresholds. The majority of reads (~95%) have the best mapping score (mapq = 70) assigned by the alignment program, suggesting high confidence of mapping results. (B) Base-level mismatch information, where M = matched bases; I = inserted bases; D = deleted bases; and S = soft clipping bases due to mismatches. (C) Duplication mapping (multiple reads mapped to the same genomic location). The x-axis represents the number of reads sharing the same mapping locations and the y-axis represents the total number of such locations. (D) Insertion size. When the BAM file includes information about paired-end reads, bamchop also summarizes the distribution of the distance between the mapped locations of the read pairs, which is known as insertion size. Insertion size equals to the size of a DNA fragment in sequencing library to be sequenced in pair.

Figure 7

Base frequency. (A) Frequency of N (uncalled) bases due to low quality or ambiguity. (B) Expected versus observed base frequency. (C) Distribution of per-read GC percentages. (D) Position-specific base frequencies at the beginning of reads, which shows a bias in favor of bases G and A to start the sequencing with.

Figure 8

The runtime of bamchop depends on the total number of mapped reads in each BAM file. Diamonds represent the BAM files described in Table 1. The basic runtime of bamchop is about 11 minutes and each 100 million extra reads requires about 10.5 more minute to finish.