生信（十）利用kseq.h和regex.h实现类似grep查找fastq reads功能的示例（C语言）

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

引出问题

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

在《生信（八）zlib库操作fq-gz文件》一文中，我们分享过一个例子：

如何输出第一行（name行）结尾是ACCGAATG的所有reads？

当时我们提到可以利用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编写了一段代码，可以解决上述问题：

如何输出第一行（name行）结尾是ACCGAATG的所有reads？

代码运行效果如下：

更多的测试：


具体代码如下：

#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