2FAST2Q: a general-purpose sequence search and counting program for FASTQ files

Installation and code availability

All 2FAST2Q executable files can be downloaded from zenodo: https://zenodo.org/record/5410822. The code, usage instructions, and test datasets are available on GitHub: https://github.com/veeninglab/2FAST2Q. 2FAST2Q is also a Python package, and can be accessed on PyPI: https://pypi.org/project/fast2q/. When using the executable version on MS Windows or MacOS, no further installation is required and a double click on the executable should suffice. For a more in depth description, please see the online tutorial on https://veeninglab.com/2fast2q. 2FAST2Q is fully implemented in Python3.

Usage considerations

All indicated 2FAST2Q running times were performed on a desktop PC with a 12 core 3.7 GHz processor, and 32GB of RAM. However, 2FAST2Q runs on any up-to-date desktop or laptop. When using 2FAST2Q without mismatch search (perfect alignment only), sample processing should be in the order of seconds or minutes (after file decompression). When using the mismatch search, it is possible for 2FAST2Q analysis to take several minutes per sample. When processing more than one sample, 2FAST2Q will automatically parallelize all analyses by distributing each sample per available processor core.

2FAST2Q fast sequence mismatch search function was possible due to the use of Python numpy (Harris et al., 2020) and numba (Lam, Pitrou & Seibert, 2015) modules. An advanced and in-depth tutorial on 2FAST2Q parameters is available on GitHub and PyPI.

2FAST2Q algorithm

When initialized in standard feature count mode, 2FAST2Q will automatically handle all compressed or uncompressed FASTQ files, and create a hash table for all supplied sequence features. 2FAST2Q will then forward all samples for parallel processing, which can be monitored via progress bars (Fig. S1). Each FASTQ file is sequentially read, saving RAM space. The individually loaded reads are submitted to trimming based on the indicated parameters, either using a fixed position, or a dynamic search. The first assumes the presence of a fixed feature length in the same location for all reads. The second requires one or two search sequences. When one sequence (either up or downstream) is provided, 2FAST2Q will search the read until the sequence is found, and return the predetermined sized feature (again, either up or downstream). When two sequences are used, 2FAST2Q will return any feature within the found search sequences. The location and feature length parameter can thus be ignored in this latter scenario. A sequence mismatch search can also be performed.

Following read trimming, the Phred-score corresponding to each nucleotide of the trimmed sequence is considered. If any of the scores is below the indicated parameter threshold, the read is discarded.

If the read passes quality control, alignment against the input features is finally attempted. Depending on the user input, any kind of feature alignment is performed using either mismatch search or not. By default, 2FAST2Q will always first check for a perfect match. Perfect matching uses hashing, directly comparing all features to the read sequence using hashing runtime complexity. When dealing with mismatches, 2FAST2Q will perform sequence search based on a faster custom made search algorithm. At first, all feature/search sequences are converted to their numerical binary form, subsequently reducing them to integer8 format using numpy. Sequence mismatches are counted by tracking the non-zero result positions of subtracting both sequences. 2FAST2Q mismatch search is therefore based on a Hamming distance calculation. As simple numpy constructs, arithmetic operations can be easily processed using the Python Numba module njit decorator. Therefore, all 2FAST2Q search functions are pre-compiled and effectively run at much faster speeds (Fig. S2). All read sequences searches, and features mismatch alignments are performed using this approach, allowing all search operations to run faster than standard Python code. Moreover, reads that fail to safely align, within the given parameters, to any of the provided features, are stored and used for quick hashed based comparison. The same is performed for reads that align with mismatches. By performing the much faster hashed comparison, this feature avoids the slower de novo mismatch search for previously seen same sequence reads. Runtime is thus decreased, paradoxically maintaining sample processing time as file size increase. “Already seen read” hashing is especially useful with datasets comprising dozens of different independent samples from the same sequencing run (see results). In this case, the generated failed/passed read hash tables for each sample are compiled and used as a seed to the next batch of samples. Each new sample thus takes advantage of the already processed reads in a previous sample, avoiding reprocessing the exact same read several times.

A Python dictionary with a class feature count is used to keep track of all found aligned sequences. When no feature file is provided (i.e., when running in “Extractor+Counter” mode), all found read sequences are returned and counted. Each FASTQ file will originate a unique output file. At the end of the analysis, all samples files are compiled into a single file, which can be readily used for downstream applications.

Developing 2FAST2Q

A major goal when doing targeted (amplicon) sequencing is to know the abundance of each target within a sample. To that end, we wrote the Python-based tool called 2FAST2Q (Fig. 1). 2FAST2Q is able to efficiently extract, align, filter, and count DNA sequences from standard FASTQ files in a single step. 2FAST2Q also performs mismatch sequence searching, nucleotide Phred score quality filtering, dynamic sequence search and trimming (including double sequence search), and automatically loads and detects FASTQ (.gz compressed or not) files. The program also exists as an easy-to-use intuitive executable version for MS Windows, macOS, and Linux, requiring no installation. Alternatively, 2FAST2Q is also available as a Python3 package in the PyPI repository, and can be installed with the “pip install fast2q” command. As input, 2FAST2Q requires only a FASTQ file, and, when reference feature sequences exist (i.e., sgRNAs, barcodes), a .csv file with all the lookup DNA sequences. As an output, 2FAST2Q returns an ordered .csv file with all the raw feature counts per condition, as well as quality control statistics (Figs. 1B–1C). 2FAST2Q contrasts with other current methods by being easy to setup and intuitive to use (Fig. 1A), while simultaneously maintaining advanced configuration settings such as efficient mismatched sequence searching, and quality filtering. 2FAST2Q is thus able of going beyond traditional CRISPRi experimental setups, handling any kind of feature extraction, known or unknown, from FASTQ files.

Counting features using 2FAST2Q

An important feature of performing CRISPRi-seq or RB-seq is to obtain reliable counts of each sgRNA or barcode, for any experimental condition. When using 2FAST2Q in “counting mode”, (i.e., for CRISPRi-seq, or sequence barcode counting), it can be used to quickly obtain an absolute feature sequence count from FASTQ files. Moreover, it might also be of interest to extract all features existing before/in-between/after a given sequence. 2FAST2Q has an “extract and count mode” for this occasion, where the program doesn’t require the input of any feature sequences, and will retrieve the count of all found read sequences. In both instances, the program can search for any feature by either specifying a starting read position, or by providing upstream and/or downstream constant search sequences. The feature length must be specified, except in the latter, where variable sized sequences can be retrieved and/or aligned to (Fig. 2).

Benchmarking 2FAST2Q

2FAST2Q was initially benchmarked against a published CRISPRi-seq dataset comprising 479M reads dispersed over 118 FASTQ files (Rousset et al., 2021). In this study, Rousset et al. (2021) examined which genes are essential in Escherichia coli under different environmental conditions using CRISPRi-seq. 2FAST2Q was used to find an alignment and count the occurrence of each feature from a list of 11,629 sgRNAs across all the 118 files. When only considering perfect alignments between a feature and a read, 2FAST2Q was able to output the final compiled sgRNA count table in 7 min on a personal desktop computer (33s per sample distributed over 10 parallel cores). For comparison, the same files and parameters were also input into MAGeCK. Despite its faster individual file processing speed, its lack of inbuilt sample multiprocessing resulted in a total run time of 23 min. Moreover, MAGeCK fails to return an organized file for all combined samples, leaving the user with the individual count files for each sample (118 in this case). MAGeCK also requires explicit indication of all the FASTQ files to be processed, a time consuming step which 2FAST2Q performs automatically, unless indicated. When comparing the read counts returned from both 2FAST2Q and MAGeCK, a perfect correlation (r = 1) was observed for all features (Fig. S3), indicating similar read counting accuracy.

When allowing for one mismatch in the sgRNA search count, the total run time only increased by 2min, to 9min. Under these program conditions, this corresponds to a more than 40x speed improvement over the use of similar purpose standard search functions, such as the Python regex module match function (Python Software Foundation). For mismatch searching, MAGeCK requires the use of Bowtie2, and respective setup, and thus was not used for further benchmarking.

Using the same dataset published by Rousset et al. (2021) as benchmark data, we assessed the impact of different initial 2FAST2Q parameters on both absolute feature counts, and on downstream data analysis. When not using any Phred-score filtering (Q ≥0), and not allowing for any mismatches, we were able to fully recapitulate the reported total read counts/sgRNA for all conditions (Fig. 3E) (Tables S1 and S2). However, high-quality read length has been reported to improve Illumina sequencing results interpretation (Bokulich et al., 2013). We therefore implemented a filtering for nucleotide wise Phred-scores (Q), where all the sequenced nucleotide scores corresponding to the found feature read location are required to be above an indicated threshold. As expected, filtering using Q ≥30, indicating a 0.1% probability of a nucleotide sequencing mistake, lowers the amount of reads/sgRNA. In some cases, by more than 1 order of magnitude (Fig. 3G). However, when considering the millions of reads generated by a typical sequencing experiment, the presence of mismatches in high quality reads is a likely event (any length of 20 nucleotides with Q ≥30 have, at most, a 2% chance of having a mismatch: 0.001 * 20 = 0.02). We therefore implemented feature mismatch search where a read is considered valid if it unambiguously aligns to a single feature for any number of considered mismatches, thus retrieving more high quality reads, especially from lower overall quality sequencing runs. Allowing for mismatches expectedly increased the number of reads/feature (Figs. 3A, 3C, 3I and 4A), without sacrificing total run time (Fig. 4B) (Tables S3 and S4). As an extreme benchmark case, we allowed for the same number of mismatches as the feature length (20 bp) (Figs. 3I–3J). In practice, these parameters mapped any read to its closest feature, meaning the sequence that unambiguously differs the least from the read. This is performed by an inbuilt safety mechanism, where if more than one feature possible matches the read at the lowest amount of allowed mismatches (i.e., 1), the read is always discarded, but otherwise kept. In regards to the Rousset et al. (2021) dataset, which is on average of high quality, these parameters recovered on average 3% more reads/sgRNA (Fig. 4A). However, it is conceivable that the use and outcome of these parameters varies depending on the experimental setup and user requirements, requiring careful consideration before proceeding to downstream data analysis. In here, we report only on the possibilities of 2FAST2Q functionalities.

Higher stringency parameters can aid in biological discovery

We used the Jupyter notebook analysis pipeline published by Rousset et al. (2021) to assay how these different read processing scenarios impact downstream analysis. Using the different read count tables directly outputted by 2FAST2Q, we calculated and compared the median gene scores as defined by Rousset et al. (2021) (essentially, the median of the log₂ fold change for each feature in all experimental replicates) for the LB medium and gut microbiota medium (GMM) conditions. Using more stringent criteria than Rousset et al. (2021) (a gene is considered significant if it has an absolute gene fold change ≥ 4, instead of ≥ 3.5), we compared how different Phred-scores and mismatch filtering criteria influenced downstream analysis, namely how these criteria influence gene score calculations, and thus gene essentiality (Figs. 3B, 3D, 3F, 3H and 3J).

We observed a higher stringency for the 2FAST2Q parameters of one mismatch, and base pair quality filtering of ≥ 30 (Fig. 3B), with fewer genes being considered essential for any given condition with these criteria than with the criteria that recapitulate the published data (0 mismatches allowed, and no Phred-score consideration) (Fig. 3F). As expected, different read filtering criteria resulted in fold change differences, and consequently in differences in the genes considered essential for these conditions. What criteria to use would depend on the specifics of each individual experiment. The default 2FAST2Q parameters uses a Phred-score of ≥ 30, while allowing up to one mismatch (representing for any 20 basepair (bp) sequence, a 5% bp deviation error with, at maximum, a 2% chance of any nucleotide being wrongly sequenced). As shown in Fig. 3, although the default setting of 2FAST2Q give slightly fewer significant hits, they were all also reported by Rousset et al. (2021). It is also conceivable for a user to be interested in aligning all reads to their closest matching feature. This is possible by setting the total amount of mismatches to the same length of the feature. In this case, the 2FAST2Q inbuilt alignment safety mechanism will prevent ambiguous read alignments from being considered. Once again, we intend only to demonstrate the range of uses of 2FAST2Q. Ultimately, the biological relevance of which parameters to choose is left upon the user.

2FAST2Q dynamically performs FASTQ feature extraction

Under certain experimental setups, the extraction of features from FASTQ files might require the use of a dynamic trimming and search function (i.e., when the location and/or size of the feature differs from read to read) (Fig. 2). In this case, a delimiting search sequence of any length (up and/or downstream of the feature) can be provided. Similar to feature mismatch search, an arbitrary number of mismatches can also be indicated for the search sequence-based trimming, as well as a minimum Phred-score. 2FAST2Q will search each read for the indicated sequences, returning the correctly trimmed read for further processing, and bypassing the need for more complex tools such as Trimmomatic and Bowtie2. As a proof of concept, we used a published CRISPRi-seq dataset by Wei et al. (2021), where dynamic read trimming was required. In this study, a CRISPRi screen was performed using Vero-E6 cells (kidney epithelial cells from an African green monkey) infected with SARS-CoV-2 to identify host genes important for viral replication (Wei et al., 2021). In this dataset, the location of each feature was at a variable location within the read. 2FAST2Q dynamic trimming allowed each read to be independently trimmed based on the relative location of the found search sequences, thus always returning the correct feature location. Using this method, we submitted six FASTQ files (SRR14668185 - SRR14668190) for 2FAST2Q processing. As search sequence we used a 10bp upstream constant sequence (CGAAACACCG), allowing for one mismatch search error in this sequence. We used the provided list of 84,953 sgRNA sequence features, and ran 2FAST2Q (Q ≥30, 0mismatches). 2FAST2Q simultaneously processed all six samples, comprising 324M reads, within 8 min on a standard desktop PC (Tables S5 and S6). This result corresponds to a slowdown of only 22% (speed comparisons were determined using processed reads/second) when compared with the non-dynamic feature extraction process, such as the one we used for the same parameter 2FAST2Q run with the Rousset et al. (2021) dataset.

Recently, Bosch et al. (2021) published a CRISPRi-seq experimental setup with variable length sgRNAs. In this case, both the trimming of each read and the length of each sgRNA need to be considered read by read. This is a feature, to our knowledge, beyond easy implementation in any of the programs mentioned in this work. Once again, 2FAST2Q was also able to extract, count, and align all the found features in a Mycobacterium tuberculosis dataset (SRR13734827), to the provided 96,700 long sgRNA file, albeit using 2 delimiting constant search sequences (upstream: GTACAAAAAC; downstream: TCCCAGATTA), while allowing for one mismatch in each. The returned variable length sequence between the two constant search sequences was used for perfect match alignment against the sgRNAs (Tables S7 and S8). When compared with the non-dynamic extraction process, a slowdown of 44% was observed, in line with what was observed for the Wei et al. dataset.

As the sequence search algorithm uses a similar process to the one used for feature alignment mismatch, a similar speed improvement over standard Python functions is also obtained. Together, these benchmarks demonstrate that 2FAST2Q is a versatile and quick computational tool that can extract relevant features and counts from FASTQ files.