从R中的vcountPattern中提取正命中序列答案

【问题标题】：Extract sequence of positive hit from vcountPattern in R从R中的vcountPattern中提取正命中序列
【发布时间】：2015-11-19 15:47:36
【问题描述】：

我进行了小 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”。所以我相信这基本上是一个与每次读取的命中数相对应的向量。

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

感谢任何帮助！

【问题讨论】：

仅供参考：用户正在使用 Bioconducter 包“ShortRead”darrenjw.wordpress.com/2010/11/29/…。除非您可以给我们一个玩具 fq 文件，否则不容易复制此代码。序列分析知识在这里很有用。
亲爱的福尔摩斯，我准备了一个玩具 fastq 文件，你可以从这里下载 link。当我使用这个 fastq 文件尝试我的脚本时，有 3 个正面命中。基本上我只想从 fastq 文件中提取正面命中。我原来的 fastq 文件的大小比这个大 200 倍。
别管福尔摩斯，我查看了您提供的链接，并从中得到了答案。 sread(reads[hit]) 解决了问题

【解决方案1】：

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

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

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

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

【讨论】：

谢谢康拉德，我认为你的方法应该有效，但它没有。基本上它不会选择正面命中，而是显示所有 fastq 序列，例如“查看 31 个字母的 DNAString 主题主题：GCATTGGTGGATCAGTGGTAGAATTCTCGCC 视图：无”。此类条目的数量与原始 fastq 中的序列数量相同。
我准备了玩具fastq文件，这是我的fastq文件的一小部分。如果你试试这个，我有 3 个正面的命中。基本上我想从这个fastq文件中提取三个正命中的序列信息，并丢弃大部分其他负命中。您可以在以下链接中下载玩具 fastq 文件link
亲爱的康拉德，我查看了福尔摩斯提供的链接，并从中得到了答案。 sread(reads[hit]) 将仅提取正序列读取。但我感谢您的帮助！