来自vcountPattern的正确命中的提取序列R

我进行了小RNA测序并尝试分析结果fastq文件。来自vcountPattern的正确命中的提取序列R

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

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

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

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

我可以在这里做的是数数的读取包含 “TCTGCATTTAAGGCAAGTT”

sum (hit)

[1] 11500

因此，有11500序列读取包含“TCTGCATTTAAGGCAAGTT”

但在此之上，w我想要的帽子是从fastq文件中提取对应于11500次读取的实际序列。

我该如何做到这一点？

hit

如果我只是这样做，它会给出一堆'0'，少量'1'，很少'2'。所以我相信这基本上是一个对应于每次阅读中点击次数的向量。

我试图使用这些信息提取序列信息，但无法实现。

任何帮助被赞赏！

来源

2015-08-25 gdy

供参考：用户正在使用Bioconducter软件包“ShortRead”https://darrenjw.wordpress.com/2010/11/29/a-quick-introduction-to-the-bioconductor-shortread-package-for-the-分析-的-NGS-数据/。除非你能给我们一个玩具fq文件，否则不容易复制这些代码。序列分析知识在这里很有用。 – Sean

亲爱的霍姆斯，我准备了一个玩具fastq文件，您可以从这里下载[link]（https://drive.google.com/file/d/0ByEbUQPY_T_oci1fbDFHSHQ4WUk/view?usp=sharing）。当我使用这个fastq文件尝试我的脚本时，有3个正面点击。基本上我只想从fastq文件中提取正面的点击。我的原始fastq文件的大小比这大200倍。 – gdy

不要介意福尔摩斯，我看着你提供的链接，我已经从中得到了答案。（读[hit]）解决了问题 – gdy

顾名思义，vcountPattern只有计数模式匹配。它不会提供你的位置。使用vmatchPattern。不幸的是，这个功能不支持with.indels = TRUE（还？） - 这既烦人又难以理解。

但是，您可以改用matchPattern。由于matchPattern只能对一个序列进行操作，而不是一组，你需要的功能，手动应用到XStringSet：

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

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

来源

2015-08-25 19:53:03

谢谢康拉德，我认为你的方法应该可以工作，但它没有。基本上它不会选择正面命中，而是显示所有fastq序列，例如“关于31个字母的DNAString主题的视图主题：GCATTGGTGGATCAGTGGTAGAATTCTCGCC views：NONE”。这种条目的数量与原始fastq的序列数量相同。 – gdy

我准备了玩具fastq文件，这是我的fastq文件的一小部分。如果你尝试这个，我有3个积极的命中。基本上我想从这个fastq文件中提取三个正面点击的序列信息，并丢弃大多数其他负面点击。你可以在下面的链接下载玩具fastq文件[链接]（https://drive.google.com/file/d/0ByEbUQPY_T_oci1fbDFHSHQ4WUk/view?usp=sharing） – gdy

亲爱的康拉德，我看了一下福尔摩斯提供的链接和我已经从中得到了答案。 sread（读取[hit]）将仅提取正序列读数。但我感谢你的帮助！ – gdy

来自vcountPattern的正确命中的提取序列R

回答

相关问题