r - 从R中的vcountPattern中提取正命中序列

Question

我进行了小RNA测序并尝试分析结果fastq文件。

首先，我使用 ShortRead 包将 fastq 文件导入 R 并转换为 DNAstringSet

reads <- readFastq("test.fq")
seq <- sread(reads)

为了查找包含特定序列字符串的读取，我使用了来自 Biostrings 库的 vcountPattern。为了我的分析目的，我必须允许突变和插入缺失。

hit <-vcountPattern("TCTGCATTTAAGGCAAGTT", seq, max.mismatch=5, with.indels=TRUE)

我可以从这里做的是计算包含“TCTGCATTTAAGGCAAGTT”的读取次数

sum (hit)

它返回

[1] 11500

所以有 11500 个包含“TCTGCATTTAAGGCAAGTT”的序列读取

但最重要的是，我想要的是从 fastq 文件中提取对应于 11500 读取的实际序列。

我怎样才能做到这一点？

hit

如果我只是这样做，它会给出一堆“0”，少量的“1”，很少的“2”。所以我相信这基本上是一个与每次读取的命中数相对应的向量。

我尝试使用此信息提取序列信息，但无法实现。

任何帮助表示赞赏！

score 1 · Accepted Answer

顾名思义，vcountPattern只计算模式匹配。它不为您提供位置。为此使用vmatchPattern。不幸的是，这个函数不支持with.indels = TRUE（还没有？）——这既烦人又有点难以理解。¹

但是，您可以matchPattern改用。由于matchPattern仅对单个序列而不是集合进行操作，因此您需要手动将该函数应用于XStringSet：

hits = lapply(seq, matchPattern,
              pattern = "TCTGCATTTAAGGCAAGTT",
              max.mismatch = 5, with.indels = TRUE)

¹表面上的原因可能vmatchPattern是使用与不同的算法实现matchPattern，并且该算法不支持 indels。然而，没有一个很好的理由不简单地为lapply我们上面使用的提供一个包装器。

1 回答 1