Answers.Body: <p>Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:</p>
  
  <pre><code>zcat file.fastq.gz | grep CGC | 
      awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
           sort -rnk2   
  </code></pre>
  
  <p>I tested this on a 2.6G fastq file (exome data) and it took:</p>
  
  <pre><code>real    1m18.987s
  user    1m41.392s
  sys     0m10.684s
  </code></pre>
  
  <p>On a 66G WGS sequence file, it took:</p>
  
  <pre><code>real    34m6.771s
  user    40m16.452s
  sys     4m10.176s
  </code></pre>
  
  <p>The extra <code>grep CGC</code> isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the <code>grep</code>:</p>
  
  <pre><code>zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
      sort -rnk2   
  </code></pre>
  
  <hr>
  
  <p>Following the OP's suggestion, I also tested with <code>pigz</code> and, it turns out, pigz can fly!</p>
  
  <pre><code>time unpigz -c exome.fastq.gz | grep CGC | 
      awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
           sort -rnk2
  
  
  real    0m44.012s
  user    1m19.544s
  sys     0m14.900s
  </code></pre>
  
  <p>And for the 66G WGS:</p>
  
  <pre><code>real    17m29.935s
  user    29m23.816s
  sys     5m17.856s
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-22T13:09:30.703Z
Answers.Text: Brilliant, sometimes awk is all you need! Would pigz speed-up the process?
Answers.UserDisplayName: 719016
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-22T14:28:22.810Z
Answers.Text: @719016 good idea! And yes, it certainly does speed it up, see updated answer.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-22T15:45:42.347Z
Answers.Text: Check out bioawk https://github.com/lh3/bioawk . This prevents counting CGC located in other lines. I wonder if it is faster than your current solutions?
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-22T16:10:46.907Z
Answers.Text: @conchoecia I doubt it because I am taking advantage of `awk`'s internal splitting by setting the field separator (`FS`) to `CGC` and using the `%` function to ignore irrelevant lines. Using `bioawk` with `-c fastx` adds a whole level of needless complexity since that will also look at the name, quality and comment lines as well. It also means I can't just use `FS` but would need to do something much splower like `split()`. And, indeed, I just tested on the exome file and stopped it after 7 minutes since it was taking too long already.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-22T12:36:04.043Z
Answers.LastActivityDate: 2018-02-22T14:28:07.530Z
Answers.LastEditDate: 2018-02-22T14:28:07.530Z
Answers.UserDisplayName: terdon
Body: <p>What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?</p>
  
  <p>I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read? </p>
  
  <p>E.g.</p>
  
  <pre><code>countingtool --pattern CGC file_R1.fastq.gz
  </code></pre>
  
  <p>Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)</p>
  
  <pre><code>CGC    0     51341234
  CGC    1     13948130
  CGC    2      1234344
  CGC    3        51344
  CGC    4         1343
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-22T11:03:00.920Z
Comments.Text: Are those numbers inclusive? I mean, are the 1343 reads that have 4 `CGC` included in the `13948130` that have 1? Are the reads with >N included in the count of reads with N?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-22T11:21:30.460Z
Comments.Text: should ```CGCGC``` as one occurrence or two?
Comments.UserDisplayName: heathobrien
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-22T12:12:14.307Z
Comments.Text: number of times the pattern is seen (exclusive numbers, no overlap)
Comments.UserDisplayName: 719016
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-22T15:01:50.467Z
Comments.Text: Are reads sequenced from a "stranded" library? If not, you need to consider the reverse complement. I would choose the strand, for each read, that gives the most count.
Comments.UserDisplayName: user172818
Tags: <fastq><benchmarking><seqtk><cutadapt>
Title: efficient counting of dinucleotides/trinucleotides on fastq reads?
ThreadText: # efficient counting of dinucleotides/trinucleotides on fastq reads?
  
  (Asked by: username_0 on 2018-02-22)
  
  What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?
  
  I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read?
  
  E.g.
  
  ```
  countingtool --pattern CGC file_R1.fastq.gz
  
  ```
  
  Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)
  
  ```
  CGC    0     51341234
  CGC    1     13948130
  CGC    2      1234344
  CGC    3        51344
  CGC    4         1343
  
  ```
  
  **Tags:** fastq, benchmarking, seqtk, cutadapt
  
  > **Comments:**
  > **username_1**: Are those numbers inclusive? I mean, are the 1343 reads that have 4 `CGC` included in the `13948130` that have 1? Are the reads with >N included in the count of reads with N?  
  > **username_2**: should ```CGCGC``` as one occurrence or two?  
  > **username_0**: number of times the pattern is seen (exclusive numbers, no overlap)  
  > **username_3**: Are reads sequenced from a "stranded" library? If not, you need to consider the reverse complement. I would choose the strand, for each read, that gives the most count.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-02-22. Score: 1 | Accepted answer
  
  Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:
  
  ```
  zcat file.fastq.gz | grep CGC | 
      awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
           sort -rnk2   
  
  ```
  
  I tested this on a 2.6G fastq file (exome data) and it took:
  
  ```
  real    1m18.987s
  user    1m41.392s
  sys     0m10.684s
  
  ```
  
  On a 66G WGS sequence file, it took:
  
  ```
  real    34m6.771s
  user    40m16.452s
  sys     4m10.176s
  
  ```
  
  The extra `grep CGC` isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the `grep`:
  
  ```
  zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
      sort -rnk2   
  
  ```
  
  ---
  
  Following the OP's suggestion, I also tested with `pigz` and, it turns out, pigz can fly!
  
  ```
  time unpigz -c exome.fastq.gz | grep CGC | 
      awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
           sort -rnk2
  
  real    0m44.012s
  user    1m19.544s
  sys     0m14.900s
  
  ```
  
  And for the 66G WGS:
  
  ```
  real    17m29.935s
  user    29m23.816s
  sys     5m17.856s
  
  ```
  
  > **Comments:**
  > **username_0**: Brilliant, sometimes awk is all you need! Would pigz speed-up the process?  
  > **username_1**: @username_0 good idea! And yes, it certainly does speed it up, see updated answer.  
  > **username_4**: Check out bioawk https://github.com/lh3/bioawk . This prevents counting CGC located in other lines. I wonder if it is faster than your current solutions?  
  > **username_1**: @username_4 I doubt it because I am taking advantage of `awk`'s internal splitting by setting the field separator (`FS`) to `CGC` and using the `%` function to ignore irrelevant lines. Using `bioawk` with `-c fastx` adds a whole level of needless complexity since that will also look at the name, quality and comment lines as well. It also means I can't just use `FS` but would need to do something much splower like `split()`. And, indeed, I just tested on the exome file and stopped it after 7 minutes since it was taking too long already.  
  
  ---
  
  