在 DNA 序列中搜索字符串答案

【问题标题】：Search for characters string in a DNA sequence在 DNA 序列中搜索字符串
【发布时间】：2017-04-15 14:42:16
【问题描述】：

我正在尝试查看基因序列中的某些核苷酸模式。我刚刚完成了read.table 来获取它，但我也尝试将它转换为向量和数据帧。

我如何搜索一个模式（例如AACG）甚至只是一个核苷酸字符？我已经尝试过grep 和%in%，但它们返回的是空结果。这可能是我忽略的相对简单的事情。

这就是我将数据输入程序的方式。这是一个巨大的文件； 20,347 个字母，全部为 ACTG。

data <- read.table(MTHFR.txt)

我一直在尝试以这种方式将其放入字符向量中；

data.cv <- as.character(data)

但这会创建一个似乎是行号的列表，而不是核苷酸序列。

数据在online here可用。截至目前，这是数据的负责人：

head(data)
                                                                      V1
1 ATGACGATAAAGGCACGGCCTCCAACGAGACCTGTGGGCACGGCCATGTTGGGGGCGGGGCTTCCGGTCA
2 CCCGCGCCGGTGGTTTCCGCCCTGTAGGCCCGCCTCTCCAGCAACCTGACACCTGCGCCGCGCCCCTTCA
3 CTGCGTTCCCCGCCCCTGCAGCGGCCACAGTGGTGCGGCCGGCGGCCGAGCGTTCTGAGTCACCCGGGAC
4 TGGAGGGTGAGTGACGGCGAGGCCGGGGTCGCCGGGAGGGAGATCCTGGAGCCGGCAAACAACCTCCCGG
5 GGGCAAGGACGTGCTTGTGGGCGGGGAGCGCTGGAGGCCGGCCTGCCTCTCTTCTTGGGGGGGGCTGCCG
6 CCTCCCTTGCGCACCCTTCGCGGGATTAGTGTAACTCCCAATGGCTACCACTTCCAGCGACCGCCAACCC

【问题讨论】：

你不能做as.character(data) - 你需要做as.character(data$variable)
预期输出是什么？
你下载数据的时候选择了什么格式？我在您提到的页面上看到了很多选择。到目前为止，我可以比较我对 R 的有限经验和我对 python 的长期经验，我认为我宁愿使用 python（可能使用 biopython）而不是 R 来读取序列并在其中搜索模式。

标签： r regex vector character bioinformatics

【解决方案1】：

对于大多数与序列相关的生物信息学任务，您确实需要熟悉 Bioconductor 项目中一些更常见的软件包。许多常见任务已经实施了非常快的解决方案。

Biostrings 有 DNAString 和 DNAStringSet 等类，用于有效地存储和操作 DNA 字符串，AA 和 RNA 有相应的类。包括用于搜索、反向补充等的各种功能。听起来您已经导入了数据，但另一种方法是使用 readDNAStringSet() 函数。

library(Biostrings)

data <- 'ATGACGATAAAGGCACGGCCTCCAACGAGACCTGTGGGCACGGCCATGTTGGGGGCGGGGCTTCC'
dna <- DNAString(data)

matchPattern('GGG', dna)

  Views on a 65-letter DNAString subject
subject: ATGACGATAAAGGCACGGCCTCCAACGAGACCTGTGGGCACGGCCATGTTGGGGGCGGGGCTTCC
views:
    start end width
[1]    36  38     3 [GGG]
[2]    51  53     3 [GGG]
[3]    52  54     3 [GGG]
[4]    53  55     3 [GGG]
[5]    57  59     3 [GGG]
[6]    58  60     3 [GGG]

countPattern('GGG', dna)
[1] 6

countPattern('GGA', reverseComplement(dna)) #number of occurrances of 'TCC' in forward strand
[1] 2

【讨论】：

【解决方案2】：

我提出了一个解决方案，使用 biopython 读取序列并获取其反向补码，然后使用简单的算法来获取简单的已知 k-mer 的位置（如果您想要更复杂的东西，biopython 具有@987654321 的功能@)。

从文件中读取序列（假设你有 fasta 格式）：

from Bio import SeqIO
seq_record = SeqIO.read("my_sequence.fa", format="fasta")

制作正反补码的大写（以防万一）版本：

fwd = str(seq_record.seq.upper())
rev = str(seq_record.seq.reverse_complement().upper())

查找模式出现的位置（位置将在基于 0 的坐标中）：

pattern = "ACTG"
k = len(pattern)

positions_in_fwd = [i for i in range(1 + len(fwd) - k) if fwd[i:i+k] == pattern]
positions_in_rev = [i for i in range(1 + len(rev) - k) if rev[i:i+k] == pattern]

（根据您给出的序列和模式，我在序列中找到了 24 个位置，在反向补码中找到了 20 个位置。）

【讨论】：

【解决方案3】：

“寻找某些模式”有点含糊。您是否尝试提取模式？您是否试图找出它在文本中出现的间隔？我将尝试假设两者，但添加任何可以帮助指定目标的内容。

library(stringi)
library(magrittr)
# Data from site you provided was stored to "t.txt" on
# my machine so starting there
a <- readLines('t.txt')

数据信息

 > summary(a)
    Length     Class      Mode 
      292   character   character

查看数据集的头部

 > head(a,2)
    [1] "GTCAAGTTTTTTTGTTTATTTTTGAGACAGAGTCTGGCTCAATTGCCCAGGCTGAAGCAGAGGAGTGATC"
    [2] "TCAGCTCACTGCAACCTCTGCCTCCCGGGTTCAAGTGATTCTCCCGCCTCAGCTTCCTGAGTAGCTGGGA"
    > sum(nchar(a))
    [1] 20374

现在我们有了数据，让我们提取“AACG”模式

> aa <- stri_extract_all_regex(a, 'AACG', 
                         omit_no_match = F, simplify = T) %>% 
unlist %>% as.character() %>% (function(x)x[!is.na(x)])

 > aa 
[1] "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG"
[9] "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG"
[17] "AACG"

将数据集转化为一个连续的字符串：

a_flat <- paste0(a, collapse = "")

而不是提取，我们可以找到它在文本中出现的位置并变成一个数据框

bb <- as.data.frame(stri_locate_all_regex(a_flat, "AACG")[[1]])

这给我们的是模式出现的位置。

> bb
   start   end
1    807   810
2   1244  1247
3   1748  1751
4   1791  1794
5   2306  2309
6   3560  3563
7   4217  4220
8   4927  4930
9   6504  6507
10  8668  8671
11  9827  9830
12 10333 10336
13 11446 11449
14 12779 12782
15 13619 13622
16 16604 16607
17 16659 16662
18 19200 19203
19 20181 20184
20 20228 20231

我们可以使用这些位置将扁平化的字符串拆分为我们想要的内容

 > sapply(1:nrow(bb), function(i){
    stri_sub(a_flat, bb[i,'start'], bb[i,'end'])
})
[1] "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG"
[9] "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG" "AACG"
[17] "AACG" "AACG" "AACG" "AACG"

希望这对您有所启发

【讨论】：