diff --git a/README.md b/README.md
index 5bd9ac1e7d0ddb47f1b98dbea90362add92683a2..a26a27c5ea5154b2996cdd54d548aa1ffc011d3e 100644
--- a/README.md
+++ b/README.md
@@ -3,17 +3,21 @@
_fastq_info_ is a command line program written in [Bash](https://www.gnu.org/software/bash/)/[AWK](https://en.wikipedia.org/wiki/AWK) for quickly estimating several standard descriptive statistics from FASTQ-formatted High-Throughput Sequencing (HTS) read files.
Estimated statistics per FASTQ file are:
-  ▹   HTS read and base numbers,
+  ▹   numbers of HTS reads and bases,
-  ▹   HTS read length distribution,
+  ▹   distribution of HTS read lengths,
  ▹   GC-content per HTS read position,
-  ▹   Phred score distribution (global and for each HTS read position),
+  ▹   distribution of Phred scores (Ewing and Green 1998), and corresponding quartiles,
-  ▹   average Phred score per HTS read distribution.
+  ▹   distribution of Phred scores per HTS read position, and corresponding quartiles,
-Different file compression formats can be handled.
+  ▹   distribution of the average Phred score per HTS read, and corresponding quartiles,
+
+  ▹   distribution of the number of sequencing error(s) per HTS read (Edgar and Flyvbjerg 2015), and corresponding quartiles.
+
+Different file compression formats can be handled (e.g. gzip, bzip2, DSRC 2.0).
Several output result formats are available (e.g. reduced/full tables, tab-delimited).
## Installation and execution
@@ -42,7 +46,7 @@ and run it with the following command line model:
_fastq_info_ requires an AWK interpreter in the `$PATH`, which is always the case for most Linux distributions.
By default, _fastq_info_ first considers [gawk](https://www.gnu.org/software/gawk/) (GNU awk, generally available on recent Linux distributions); otherwise the basic command `awk` in the `$PATH` is used.
However, alternative implementations of AWK can be specified using option `-a` (see Usage section).
-In practice, as two POSIX codes are implemented (i.e. original One-True-AWK and extended GNU AWK), _fastq_info_ is able to detect and deal with almost all AWK interpreters (e.g. [nawk](https://github.com/onetrueawk/), [mawk](https://invisible-island.net/mawk/), [goawk](https://github.com/benhoyt/goawk)).
+In practice, as two different POSIX codes are implemented (i.e. original One-True-AWK and extended GNU AWK), _fastq_info_ is able to detect and deal with almost all AWK interpreters (e.g. [nawk](https://github.com/onetrueawk/), [mawk](https://invisible-island.net/mawk/), [goawk](https://github.com/benhoyt/goawk)).
##### About compressed FASTQ files
@@ -122,9 +126,9 @@ Run _fastq_info_ without option to read the following documentation:
* In outputted results, every empty entry is indicated by a dot instead of zero.
-* Tab-delimited option `-v t` enables to output only several statistics: numbers of HTS reads and bases (NR and NB, respectively), average HTS read length (AL), the three quartiles of the global Phred score distribution (BQ1, BQ2, BQ3) and the three quartiles of the average Phred score per HTS read distribution (RQ1, RQ2, RQ3). For detailed distributions per HTS read position and/or Phred score value, use options `-v r` or `-v f`.
+* Tab-delimited option `-v t` enables to output only several statistics: numbers of HTS reads and bases (NR and NB, respectively), average HTS read length (AL), the three quartiles of the distribution of the Phred scores (BQ1, BQ2, BQ3), the three quartiles of the distribution of the average Phred score per HTS read (RQ1, RQ2, RQ3) and the three quartiles of the distribution of the (expected) number of sequencing error(s) per HTS read (EQ1, EQ2, EQ3). For detailed distributions per HTS read position and/or Phred score value, use options `-v r` or `-v f`. However it is worth noticing that _fastq_info_ runs twice faster when using option `-v t`.
-* Specific AWK interpreters can be used via the option `-a` (either a name within the `$PATH` or the full path to a binary). _fastq_info_ was successfully run together with [gawk](https://www.gnu.org/software/gawk/), [nawk](https://github.com/onetrueawk/), [mawk](https://invisible-island.net/mawk/), and [goawk](https://github.com/benhoyt/goawk). However, faster running times were generally observed using [gawk](https://www.gnu.org/software/gawk/) versions ≥ 4.0 ([mawk](https://invisible-island.net/mawk/) remains acceptable, but [goawk](https://github.com/benhoyt/goawk) is not recommended).
+* Specific AWK interpreters can be used via the option `-a` (to set either a name within the `$PATH` or the full path to a binary). _fastq_info_ was successfully run together with [gawk](https://www.gnu.org/software/gawk/), [nawk](https://github.com/onetrueawk/), [mawk](https://invisible-island.net/mawk/), and [goawk](https://github.com/benhoyt/goawk). However, faster running times were generally observed using [gawk](https://www.gnu.org/software/gawk/) versions ≥ 4.0 ([mawk](https://invisible-island.net/mawk/) remains acceptable, but [goawk](https://github.com/benhoyt/goawk) is not recommended).
* Option `-c` can be useful to obtain a check list of the required/expected binaries available in the `$PATH`, as well as their respective version (especially for the AWK interpreter).
@@ -146,7 +150,7 @@ wget -q ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666*.fastq.gz
The following command line runs `fastq_info.sh` to analyze the second (i.e. R2) downloaded file :
```bash
-./fastq_info.sh SRR001666_2.fastq.gz
+fastq_info.sh SRR001666_2.fastq.gz
```
leading to the following standard output:
@@ -154,166 +158,172 @@ leading to the following standard output:
```
##File: SRR001666_2.fastq.gz
-#no.reads: 7047668
-#no.bases: 253716048
-#avg.lgt: 36.0
-------------------------------
-pos Lfreq GCfreq Q1 Q2 Q3
------ ------ ------ -- -- --
-1 . 47.61 40 40 40
-2 . 47.18 40 40 40
-3 . 49.13 40 40 40
-4 . 50.43 40 40 40
-5 . 49.95 40 40 40
-6 . 51.04 40 40 40
-7 . 50.40 40 40 40
-8 . 50.22 40 40 40
-9 . 50.63 40 40 40
-10 . 50.74 40 40 40
-11 . 50.72 40 40 40
-12 . 51.04 39 40 40
-13 . 50.77 35 40 40
-14 . 50.69 33 40 40
-15 . 51.01 31 40 40
-16 . 51.08 28 40 40
-17 . 51.07 25 40 40
-18 . 51.44 24 40 40
-19 . 51.40 22 39 40
-20 . 51.25 20 35 40
-21 . 51.52 19 33 40
-22 . 51.46 18 31 40
-23 . 51.24 16 28 40
-24 . 51.53 15 26 40
-25 . 51.55 13 24 40
-26 . 51.46 12 22 38
-27 . 51.75 12 21 36
-28 . 51.72 11 19 33
-29 . 51.45 10 18 32
-30 . 51.85 10 17 29
-31 . 51.43 9 16 28
-32 . 51.36 9 15 26
-33 . 51.74 8 14 24
-34 . 51.73 7 13 23
-35 . 51.63 7 12 22
-36 100.00 51.70 6 12 22
------ ------ ------ -- -- --
- Q1 Q2 Q3
- -- -- --
-bases 18 40 40
-reads 26 30 34
-------------------------------
+#no.reads(NR): 7047668
+#no.bases(NB): 253716048
+#avg.lgt(AL): 36.0
+------------------------------------
+n Lfreq GCfreq Efreq Q1 Q2 Q3
+----- ------ ------ ------ -- -- --
+0 . . 47.17 . . .
+1 . 47.61 26.28 40 40 40
+2 . 47.18 12.42 40 40 40
+3 . 49.13 5.97 40 40 40
+4 . 50.43 2.98 40 40 40
+5 . 49.95 1.64 40 40 40
+6 . 51.04 1.00 40 40 40
+7 . 50.40 0.63 40 40 40
+8 . 50.22 0.45 40 40 40
+9 . 50.63 0.32 40 40 40
+10 . 50.74 0.25 40 40 40
+11 . 50.72 0.17 40 40 40
+12 . 51.04 0.13 39 40 40
+13 . 50.77 0.11 35 40 40
+14 . 50.69 0.09 33 40 40
+15 . 51.01 0.06 31 40 40
+16 . 51.08 0.05 28 40 40
+17 . 51.07 0.03 25 40 40
+18 . 51.44 0.02 24 40 40
+19 . 51.40 0.01 22 39 40
+20 . 51.25 0.01 20 35 40
+21 . 51.52 0.01 19 33 40
+22 . 51.46 0.01 18 31 40
+23 . 51.24 0.01 16 28 40
+24 . 51.53 0.01 15 26 40
+25 . 51.55 0.01 13 24 40
+26 . 51.46 0.01 12 22 38
+27 . 51.75 0.01 12 21 36
+28 . 51.72 0.01 11 19 33
+29 . 51.45 0.01 10 18 32
+30 . 51.85 0.01 10 17 29
+31 . 51.43 0.01 9 16 28
+32 . 51.36 0.01 9 15 26
+33 . 51.74 0.01 8 14 24
+34 . 51.73 0.02 7 13 23
+35 . 51.63 0.04 7 12 22
+36 100.00 51.70 0.06 6 12 22
+----- ------ ------ ------ -- -- --
+ Q1 Q2 Q3
+ -- -- --
+Phred.per.base(B) 18 40 40
+avg.Phred.per.read(R) 26 30 34
+error.per.read(E) 0 1 2
+------------------------------------
```
-The first part of the outputted table is made up by one row per HTS read position (column `pos`).
-For each `pos` value (varying from 1 to the largest observed HTS read length), the corresponding row indicates the percentage of HTS read of length being equal to `pos` (column `Lfreq`), the percentage of observed GC bases (column `GCfreq`), and the 1st, 2nd and 3rd quartiles of observed Phred scores (columns `Q1`, `Q2` and `Q3`, respectively).
-The bottom part of the table summarizes the global Phred score distribution (row `bases`: three quartiles `Q1`, `Q2` and `Q3`), and the average Phred score per HTS read distribution (last row `reads`: three quartiles `Q1`, `Q2` and `Q3`).
+The first part of the outputted table is made up by one row per HTS read position (column `n`).
+For each `n` value (varying from 0 to the largest observed HTS read length), the corresponding row indicates the percentage of HTS reads of length being equal to `n` (column `Lfreq`), the percentage of observed GC bases (column `GCfreq`), the percentage of HTS reads with `n` sequencing error(s) (column `Efreq`), and the 1st, 2nd and 3rd quartiles of observed Phred scores at position `n` (columns `Q1`, `Q2` and `Q3`, respectively).
+The bottom part of the table summarizes the global Phred score distribution (first row `Phred.per.base(B)`: three quartiles `Q1`, `Q2` and `Q3`), the average Phred score per HTS read distribution (middle row `avg.Phred.per.read(R)`: three quartiles `Q1`, `Q2` and `Q3`), and the (expected) number of error(s) per HTS read distribution (last row `error.per.read(E)`).
The above example therefore shows that the majority of Phred scores are decreasing below _Q_ = 20 at positions 28-36 (i.e. the median Phred score _Q_<sub>2</sub> is lower than 20 as of HTS read position 28).
-At least 25% of all sequenced bases are associated to Phred scores < 19 (i.e. first quartile _Q_<sub>1</sub> = 18 in row `bases`), but at least 50% of the HTS reads have an average Phred score > 29 (median _Q_<sub>2</sub> = 30 in row `reads`)
+At least 25% of all sequenced bases are associated to Phred scores < 19 (i.e. first quartile _Q_<sub>1</sub> = 18 in row `B`), but at least 50% of the HTS reads have an average Phred score > 29 (median _Q_<sub>2</sub> = 30 in row `R`).
+At least 2 sequencing errors are expected within 25% of the HTS reads (third quartile _Q_<sub>3</sub> = 2 in last row `E`).
#### Advanced usage
-For more details (i.e. one supplementary column for each observed Phred score), a full table can be outputted using option `-v f`:
+For more details (i.e. one supplementary column for each observed Phred score _Q_), a full table can be outputted using option `-v f`:
```bash
-./fastq_info.sh -v f SRR001666_2.fastq.gz
+fastq_info.sh -v f SRR001666_2.fastq.gz
```
```
##File: SRR001666_2.fastq.gz
-#no.reads: 7047668
-#no.bases: 253716048
-#avg.lgt: 36.0
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-pos Lfreq GCfreq Q1 Q2 Q3 Q= 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
------ ------ ------ -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
-1 . 47.61 40 40 40 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 94.3
-2 . 47.18 40 40 40 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 94.2
-3 . 49.13 40 40 40 0.2 0.1 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 92.3
-4 . 50.43 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 92.3
-5 . 49.95 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.3 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 91.5
-6 . 51.04 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.3 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.3 0.4 0.4 0.4 0.4 0.4 90.6
-7 . 50.40 40 40 40 0.2 0.1 0.0 0.1 0.2 0.2 0.1 0.1 0.1 0.2 0.3 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 87.4
-8 . 50.22 40 40 40 0.2 0.1 0.0 0.1 0.2 0.2 0.1 0.1 0.1 0.2 0.4 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 86.8
-9 . 50.63 40 40 40 0.2 0.1 0.0 0.1 0.2 0.3 0.2 0.2 0.2 0.2 0.5 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 82.7
-10 . 50.74 40 40 40 0.2 0.1 0.1 0.1 0.3 0.3 0.2 0.2 0.2 0.2 0.6 0.3 0.3 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 80.6
-11 . 50.72 40 40 40 0.3 0.1 0.1 0.2 0.3 0.3 0.2 0.2 0.2 0.3 0.6 0.3 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 78.6
-12 . 51.04 39 40 40 0.2 0.1 0.1 0.2 0.3 0.4 0.2 0.3 0.3 0.3 0.7 0.4 0.5 0.5 0.5 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 1.0 0.9 1.0 1.0 74.4
-13 . 50.77 35 40 40 0.2 0.2 0.1 0.2 0.3 0.4 0.2 0.3 0.3 0.4 0.8 0.5 0.5 0.6 0.6 0.7 0.7 0.8 0.8 0.8 0.8 0.9 0.9 0.9 1.0 1.0 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 70.1
-14 . 50.69 33 40 40 0.2 0.2 0.1 0.2 0.4 0.5 0.3 0.3 0.3 0.4 0.9 0.5 0.6 0.6 0.7 0.7 0.8 0.8 0.8 0.9 0.9 1.0 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.2 1.1 1.2 1.1 1.2 1.1 1.1 1.1 67.8
-15 . 51.01 31 40 40 0.2 0.2 0.1 0.2 0.4 0.5 0.3 0.3 0.4 0.5 1.1 0.6 0.7 0.7 0.8 0.8 0.9 0.9 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 64.9
-16 . 51.08 28 40 40 0.2 0.3 0.1 0.3 0.5 0.6 0.4 0.5 0.5 0.6 1.3 0.8 0.8 0.9 0.9 1.0 1.1 1.1 1.1 1.2 1.2 1.2 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.4 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.2 1.2 59.9
-17 . 51.07 25 40 40 0.2 0.3 0.2 0.3 0.6 0.8 0.5 0.6 0.6 0.7 1.6 0.9 1.0 1.1 1.1 1.2 1.2 1.3 1.3 1.4 1.4 1.4 1.4 1.5 1.5 1.5 1.4 1.5 1.5 1.5 1.4 1.5 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.3 54.9
-18 . 51.44 24 40 40 0.4 0.3 0.2 0.4 0.6 0.8 0.5 0.6 0.7 0.7 1.7 1.0 1.1 1.1 1.2 1.3 1.3 1.4 1.4 1.5 1.5 1.5 1.6 1.5 1.5 1.6 1.6 1.6 1.5 1.6 1.6 1.5 1.5 1.5 1.4 1.4 1.4 1.4 1.3 1.3 52.2
-19 . 51.40 22 39 40 0.2 0.3 0.2 0.4 0.7 0.9 0.6 0.7 0.7 0.9 2.0 1.1 1.2 1.3 1.4 1.4 1.5 1.5 1.6 1.6 1.6 1.7 1.6 1.7 1.6 1.7 1.7 1.6 1.6 1.6 1.6 1.6 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 48.8
-20 . 51.25 20 35 40 0.2 0.4 0.2 0.5 0.9 1.1 0.7 0.8 0.9 1.0 2.4 1.3 1.4 1.5 1.6 1.6 1.7 1.7 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.7 1.7 1.6 1.6 1.6 1.5 1.5 1.4 1.4 1.4 1.3 1.3 43.9
-21 . 51.52 19 33 40 0.3 0.4 0.3 0.6 1.0 1.3 0.8 0.9 1.0 1.1 2.6 1.5 1.6 1.6 1.8 1.8 1.8 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.8 1.8 1.8 1.8 1.7 1.7 1.6 1.6 1.5 1.5 1.4 1.4 1.4 1.3 1.3 40.8
-22 . 51.46 18 31 40 0.2 0.5 0.3 0.7 1.1 1.4 0.9 1.0 1.2 1.3 3.0 1.7 1.7 1.8 1.9 2.0 2.0 2.0 2.1 2.0 2.1 2.1 2.0 2.0 2.0 1.9 1.9 1.9 1.8 1.8 1.7 1.6 1.6 1.5 1.5 1.4 1.4 1.3 1.3 1.2 37.3
-23 . 51.24 16 28 40 0.3 0.6 0.4 0.8 1.2 1.7 1.0 1.2 1.3 1.5 3.4 1.9 2.0 2.1 2.1 2.2 2.2 2.3 2.3 2.2 2.2 2.2 2.1 2.1 2.1 2.0 1.9 1.8 1.8 1.7 1.7 1.6 1.6 1.5 1.4 1.4 1.3 1.3 1.2 1.2 33.3
-24 . 51.53 15 26 40 0.3 0.7 0.4 0.9 1.4 1.9 1.2 1.4 1.6 1.7 3.9 2.1 2.2 2.3 2.4 2.4 2.4 2.4 2.4 2.4 2.3 2.3 2.2 2.2 2.1 2.1 1.9 1.9 1.8 1.7 1.6 1.6 1.5 1.5 1.4 1.3 1.2 1.2 1.2 1.1 29.5
-25 . 51.55 13 24 40 0.2 0.8 0.5 1.0 1.7 2.3 1.4 1.6 1.8 2.0 4.4 2.4 2.5 2.6 2.6 2.6 2.6 2.6 2.5 2.5 2.4 2.3 2.2 2.2 2.1 2.0 1.9 1.8 1.7 1.7 1.6 1.5 1.4 1.4 1.3 1.2 1.2 1.1 1.1 1.0 26.0
-26 . 51.46 12 22 38 0.2 1.0 0.6 1.2 2.0 2.6 1.6 1.8 2.0 2.2 4.9 2.6 2.7 2.7 2.7 2.7 2.7 2.7 2.6 2.5 2.4 2.4 2.2 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.5 1.4 1.3 1.2 1.2 1.1 1.1 1.0 1.0 23.4
-27 . 51.75 12 21 36 0.2 1.2 0.7 1.4 2.2 2.9 1.8 2.0 2.2 2.4 5.2 2.8 2.8 2.9 2.8 2.8 2.8 2.7 2.7 2.6 2.4 2.3 2.3 2.1 2.0 2.0 1.8 1.7 1.6 1.6 1.5 1.4 1.3 1.2 1.2 1.1 1.1 1.0 1.0 0.9 21.3
-28 . 51.72 11 19 33 0.2 1.3 0.8 1.6 2.5 3.3 2.0 2.2 2.4 2.6 5.7 3.0 3.0 3.0 3.0 2.9 2.8 2.7 2.7 2.5 2.4 2.3 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.2 1.1 1.0 1.0 0.9 0.9 0.8 19.1
-29 . 51.45 10 18 32 0.4 1.6 0.9 1.8 2.8 3.6 2.2 2.4 2.6 2.8 6.0 3.1 3.1 3.1 3.1 2.9 2.9 2.8 2.7 2.6 2.4 2.3 2.2 2.1 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.3 1.2 1.1 1.0 1.0 0.9 0.9 0.8 0.8 17.3
-30 . 51.85 10 17 29 0.3 1.8 1.1 2.1 3.2 4.1 2.5 2.7 2.8 3.0 6.4 3.3 3.3 3.2 3.1 3.0 2.9 2.8 2.7 2.5 2.4 2.2 2.1 2.0 1.8 1.8 1.6 1.5 1.4 1.4 1.2 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 15.4
-31 . 51.43 9 16 28 0.2 2.1 1.3 2.4 3.5 4.6 2.6 2.9 3.1 3.2 6.8 3.4 3.4 3.3 3.2 3.1 2.9 2.8 2.6 2.5 2.3 2.2 2.0 1.9 1.8 1.7 1.5 1.4 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 0.6 13.8
-32 . 51.36 9 15 26 0.3 2.5 1.4 2.6 3.9 4.9 2.9 3.1 3.4 3.4 7.1 3.5 3.4 3.3 3.2 3.1 2.9 2.8 2.6 2.4 2.2 2.1 2.0 1.8 1.7 1.6 1.4 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 0.6 0.6 12.7
-33 . 51.74 8 14 24 0.2 2.9 1.7 3.0 4.3 5.5 3.1 3.4 3.6 3.7 7.3 3.6 3.5 3.4 3.2 3.0 2.8 2.7 2.5 2.3 2.2 2.0 1.9 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.9 0.8 0.7 0.7 0.7 0.6 0.6 0.5 11.1
-34 . 51.73 7 13 23 0.2 3.3 1.9 3.3 4.8 5.9 3.4 3.6 3.8 3.9 7.7 3.7 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.1 1.9 1.7 1.6 1.5 1.4 1.2 1.2 1.1 1.0 0.9 0.9 0.8 0.7 0.7 0.7 0.6 0.6 0.5 0.5 9.7
-35 . 51.63 7 12 22 0.3 3.9 2.1 3.7 5.2 6.4 3.5 3.8 3.9 4.0 7.8 3.7 3.5 3.4 3.2 3.0 2.7 2.5 2.3 2.1 2.0 1.8 1.7 1.5 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 8.8
-36 100.00 51.70 6 12 22 0.2 4.1 2.4 4.0 5.5 6.7 3.6 3.8 3.9 3.9 7.6 3.6 3.5 3.3 3.0 2.8 2.6 2.4 2.2 2.1 1.9 1.8 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.5 0.4 9.0
------ ------ ------ -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
- Q1 Q2 Q3 Q= 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
- -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
-bases 18 40 40 0.2 0.9 0.5 1.0 1.5 1.9 1.1 1.2 1.3 1.4 2.9 1.5 1.6 1.6 1.6 1.5 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.2 1.1 1.1 1.1 1.0 1.0 1.0 1.0 0.9 0.9 0.9 0.9 0.8 0.8 51.0
-reads 26 30 34 0.1 0.0 0.0 0.0 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.6 0.8 1.0 1.3 1.6 2.2 2.7 3.5 4.2 5.1 5.8 6.5 7.0 7.3 7.3 7.1 6.7 6.2 5.5 4.8 3.9 2.9 1.9 0.7
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+#no.reads(NR): 7047668
+#no.bases(NB): 253716048
+#avg.lgt(AL): 36.0
+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+n Lfreq GCfreq Efreq Q1 Q2 Q3 Q= 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
+----- ------ ------ ------ -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+0 . . 47.17 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
+1 . 47.61 26.28 40 40 40 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 94.3
+2 . 47.18 12.42 40 40 40 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 94.2
+3 . 49.13 5.97 40 40 40 0.2 0.1 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 92.3
+4 . 50.43 2.98 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 92.3
+5 . 49.95 1.64 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.3 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 91.5
+6 . 51.04 1.00 40 40 40 0.2 0.1 0.0 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.3 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.3 0.4 0.4 0.4 0.4 0.4 90.6
+7 . 50.40 0.63 40 40 40 0.2 0.1 0.0 0.1 0.2 0.2 0.1 0.1 0.1 0.2 0.3 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 87.4
+8 . 50.22 0.45 40 40 40 0.2 0.1 0.0 0.1 0.2 0.2 0.1 0.1 0.1 0.2 0.4 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 86.8
+9 . 50.63 0.32 40 40 40 0.2 0.1 0.0 0.1 0.2 0.3 0.2 0.2 0.2 0.2 0.5 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 82.7
+10 . 50.74 0.25 40 40 40 0.2 0.1 0.1 0.1 0.3 0.3 0.2 0.2 0.2 0.2 0.6 0.3 0.3 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 80.6
+11 . 50.72 0.17 40 40 40 0.3 0.1 0.1 0.2 0.3 0.3 0.2 0.2 0.2 0.3 0.6 0.3 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 78.6
+12 . 51.04 0.13 39 40 40 0.2 0.1 0.1 0.2 0.3 0.4 0.2 0.3 0.3 0.3 0.7 0.4 0.5 0.5 0.5 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 1.0 0.9 1.0 1.0 74.4
+13 . 50.77 0.11 35 40 40 0.2 0.2 0.1 0.2 0.3 0.4 0.2 0.3 0.3 0.4 0.8 0.5 0.5 0.6 0.6 0.7 0.7 0.8 0.8 0.8 0.8 0.9 0.9 0.9 1.0 1.0 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 70.1
+14 . 50.69 0.09 33 40 40 0.2 0.2 0.1 0.2 0.4 0.5 0.3 0.3 0.3 0.4 0.9 0.5 0.6 0.6 0.7 0.7 0.8 0.8 0.8 0.9 0.9 1.0 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.2 1.1 1.2 1.1 1.2 1.1 1.1 1.1 67.8
+15 . 51.01 0.06 31 40 40 0.2 0.2 0.1 0.2 0.4 0.5 0.3 0.3 0.4 0.5 1.1 0.6 0.7 0.7 0.8 0.8 0.9 0.9 1.0 1.0 1.0 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 64.9
+16 . 51.08 0.05 28 40 40 0.2 0.3 0.1 0.3 0.5 0.6 0.4 0.5 0.5 0.6 1.3 0.8 0.8 0.9 0.9 1.0 1.1 1.1 1.1 1.2 1.2 1.2 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.4 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.2 1.2 59.9
+17 . 51.07 0.03 25 40 40 0.2 0.3 0.2 0.3 0.6 0.8 0.5 0.6 0.6 0.7 1.6 0.9 1.0 1.1 1.1 1.2 1.2 1.3 1.3 1.4 1.4 1.4 1.4 1.5 1.5 1.5 1.4 1.5 1.5 1.5 1.4 1.5 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.3 54.9
+18 . 51.44 0.02 24 40 40 0.4 0.3 0.2 0.4 0.6 0.8 0.5 0.6 0.7 0.7 1.7 1.0 1.1 1.1 1.2 1.3 1.3 1.4 1.4 1.5 1.5 1.5 1.6 1.5 1.5 1.6 1.6 1.6 1.5 1.6 1.6 1.5 1.5 1.5 1.4 1.4 1.4 1.4 1.3 1.3 52.2
+19 . 51.40 0.01 22 39 40 0.2 0.3 0.2 0.4 0.7 0.9 0.6 0.7 0.7 0.9 2.0 1.1 1.2 1.3 1.4 1.4 1.5 1.5 1.6 1.6 1.6 1.7 1.6 1.7 1.6 1.7 1.7 1.6 1.6 1.6 1.6 1.6 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 48.8
+20 . 51.25 0.01 20 35 40 0.2 0.4 0.2 0.5 0.9 1.1 0.7 0.8 0.9 1.0 2.4 1.3 1.4 1.5 1.6 1.6 1.7 1.7 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.7 1.7 1.6 1.6 1.6 1.5 1.5 1.4 1.4 1.4 1.3 1.3 43.9
+21 . 51.52 0.01 19 33 40 0.3 0.4 0.3 0.6 1.0 1.3 0.8 0.9 1.0 1.1 2.6 1.5 1.6 1.6 1.8 1.8 1.8 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.8 1.8 1.8 1.8 1.7 1.7 1.6 1.6 1.5 1.5 1.4 1.4 1.4 1.3 1.3 40.8
+22 . 51.46 0.01 18 31 40 0.2 0.5 0.3 0.7 1.1 1.4 0.9 1.0 1.2 1.3 3.0 1.7 1.7 1.8 1.9 2.0 2.0 2.0 2.1 2.0 2.1 2.1 2.0 2.0 2.0 1.9 1.9 1.9 1.8 1.8 1.7 1.6 1.6 1.5 1.5 1.4 1.4 1.3 1.3 1.2 37.3
+23 . 51.24 0.01 16 28 40 0.3 0.6 0.4 0.8 1.2 1.7 1.0 1.2 1.3 1.5 3.4 1.9 2.0 2.1 2.1 2.2 2.2 2.3 2.3 2.2 2.2 2.2 2.1 2.1 2.1 2.0 1.9 1.8 1.8 1.7 1.7 1.6 1.6 1.5 1.4 1.4 1.3 1.3 1.2 1.2 33.3
+24 . 51.53 0.01 15 26 40 0.3 0.7 0.4 0.9 1.4 1.9 1.2 1.4 1.6 1.7 3.9 2.1 2.2 2.3 2.4 2.4 2.4 2.4 2.4 2.4 2.3 2.3 2.2 2.2 2.1 2.1 1.9 1.9 1.8 1.7 1.6 1.6 1.5 1.5 1.4 1.3 1.2 1.2 1.2 1.1 29.5
+25 . 51.55 0.01 13 24 40 0.2 0.8 0.5 1.0 1.7 2.3 1.4 1.6 1.8 2.0 4.4 2.4 2.5 2.6 2.6 2.6 2.6 2.6 2.5 2.5 2.4 2.3 2.2 2.2 2.1 2.0 1.9 1.8 1.7 1.7 1.6 1.5 1.4 1.4 1.3 1.2 1.2 1.1 1.1 1.0 26.0
+26 . 51.46 0.01 12 22 38 0.2 1.0 0.6 1.2 2.0 2.6 1.6 1.8 2.0 2.2 4.9 2.6 2.7 2.7 2.7 2.7 2.7 2.7 2.6 2.5 2.4 2.4 2.2 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.5 1.4 1.3 1.2 1.2 1.1 1.1 1.0 1.0 23.4
+27 . 51.75 0.01 12 21 36 0.2 1.2 0.7 1.4 2.2 2.9 1.8 2.0 2.2 2.4 5.2 2.8 2.8 2.9 2.8 2.8 2.8 2.7 2.7 2.6 2.4 2.3 2.3 2.1 2.0 2.0 1.8 1.7 1.6 1.6 1.5 1.4 1.3 1.2 1.2 1.1 1.1 1.0 1.0 0.9 21.3
+28 . 51.72 0.01 11 19 33 0.2 1.3 0.8 1.6 2.5 3.3 2.0 2.2 2.4 2.6 5.7 3.0 3.0 3.0 3.0 2.9 2.8 2.7 2.7 2.5 2.4 2.3 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.2 1.1 1.0 1.0 0.9 0.9 0.8 19.1
+29 . 51.45 0.01 10 18 32 0.4 1.6 0.9 1.8 2.8 3.6 2.2 2.4 2.6 2.8 6.0 3.1 3.1 3.1 3.1 2.9 2.9 2.8 2.7 2.6 2.4 2.3 2.2 2.1 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.3 1.2 1.1 1.0 1.0 0.9 0.9 0.8 0.8 17.3
+30 . 51.85 0.01 10 17 29 0.3 1.8 1.1 2.1 3.2 4.1 2.5 2.7 2.8 3.0 6.4 3.3 3.3 3.2 3.1 3.0 2.9 2.8 2.7 2.5 2.4 2.2 2.1 2.0 1.8 1.8 1.6 1.5 1.4 1.4 1.2 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 15.4
+31 . 51.43 0.01 9 16 28 0.2 2.1 1.3 2.4 3.5 4.6 2.6 2.9 3.1 3.2 6.8 3.4 3.4 3.3 3.2 3.1 2.9 2.8 2.6 2.5 2.3 2.2 2.0 1.9 1.8 1.7 1.5 1.4 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 0.6 13.8
+32 . 51.36 0.01 9 15 26 0.3 2.5 1.4 2.6 3.9 4.9 2.9 3.1 3.4 3.4 7.1 3.5 3.4 3.3 3.2 3.1 2.9 2.8 2.6 2.4 2.2 2.1 2.0 1.8 1.7 1.6 1.4 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.8 0.7 0.7 0.6 0.6 12.7
+33 . 51.74 0.01 8 14 24 0.2 2.9 1.7 3.0 4.3 5.5 3.1 3.4 3.6 3.7 7.3 3.6 3.5 3.4 3.2 3.0 2.8 2.7 2.5 2.3 2.2 2.0 1.9 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.9 0.8 0.7 0.7 0.7 0.6 0.6 0.5 11.1
+34 . 51.73 0.02 7 13 23 0.2 3.3 1.9 3.3 4.8 5.9 3.4 3.6 3.8 3.9 7.7 3.7 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.1 1.9 1.7 1.6 1.5 1.4 1.2 1.2 1.1 1.0 0.9 0.9 0.8 0.7 0.7 0.7 0.6 0.6 0.5 0.5 9.7
+35 . 51.63 0.04 7 12 22 0.3 3.9 2.1 3.7 5.2 6.4 3.5 3.8 3.9 4.0 7.8 3.7 3.5 3.4 3.2 3.0 2.7 2.5 2.3 2.1 2.0 1.8 1.7 1.5 1.4 1.3 1.2 1.1 1.0 1.0 0.9 0.8 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 8.8
+36 100.00 51.70 0.06 6 12 22 0.2 4.1 2.4 4.0 5.5 6.7 3.6 3.8 3.9 3.9 7.6 3.6 3.5 3.3 3.0 2.8 2.6 2.4 2.2 2.1 1.9 1.8 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.5 0.4 9.0
+----- ------ ------ ------ -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+ Q1 Q2 Q3 Q= 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
+ -- -- -- ------ ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+Phred.per.base(B) 18 40 40 0.2 0.9 0.5 1.0 1.5 1.9 1.1 1.2 1.3 1.4 2.9 1.5 1.6 1.6 1.6 1.5 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.2 1.1 1.1 1.1 1.0 1.0 1.0 1.0 0.9 0.9 0.9 0.9 0.8 0.8 51.0
+avg.Phred.per.read(R) 26 30 34 0.1 0.0 0.0 0.0 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.6 0.8 1.0 1.3 1.6 2.2 2.7 3.5 4.2 5.1 5.8 6.5 7.0 7.3 7.3 7.1 6.7 6.2 5.5 4.8 3.9 2.9 1.9 0.7
+error.per.read(E) 0 1 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#subsampling.rate: 0.111
```
#### Tab-delimited outputs
-To help with reading, the main statistics for all files can be summarized in tab-delimited format using option `-v t`:
+To help with reading (and to observe faster running times), the main Phred-derived statistics for all files can be summarized in tab-delimited format using option `-v t`:
```bash
-./fastq_info.sh -v t SRR001666*.fastq.gz
+fastq_info.sh -v t SRR001666*.fastq.gz
```
```
-#File NR NB AL BQ1 BQ2 BQ3 RQ1 RQ2 RQ3
-SRR001666_1.fastq.gz 7047668 253716048 36.0 30 40 40 32 35 37
-SRR001666_2.fastq.gz 7047668 253716048 36.0 18 40 40 26 30 34
+#File NR NB AL BQ1 BQ2 BQ3 RQ1 RQ2 RQ3 EQ1 EQ2 EQ3
+SRR001666_1.fastq.gz 7047668 253716048 36.0 30 40 40 32 35 37 0 0 0
+SRR001666_2.fastq.gz 7047668 253716048 36.0 18 40 40 26 30 34 0 1 2
```
-This simple output format enables to easily read every file name (`#File`), no. HTS reads (`NR`) and bases (`NB`), average HTS read length (`AL`), as well as the three quartiles of the global Phred score distribution (`BQ1`, `BQ2`, `BQ3`) and of the average Phred score per HTS read distribution (`RQ1`, `RQ2`, `RQ3`).
+This simple output format enables to easily read every file name (`#File`), no. HTS reads (`NR`) and bases (`NB`), average HTS read length (`AL`), as well as the three quartiles of the three Phred-based distributions, i.e. Phred score per base (`BQ1`, `BQ2`, `BQ3`), average Phred score per HTS read (`RQ1`, `RQ2`, `RQ3`), and expected number of sequencing error(s) per HTS read (`EQ1`, `EQ2`, `EQ3`).
-The above example clearly shows that the overall sequencing error rate is lower in file *SRR001666\_1.fastq.gz* than in file *SRR001666\_2.fastq.gz*.
+The above example clearly shows that the overall sequencing error rate is lower in file *SRR001666\_1.fastq.gz* than in file *SRR001666\_2.fastq.gz*, therefore leading to many more HTS reads without sequencing error in the former FASTQ file.
-#### Note on the subsampling rate
+#### Notes on the subsampling rate
By default, all distributions are estimated from a subset of all input FASTQ blocks to obtain fast running times (the subsampling rate is indicated when using option `-v f`).
-In almost all cases, default subsampling rate (i.e. ~10% with option `-s 5`) is sufficient to efficiently approximate the different distributions (i.e. HTS read lengths, GC-content, Phred scores).
+In almost all cases, default subsampling rate (i.e. ~10% with option `-s 5`) is sufficient to efficiently approximate the different distributions (i.e. HTS read lengths, GC-content, Phred scores, expected numbers of sequencing errors).
For example, the below command line uses all FASTQ blocks from each input file (i.e. option `-s 1`):
```bash
-./fastq_info.sh -s 1 -v t SRR001666*.fastq.gz
+fastq_info.sh -s 1 -v t SRR001666*.fastq.gz
```
```
-#File NR NB AL BQ1 BQ2 BQ3 RQ1 RQ2 RQ3
-SRR001666_1.fastq.gz 7047668 253716048 36.0 30 40 40 32 35 37
-SRR001666_2.fastq.gz 7047668 253716048 36.0 18 40 40 26 30 34
+#File NR NB AL BQ1 BQ2 BQ3 RQ1 RQ2 RQ3 EQ1 EQ2 EQ3
+SRR001666_1.fastq.gz 7047668 253716048 36.0 30 40 40 32 35 37 0 0 0
+SRR001666_2.fastq.gz 7047668 253716048 36.0 18 40 40 26 30 34 0 1 2
```
All statistics are identical to the ones previously estimated (see above), but the overall running time was 8 times slower...
-For comparison, when used with default options, _fastq_info_ is expected to run from 1.5 to 2 times faster than [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc) to process one FASTQ file.
+For comparison, when used with default options, _fastq_info_ is expected to run from 1.5 to 2 times faster than [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc) to process one FASTQ file (and up to twice faster when using option `-v t`).
+
@@ -321,10 +331,14 @@ For comparison, when used with default options, _fastq_info_ is expected to run
Bonfield JK, Mahoney MV (2013) _Compression of FASTQ and SAM format sequencing data_. PLOS One, 8(3):e59190. [doi:10.1371/journal.pone.0059190](https://doi.org/10.1371/journal.pone.0059190).
+Edgar RC, Flyvbjerg H (2015) _Error filtering, pair assembly and error correction for next-generation sequencing reads_. Bioinformatics, 31(21):3476-3482. [doi:10.1093/bioinformatics/btv401](https://doi.org/10.1093/bioinformatics/btv401).
+
+
+Ewing D, Green P (1998) _Base-Calling of Automated Sequencer Traces Using Phred. II. Error Probabilities_. Genome Research, 8:186-194. [doi:10.1101/gr.8.3.186](https://genome.cshlp.org/content/8/3/186).
+
Jones DC, Ruzzo WL, Peng X, Katze MG (2012) _Compression of next-generation sequencing reads aided by highly efficient de novo assembly_. Nucleic Acids Research, 40(22):e171–e171. [doi:10.1093/nar/gks754](https://doi.org/10.1093/nar/gks754).
Roguski L, Deorowicz S (2014) _DSRC 2 - Industry-oriented compression of FASTQ files._ Bioinformatics, 30(15):2213-2215. [doi:10.1093/bioinformatics/btu208](https://doi.org/10.1093/bioinformatics/btu208).
-
diff --git a/fastq_info.sh b/fastq_info.sh
index 23f7d98fec0ef23e6a70f782802fa1376b84fa6f..3c5ca0333a43afbea28ea60f4acfdb367524f1ae 100755
--- a/fastq_info.sh
+++ b/fastq_info.sh
@@ -32,7 +32,12 @@
# = VERSIONS = #
# ============ #
# #
- VERSION=1.0.210325ac #
+ VERSION=2.0.210401ac #
+# + estimates the distribution of the (expected) number of error(s) per HTS read #
+# + modified output tables #
+# + faster running times with option -v t #
+# #
+# VERSION=1.0.210325ac #
# #
##############################################################################################################
@@ -171,7 +176,7 @@ if [ $# -lt 1 ]; then mandoc ; exit 1 ; fi
export LC_ALL=C;
-PHRED=33;
+OFFSET=33;
SPEED=5;
OUT="r"; # r=reduced f=full t=tab-delimited
NTHREADS=1;
@@ -182,7 +187,7 @@ ALTAWK="N.o.N.e";
while getopts p:s:t:v:a:cdh option
do
case $option in
- p) PHRED=$OPTARG ;;
+ p) OFFSET=$OPTARG ;;
s) SPEED=$OPTARG ;;
v) OUT="$OPTARG" ;;
t) NTHREADS=$OPTARG ;;
@@ -195,12 +200,12 @@ do
done
shift "$(( $OPTIND - 1 ))"
-if ! [[ $PHRED =~ ^[0-9]+$ ]]; then echo " incorrect value (option -p): $PHRED" >&2 ; exit 1 ; fi
-if [ $PHRED -lt 33 ]; then echo " Phred should be >32 (option -p): $PHRED" >&2 ; exit 1 ; fi
-if ! [[ $SPEED =~ ^[0-9]+$ ]]; then echo " incorrect value (option -s): $SPEED" >&2 ; exit 1 ; fi
-if ! [[ $NTHREADS =~ ^[0-9]+$ ]]; then echo " incorrect value (option -t): $NTHREADS" >&2 ; exit 1 ; fi
-if ! [[ "rft" =~ $OUT ]] || [ ${#OUT} -ne 1 ]; then echo " incorrect value (option -v): $OUT" >&2 ; exit 1 ; fi
-if [ "$ALTAWK" != "N.o.N.e" ]
+if ! [[ $OFFSET =~ ^[0-9]+$ ]]; then echo " incorrect value (option -p): $OFFSET" >&2 ; exit 1 ; fi
+if [ $OFFSET -lt 33 ]; then echo " Phred offsetshould be >32 (option -p): $OFFSET" >&2 ; exit 1 ; fi
+if ! [[ $SPEED =~ ^[0-9]+$ ]]; then echo " incorrect value (option -s): $SPEED" >&2 ; exit 1 ; fi
+if ! [[ $NTHREADS =~ ^[0-9]+$ ]]; then echo " incorrect value (option -t): $NTHREADS" >&2 ; exit 1 ; fi
+if ! [[ "rft" =~ $OUT ]] || [ ${#OUT} -ne 1 ]; then echo " incorrect value (option -v): $OUT" >&2 ; exit 1 ; fi
+if [ "$ALTAWK" != "N.o.N.e" ]
then
if [ ! $(command -v $ALTAWK) ]; then echo " not found (option -a): $ALTAWK" >&2 ; exit 1 ; fi
AWK="$ALTAWK"; POSIX92=false;
@@ -253,7 +258,7 @@ trap "disown -a ; kill -9 $(jobs -pr) &> /dev/null ; rm -f $TMP1 $TMP2 &>/dev/nu
#### FASTQ INFO ####
#### ####
#############################################################################################################
-[ "$OUT" == "t" ] && echo -e "#File\tNR\tNB\tAL\tBQ1\tBQ2\tBQ3\tRQ1\tRQ2\tRQ3" ;
+[ "$OUT" == "t" ] && echo -e "#File\tNR\tNB\tAL\tBQ1\tBQ2\tBQ3\tRQ1\tRQ2\tRQ3\tEQ1\tEQ2\tEQ3" ;
for INFILE in "$@"
do
@@ -300,173 +305,204 @@ do
fi
$DOS2UNIX && READFILE="$READFILE | tr -d '\r'" ;
+ NEEDGC=1; [ "$OUT" == "t" ] && NEEDGC=0; ## option -v t: no need to estimate GC-content per position
+
if ! $POSIX92
then ############################################################################################################################ GNU AWK v4+
$READFILE |
tee >( sed -n '4~4p' | LC_ALL=C wc -l -c > $TMP1 ) |
sed -n "2~$STEP p" |
- $AWK -v phred=$PHRED '
- BEGIN{
- FS="";c=32;while(++c<=126)c2p[sprintf("%c",c)]=c-phred
- }
- (c){ ## (c!=0)?read:phred
- ++ld[c=NF];while(c){($c~/[ATat]/)&&++at[c]||($c~/[CGcg]/)&&++gc[c];--c}
- next; ## here, c==0
- }
- {
- ++ld[c=NF];sq=0;while(c){sq+=(q=c2p[$c]);++qd[c][q];--c}
- ++qrd[int(0.5+sq/NF)]; c=(++np); ## here, c!=0
- }
- END{
- nr=NR-np;
- for(l in ld)(lmax<l+0)&&lmax=l+0;
- c=0;while(++c<=lmax)for(q in qd[c])(qmax<q+0)&&qmax=q+0;
- # (qmax<40)&&qmax=40;
- while(--c>0){
- q=-1;while(++q<=qmax)if(q in qd[c]){qbd[q]+=qd[c][q];nrm[c]+=qd[c][q]}
- anrm+=nrm[c];
- }
-
- printf"-------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
- printf"pos Lfreq GCfreq Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
- printf"----- ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- c=0;
- while(++c<=lmax){
- printf("%-5s", c); ## pos or lgt
- x=((c in ld)?sprintf("%.2f", 100*ld[c]/(nr+np)):".");printf(" %6s", x); ## lgt freq
- x=((c in gc)?sprintf("%.2f", 100*gc[c]/(at[c]+gc[c])):".");printf(" %6s ", x); ## GC% at pos c
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if(q in qd[c]){f=100*qd[c][q]/nrm[c];cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## Q freq at pos c
- (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
- }
-
- printf"----- ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- printf" Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
- printf" -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- printf"bases ";
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if(q in qbd){f=100*qbd[q]/anrm;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## Q freq
- (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
+ $AWK -v offset=$OFFSET -v ngc=$NEEDGC '
+ BEGIN{
+ FS="";c=32;while(++c<=126){c2q[sprintf("%c",c)]=(q=c-offset);q2p[q]=10^(-q/10)}
+ }
+ (c){ ## (c!=0)?read:phred
+ if(! ngc){++ld[NF];c=0;next}
+ ++ld[c=NF];
+ while(c){($c~/[ATat]/)&&++at[c]||($c~/[CGcg]/)&&++gc[c];--c}
+ next; ## here, c==0
+ }
+ {
+ ++ld[c=NF];sq=e=0;while(c){sq+=(q=c2q[$c]);e+=q2p[q];++qd[c][q];--c}
+ ++qrd[int(0.5+sq/NF)];++erd[int(e)]; c=(++np); ## here, c!=0
+ }
+ END{
+ nr=NR-np;
+ for(l in ld)(lmax<l+0)&&lmax=l+0;
+ c=0;while(++c<=lmax)for(q in qd[c])(qmax<q+0)&&qmax=q+0;
+ # (qmax<40)&&qmax=40;
+ while(--c>0){
+ q=-1;while(++q<=qmax)if(q in qd[c]){qbd[q]+=qd[c][q];nrm[c]+=qd[c][q]}
+ anrm+=nrm[c];
+ }
+
+ printf"-------------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
+ printf"n Lfreq GCfreq Efreq Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
+ printf"----- ------ ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ c=0;
+ printf("%-5s . .", c);
+ x=((c in erd)?sprintf("%.2f", 100*erd[c]/np):".");printf(" %6s", x); ## expected 0 error freq
+ out=" . . . ";q=-1;while(++q<=qmax)out=out" .";print out;
+ while(++c<=lmax){
+ printf("%-5s", c); ## n
+ x=((c in ld)?sprintf("%.2f", 100*ld[c]/(nr+np)):".");printf(" %6s", x); ## lgt freq
+ x=((c in gc)?sprintf("%.2f", 100*gc[c]/(at[c]+gc[c])):".");printf(" %6s", x); ## GC% at pos c
+ x=((c in erd)?sprintf("%.2f", 100*erd[c]/np):".");printf(" %6s ", x); ## expected error freq
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if(q in qd[c]){f=100*qd[c][q]/nrm[c];cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## Q freq at pos c
+ (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+ }
+
+ printf"----- ------ ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ printf" Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
+ printf" -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ printf"Phred.per.base(B) ";
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if(q in qbd){f=100*qbd[q]/anrm;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## Q freq
+ (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+
+ printf"avg.Phred.per.read(R) ";
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if(q in qrd){f=100*qrd[q]/np;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## read freq per avg Q
+ (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
- printf"reads ";
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if(q in qrd){f=100*qrd[q]/np;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## read freq per avg Q
- (q1<0&&cdf>=25)&&q1=q;(q2<0&&cdf>=50)&&q2=q;(q3<0&&cdf>=75)&&q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
- printf"-------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
- print((nr<np)?np:nr);
- }' > $TMP2 ;
+ printf"error.per.read(E) ";
+ q1=q2=q3=c=-1;cdf=0;
+ while(++c<=lmax)if(c in erd){cdf+=100*erd[c]/np;(q1<0&&cdf>=25)&&q1=c;(q2<0&&cdf>=50)&&q2=c;(q3<0&&cdf>=75)&&q3=c}
+ out="";q=-1;while(++q<=qmax)out=out" .";
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+ printf"-------------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
+ print((nr<np)?np:nr);
+ }' > $TMP2 ;
else ############################################################################################################################ AWK POSIX 1003.2
$READFILE |
tee >( sed -n '4~4p' | LC_ALL=C wc -l -c > $TMP1 ) |
sed -n "2~$STEP p" |
- $AWK -v phred=$PHRED '
- BEGIN{
- FS="";c=32;while(++c<=126)c2p[sprintf("%c",c)]=c-phred
- }
- (c){ ## (c!=0)?read:phred
- ++ld[c=NF];while(c){($c~/[ATat]/)&&++at[c]||($c~/[CGcg]/)&&++gc[c];--c}
- next; ## here, c==0
- }
- {
- ++ld[c=NF];sq=0;while(c){sq+=(q=c2p[$c]);++qd[c,q];--c}
- ++qrd[int(0.5+sq/NF)]; c=(++np); ## here, c!=0
- }
- END{
- nr=NR-np;
- for(e in qd){split(e, cq, SUBSEP);if(lmax<cq[1]+0)lmax=cq[1]+0;if(qmax<cq[2]+0)qmax=cq[2]+0;}
- # (qmax<40)&&qmax=40;
- c=lmax+1;
- while(--c>0){
- q=-1;while(++q<=qmax)if((c,q) in qd){qbd[q]+=qd[c,q];nrm[c]+=qd[c,q]}
- anrm+=nrm[c];
- }
+ $AWK -v offset=$OFFSET -v ngc=$NEEDGC '
+ BEGIN{
+ FS="";c=32;while(++c<=126){c2q[sprintf("%c",c)]=(q=c-offset);q2p[q]=10^(-q/10)}
+ }
+ (c){ ## (c!=0)?read:phred
+ if(! ngc){++ld[NF];c=0;next}
+ ++ld[c=NF];
+ while(c){($c~/[ATat]/)&&++at[c]||($c~/[CGcg]/)&&++gc[c];--c}
+ next; ## here, c==0
+ }
+ {
+ ++ld[c=NF];sq=e=0;while(c){sq+=(q=c2q[$c]);e+=q2p[q];++qd[c,q];--c}
+ ++qrd[int(0.5+sq/NF)];++erd[int(e)]; c=(++np); ## here, c!=0
+ }
+ END{
+ nr=NR-np;
+ for(e in qd){split(e, cq, SUBSEP);if(lmax<cq[1]+0)lmax=cq[1]+0;if(qmax<cq[2]+0)qmax=cq[2]+0}
+ # (qmax<40)&&qmax=40;
+ c=lmax+1;
+ while(--c>0){
+ q=-1;while(++q<=qmax)if((c,q) in qd){qbd[q]+=qd[c,q];nrm[c]+=qd[c,q]}
+ anrm+=nrm[c];
+ }
- printf"-------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
- printf"pos Lfreq GCfreq Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
- printf"----- ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- c=0;
- while(++c<=lmax){
- printf("%-5s", c); ## pos or lgt
- x=((c in ld)?sprintf("%.2f", 100*ld[c]/(nr+np)):".");printf(" %6s", x); ## lgt freq
- x=((c in gc)?sprintf("%.2f", 100*gc[c]/(at[c]+gc[c])):".");printf(" %6s ", x); ## GC% at pos c
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if((c,q) in qd){f=100*qd[c,q]/nrm[c];cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## Q freq at pos c
- if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
- }
-
- printf"----- ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- printf" Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
- printf" -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
- printf"bases ";
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if(q in qbd){f=100*qbd[q]/anrm;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## Q freq
- if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
+ printf"-------------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
+ printf"n Lfreq GCfreq Efreq Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
+ printf"----- ------ ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ c=0;
+ printf("%-5s . .", c);
+ x=((c in erd)?sprintf("%.2f", 100*erd[c]/np):".");printf(" %6s", x); ## expected 0 error freq
+ out=" . . . ";q=-1;while(++q<=qmax)out=out" .";print out;
+ while(++c<=lmax){
+ printf("%-5s", c); ## n
+ x=((c in ld)?sprintf("%.2f", 100*ld[c]/(nr+np)):".");printf(" %6s", x); ## lgt freq
+ x=((c in gc)?sprintf("%.2f", 100*gc[c]/(at[c]+gc[c])):".");printf(" %6s", x); ## GC% at pos c
+ x=((c in erd)?sprintf("%.2f", 100*erd[c]/np):".");printf(" %6s ", x); ## expected error freq
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if((c,q) in qd){f=100*qd[c,q]/nrm[c];cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## Q freq at pos c
+ if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+ }
+
+ printf"----- ------ ------ ------ -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ printf" Q1 Q2 Q3 Q= 0";q=0;while(++q<=qmax)printf(" %5s", q);print"";
+ printf" -- -- -- ------";q=0;while(++q<=qmax)printf" -----";print"";
+ printf"Phred.per.base(B) ";
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if(q in qbd){f=100*qbd[q]/anrm;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## Q freq
+ if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+
+ printf"avg.Phred.per.read(R) ";
+ cdf=0;out="";q1=q2=q3=q=-1;
+ while(++q<=qmax){
+ if(q in qrd){f=100*qrd[q]/np;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
+ out=out""sprintf(" %5s", x); ## read freq per avg Q
+ if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
+ }
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
- printf"reads ";
- cdf=0;out="";q1=q2=q3=q=-1;
- while(++q<=qmax){
- if(q in qrd){f=100*qrd[q]/np;cdf+=f;x=sprintf("%.1f", f)}else{x="."}
- out=out""sprintf(" %5s", x); ## read freq per avg Q
- if(q1<0&&cdf>=25)q1=q;if(q2<0&&cdf>=50)q2=q;if(q3<0&&cdf>=75)q3=q;
- }
- q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
- print q1" "q2" "q3" "out;
- printf"-------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
- print((nr<np)?np:nr);
- }' > $TMP2 ;
+ printf"error.per.read(E) ";
+ q1=q2=q3=c=-1;cdf=0;
+ while(++c<=lmax)
+ if(c in erd){cdf+=100*erd[c]/np;if(q1<0&&cdf>=25)q1=c;if(q2<0&&cdf>=50)q2=c;if(q3<0&&cdf>=75)q3=c}
+ out="";q=-1;while(++q<=qmax)out=out" .";
+ q1=sprintf("%2s", q1);q2=sprintf("%2s", q2);q3=sprintf("%2s", q3);
+ print q1" "q2" "q3" "out;
+ printf"-------------------------------------------";q=0;while(++q<=qmax)printf"------";print"";
+ print((nr<np)?np:nr);
+ }' > $TMP2 ;
fi
+
wait ;
- NR=`$AWK '{print$1}' $TMP1`;
- NB=`$AWK '{print$2}' $TMP1`;
+ NR=$($AWK '{print$1}' $TMP1);
+ NB=$($AWK '{print$2}' $TMP1);
NB=$(( $NB - $NR ));
- AL=`bc -l <<<"scale=1;$NB/$NR"`;
- SS=`sed -n '$p' $TMP2`;
+ AL=$(bc -l <<<"scale=1;$NB/$NR");
+ SS=$(sed -n '$p' $TMP2);
if [ "$OUT" == "r" ] || [ "$OUT" == "f" ]
then
echo "##File: $INFILE" ;
- echo "#no.reads: $NR" ;
- echo "#no.bases: $NB" ;
- echo "#avg.lgt: $AL" ;
+ echo "#no.reads(NR): $NR" ;
+ echo "#no.bases(NB): $NB" ;
+ echo "#avg.lgt(AL): $AL" ;
if [ "$OUT" == "r" ]
then
- sed '$d' $TMP2 | cut -c1-30 ;
+ sed '$d' $TMP2 | cut -c1-36 ;
else
sed '$d' $TMP2 ;
- echo "#subsampling.rate: "`bc -l <<<"scale=3;$SS/$NR" | sed 's/^0$/0.00/;s/^\./0./'` ;
+ echo "#subsampling.rate: "$(bc -l <<<"scale=3;$SS/$NR" | sed 's/^0$/0.00/;s/^\./0./') ;
fi
echo ;
elif [ "$OUT" == "t" ]
then
- BQ123="`tail $TMP2 | grep "^bases " $TMP2 | cut -c22-30 | xargs echo | tr ' ' '\t'`";
- RQ123="`tail $TMP2 | grep "^reads " $TMP2 | cut -c22-30 | xargs echo | tr ' ' '\t'`";
- echo -e "$INFILE\t$NR\t$NB\t$AL\t$BQ123\t$RQ123" ;
+ BQ123="$(tail $TMP2 | grep -F "(B)" $TMP2 | cut -c28-36 | tr -s ' ' | sed 's/^ //' | tr ' ' '\t')";
+ RQ123="$(tail $TMP2 | grep -F "(R)" $TMP2 | cut -c28-36 | tr -s ' ' | sed 's/^ //' | tr ' ' '\t')";
+ EQ123="$(tail $TMP2 | grep -F "(E)" $TMP2 | cut -c28-36 | tr -s ' ' | sed 's/^ //' | tr ' ' '\t')";
+ echo -e "$INFILE\t$NR\t$NB\t$AL\t$BQ123\t$RQ123\t$EQ123" ;
fi
echo -n > $TMP1 ;