☰

从 fastq 文件中检索特定序列的读取名称

tag-icon • tag-icon shell-script
从 fastq 文件中检索特定序列的读取名称

我有一个 fastq 文件。读取长度为 75bp。fastq 文件如下所示

@ERR030899.2391 HWI-BRUNOP16X_0001:4:1:16426:3738#0/1
NNCGTGTCAGTGGCCAGCAGCCCACACTGCGCATGGCTGATACTGTGGACCCCCTGGACTGGCTTTTTGGGGAGT
+
###########################################################################
@ERR030899.2392 HWI-BRUNOP16X_0001:4:1:16582:3744#0/1
NNAAAGTCCTGCGCTGCGGAGGACAGGAAGCACCCCCTCAATAGCCAGCACCCACAGTGAGCTGAGCACTTACAG
+
##'(''((((5/333**+)(10-11+1,,,,1/1/F<<<<FF:FFFFFF=FFFFFFFFFFFFFFFFFFFFFFFF#
@ERR030899.2393 HWI-BRUNOP16X_0001:4:1:16911:3745#0/1
NNGGATTTAGCGGGGTGATGCCTGTTGGGGGCCCGTGCCCTCCTACTTGGGGGGCAGGGGCTAGGCTGCAGAGGT
+
###########################################################################
@ERR030899.2394 HWI-BRUNOP16X_0001:4:1:18194:3739#0/1
NNACAAGCAATTTAGTGATAATGTCCAGAGGGCCAAGGATGCGGACCACCTTTTGCAGAACTCATATCTCGAGCA
+
##*+*)'+++FFFFFFFFFF58588==?:?FFFFFFFFFFFFFF<FFFFFFFFFF=FFFFFFFFFFFFFF6=??;

我有一个大约30bp的小核苷酸序列。

假设这是我的核苷酸序列

韓國語: 韓國語

我想要做的是在 fastq 文件中搜索此序列并提取序列所在的相应读取名称。因此，在这种情况下读取名称将是

@ERR030899.2393 HWI-BRUNOP16X_0001:4:1:16911:3745#0/1

我还想允许序列中的不匹配率为 1。

任何脚本或命令行？

谢谢

问候

答案1

awk -v seq="CTGTTGGGGGCCCGTGC" '
  NR%4 == 1 {name = $0}
  NR%4 == 2 && index($0, seq) {print name}
' filename

如果根据“不匹配率为 1”，您希望能够匹配这 30 个中的任意 29 个（例如CTG.TGGGGGCCCGTGC，那就相当复杂了。

...

呃，没那么复杂：

awk -v seq="CTGTTGGGGGCCCGTGC" '
  NR%4 == 1 {name=$0}
  NR%4 == 2 {
    if (index($0, seq))
      print "found seq \"" seq "\" in " name
    else
      for (i=1; i<=length(seq); i++) {
        patt = substr(seq, 1, i-1) "." substr(seq, i+1)
        if (match($0, patt)) {
          print "found pattern \"" patt "\" in " name
          break
        }
      }
  }
' filename

相关内容