我正在 Linux 集群上使用 bash。如果它们包含与查询序列的匹配项,我正在尝试从 .fastq 文件中提取读数。下面是一个包含三个读数的示例 .fastq 文件。
$ cat example.fastq
@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.2 2/1
CTATANTATTCTATATTTATTCTAGATAAAAGCATTCTATATTTAGCATATGTCTAGCAAAAAAAA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
我想提取包含序列 GAAATAATA 的读数。我可以使用 grep 执行此提取,如以下命令所示。
$ grep -F -B 1 -A 2 "GAAATAATA" example.fastq > MATCH.fastq
$猫匹配.fastq
@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
但是,此策略不容忍任何不匹配。例如,包含序列 GAAATGATA 的读取将被忽略。我需要这个提取来容忍查询序列中任何位置的一个不匹配。所以我的问题是我怎样才能做到这一点?是否有可用的序列比对包,其功能与 grep 类似?是否有可用的执行此类操作的 fastq 子集包?需要注意的是速度非常重要。感谢您的指导。
最佳答案
这是一个解决方案,使用 agrep
获取匹配项的记录数,并使用 awk 打印出具有某些上下文的记录(由于缺少 -A
和 -B
in agrep
):
$ agrep -1 -n "GAAATGATA" file |
awk -F: 'NR==FNR{for(i=($1-1);i<=($1+2);i++)a[i];next}FNR in a' - file
输出:
@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
关于awk - 容忍子集 .fastq 不匹配的 Grep,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/53768447/