允许与子集.fastq不匹配的Grep [英] Grep that tolerates mismatches to subset .fastq

问题描述

我正在Linux集群上使用bash.我正在尝试从.fastq文件中提取读取，如果它们包含与查询序列匹配的内容.下面是一个包含三个读取的.fastq文件示例.

I am working with bash on a linux cluster. I am trying to extract reads from a .fastq file if they contain a match to a queried sequence. Below is an example .fastq file containing three reads.

$ 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执行此提取，如以下命令所示.

I would like to extract reads containing the sequence GAAATAATA. I can perform this extraction using grep as shown in the following command.

$ grep -F -B 1 -A 2 "GAAATAATA" example.fastq > MATCH.fastq

$ cat MATCH.fastq

@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6

但是，此策略不能容忍任何不匹配.例如，包含序列GAAAT G ATA的读物将被忽略.我需要这种提取来容忍查询序列中任何位置的一个不匹配.所以我的问题是如何实现这一目标?是否存在与grep具有类似功能的序列比对软件包?有没有可用的fastq子设置包可以执行这种类型的操作?需要注意的是速度非常重要.感谢您的指导.

However, this strategy does not tolerate any mismatches. For example, a read containing the sequence GAAATGATA will be ignored. I need this extraction to tolerate one mismatch at any position in the queried sequence. So my question is how can I achieve this? Is there a sequence alignment package available with similar functionality to grep? Are there any fastq subsetting packages available that perform this type of operation? One note is that speed is very important. Thanks for your guidance.

推荐答案

以下是使用agrep来获取匹配记录数的解决方案，以及awk可以在某些上下文中打印出这些记录(由于缺少-A和agrep中的-B):

Here is a solution using agrep to get the record numbers of matches and an awk that prints out those records with some context (due to missing -Aand -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

这篇关于允许与子集.fastq不匹配的Grep的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！