2
$\begingroup$

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
Improve this question
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – terdon ♦
    Feb 22, 2018 at 11:03
  • 3
    $\begingroup$ should CGCGC as one occurrence or two? $\endgroup$
    – heathobrien
    Feb 22, 2018 at 11:21
  • $\begingroup$ number of times the pattern is seen (exclusive numbers, no overlap) $\endgroup$
    – 719016
    Feb 22, 2018 at 12:12
  • $\begingroup$ 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. $\endgroup$
    – user172818
    Feb 22, 2018 at 15:01

1 Answer 1

Reset to default
1
$\begingroup$

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
Improve this answer
$\endgroup$
4
  • $\begingroup$ Brilliant, sometimes awk is all you need! Would pigz speed-up the process? $\endgroup$
    – 719016
    Feb 22, 2018 at 13:09
  • 1
    $\begingroup$ @719016 good idea! And yes, it certainly does speed it up, see updated answer. $\endgroup$
    – terdon ♦
    Feb 22, 2018 at 14:28
  • 1
    $\begingroup$ Check out bioawk github.com/lh3/bioawk . This prevents counting CGC located in other lines. I wonder if it is faster than your current solutions? $\endgroup$
    – conchoecia
    Feb 22, 2018 at 15:45
  • 2
    $\begingroup$ @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. $\endgroup$
    – terdon ♦
    Feb 22, 2018 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.