本文给出了一个利用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
或者
代码语言：javascript复制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？
代码运行效果如下：
更多的测试：
具体代码如下：
代码语言：javascript复制#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(®ex, argv[2], 0); /* Compile regular expression */
if (reti) {
fprintf(stderr, "Could not compile regexn");
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(®ex, 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, ®ex, msgbuf, sizeof(msgbuf));
fprintf(stderr, "Regex match failed: %sn", msgbuf);
regerr = 1;
break;
}
}
kseq_destroy(seq); // destroy seq
gzclose(fp); // close the file handler
regfree(®ex); /* Free compiled regular expression */
return regerr ? 5 : 0;
}
对代码的说明：
- kseq.h中的
seq->name.s（即reads的sample name）是不包含开头的'@'符号的，所以在输出name行时要首先输出'@'符号； - reads第一行中空格后面的部分是保存在
seq->comment.s中的； - 上面的代码只能针对
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
