生信（十）利用kseq.h和regex.h实现类似grep查找f

原创：hxj7

本文给出了一个利用kseq.h和regex.h实现类似grep查找fastq reads功能的示例（C语言）。

引出问题

做生信的朋友应该都很熟悉类Unix系统中的grep命令，该命令可以快速查找并输出包含目标字符串的行。在对fastq文件进行处理时，我们有时候需要查找包含特定字符串的reads。因为一个reads包含了多行，所以grep命令不能完全适用。那有没有其它命令或者工具可以实现快速简便地实现上述查找特定reads的功能呢？就像grep快速查找行一样。

在《生信（八）zlib库操作fq-gz文件》一文中，我们分享过一个例子：
<center>如何输出第一行（name行）结尾是ACCGAATG的所有reads？</center>

image
当时我们提到可以利用zcat配合sed或者awk命令比较精确地查找特定reads：

zcat test.fq.gz | sed –n ‘h;N;N;N;x;/:ACCGAATG$/{x;p}’ | gzip –c > out3.fq.gz

或者

zcat test.fq.gz | awk –f index.awk | gzip –c > out4.fq.gz

# 具体的awk命令如下：
 #!/usr/bin/awk
{
        f[0] = $0;
        for (i = 1; i < 4; i++) {
                getline;
                f[i] = $0;
        }
        if (f[0] ~ /:ACCGAATG$/) {
                for (i = 0; i < 4; i++)
                        print f[i];
        }
}

在最近的回顾中，笔者发现上面的两种解决方式只适用于reads只占4行的情况，如果一个reads超过4行就会出错：比如下面这样一个reads就不会被输出，并且可能会导致上述sed和awk命令的运行结果出错：

@K00137:236:H7NLVBBXX:6:1126:29721:23241 1:N:0:ACCGAATG
TGGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAA
ATACTGGAGACAGTCAGCTTATTTATCAGAAAGGTTTATT
+
```eeiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
iiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiieiiiii

解决问题

为此，在解析fastq文件时需要考虑到reads占据超过4行的情况。lh3编写的kseq.h已经可以很好地处理这个问题。而类似grep那样强大的查找功能可以通过regex.h这个头文件来实现，regex.h是C语言中支持正则表达的一个库。

笔者利用kseq.h和regex.h编写了一段代码，可以解决上述问题：
<center>如何输出第一行（name行）结尾是ACCGAATG的所有reads？</center>

代码运行效果如下：

更多的测试：

具体代码如下：

#include <zlib.h> 
#include <stdio.h>
#include <regex.h>
#include <limits.h>
#include "kseq.h"  

// declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)  

/* the stk_printstr function is from https://github.com/lh3/seqtk */
static void stk_printstr(FILE *fp, const kstring_t *s, unsigned line_len) {
    if (line_len != UINT_MAX && line_len != 0) {
        int i, rest = s->l;
        for (i = 0; i < s->l; i += line_len, rest -= line_len) {
            fputc('\n', fp);
            if (rest > line_len) fwrite(s->s + i, 1, line_len, fp);
            else fwrite(s->s + i, 1, rest, fp);
        }
        fputc('\n', fp);
    } else {
        fputc('\n', fp);
        fputs(s->s, fp);
        fputc('\n', fp);
    }
}

/* the stk_printseq2 function is modified from https://github.com/lh3/seqtk */
static inline void stk_printseq2(FILE *fp, const kseq_t *s, int line_len) {
    fputc('@', fp);
    fputs(s->name.s, fp);
    if (s->comment.l) {
        fputc(' ', fp); 
        fputs(s->comment.s, fp);
    }
    stk_printstr(fp, &s->seq, line_len);
    if (s->qual.l) {
        fputc('+', fp);
        stk_printstr(fp, &s->qual, line_len);
    }
}
  
int main(int argc, char *argv[]) {  
    gzFile fp;
    kseq_t *seq;
    regex_t regex;  
    int l, reti, regerr;
    char msgbuf[100];   
    
    if (argc <= 2) {  
        fprintf(stderr, "Usage: %s <in.seq> <regexp>\n", argv[0]);
        return 1;  
    }  
    reti = regcomp(&regex, argv[2], 0);   /* Compile regular expression */
    if (reti) {
        fprintf(stderr, "Could not compile regex\n"); 
        return 3; 
    }
    fp = gzopen(argv[1], "r");   // open the file handler 
    seq = kseq_init(fp);     // initialize seq  
    for (regerr = 0; (l = kseq_read(seq)) >= 0; ) {   // read sequence 
        reti = regexec(&regex, seq->comment.s, 0, NULL, 0);  /* Execute regular expression */
        if (! reti) {  /* RegExp Match */
            stk_printseq2(stdout, seq, UINT_MAX);
        } else if (reti != REG_NOMATCH) {  /* RegExp Error */
            regerror(reti, &regex, msgbuf, sizeof(msgbuf));
            fprintf(stderr, "Regex match failed: %s\n", msgbuf);            
            regerr = 1;
            break;
        }
    }
    kseq_destroy(seq);   // destroy seq
    gzclose(fp);      // close the file handler
    regfree(&regex);  /* Free compiled regular expression */
    return regerr ? 5 : 0;
}  

对代码的说明：

1. kseq.h中的seq->name.s（即reads的sample name）是不包含开头的'@'符号的，所以在输出name行时要首先输出'@'符号；
2. reads第一行中空格后面的部分是保存在seq->comment.s中的；
3. 上面的代码只能针对seq->comment.s中的字符串进行匹配，如果需要对reads的其它部分进行匹配，可以对上述代码做适当改编。比如，要针对sample name进行匹配，可以将代码中的seq->comment.s改为seq->name.s即可。

最后给出用作测试的fastq文件，一共6个reads。
（注意：如果要做测试，请将下面的reads保存成以'\n'为行结尾的文件格式。）

@K00137:236:H7NLVBBXX:6:1126:29721:23241 1:N:0:ACCGAATG
TGGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAA
ATACTGGAGACAGTCAGCTTATTTATCAGAAAGGTTTATT
+
```eeiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
iiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiieiiiii
@K00137:236:H7NLVBBXX:6:1117:32410:45906 1:N:0:ACCGAATG
NCATTCATTATCTCAGCACCGGCATCACGCACGCGGTCTACATAACGGCCCGGCTCGGCC
ACCATCATGTGGACATCCAGAGGTTTTTCGGCAATGGTGC
+
B``eeiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiiiiiiiiiiieii`i`eii
[eiiiieVeeieiii[`ei``e[L``eiiiii`i`i`[ei
@K00137:236:H7NLVBBXX:6:1210:12642:24982 1:N:0:GTAAGCCA
GGTCTTTTGGTTATGTACGGAATTCAAAAAATTGATCAGTGTTACGCAGGTCTTGATTTG
GGATCTTTTCTAACCATTGCGGTGATTGCAGCAGGCTGCG
+
``eeeiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
@K00137:236:H7NLVBBXX:6:1105:13971:18669 1:N:0:GTAAGCCA
TATTTCAGCGTGATAACTACCTTCAAGGGCTTTCAAGTTATCCGTGTGTCGCAAGTGCGA
TATGTCAGCCGTAAGGGAGAGCCTATGCAGTTCTACTGTC
+
``eeeiiiiiiiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiieeiiiiiiiiii
iiiiiiiiiiieiieieV[eieeeiiiiiiiiiiiiiiii
@K00137:236:H7NLVBBXX:6:1105:21816:9789 1:N:0:GTAAGCCA
ACCTTCCGTCTACATCGCATCTGAACGAGAATCGTCCGCAGCTGTCGGTCATAAAGCTAT
CCATAAAAGACCTCCTCAACGAGTCCGACGGCTCGTTTTC
+
``eeeiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieiiii
@K00137:236:H7NLVBBXX:6:1110:22800:24261 1:N:0:ACCGAATG
TGCTGACCACATGTGCGCTTATCCTTATATTCACTAAAACCAATGGAATTACGCCCACTC
AAGGTAGTGTATTCATAGCCGGAATGCAGGCTGTGATTGC
+
``eeeiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii

（公众号：生信了）