1

我有一个文件（mydata.txt），其中包含许多fasta格式的外显子序列。我想为每个 DNA 序列（考虑框架）找到开始（'atg'）和停止（'taa'，'tga'，'tag'）密码子。我尝试使用matchPattern（来自BiostringsR 包的函数）来查找这些氨基酸：

例如 mydata.txt 可以是：

>a
atgaatgctaaccccaccgagtaa
>b
atgctaaccactgtcatcaatgcctaa
>c
atggcatgatgccgagaggccagaataggctaa
>d
atggtgatagctaacgtatgctag
>e
atgccatgcgaggagccggctgccattgactag

file=read.fasta(file="mydata.txt") 
matchPattern( "atg" , file)

注意：read.fasta 是seqinr包中用于导入 fasta 格式文件的函数。

但是这个命令不起作用！如何使用此功能查找每个外显子序列中的起始和终止密码子？（不移帧）

4

2 回答 2

2

我很难相信其中一个 BioC 软件包还没有做到这一点，但如果你想用基本 R 功能做到这一点，那么考虑使用gregexpr

x <- c(a='atgaatgctaaccccaccgagtaa', 
  b='atgctaaccactgtcatcaatgcctaa', 
  c='atggcatgatgccgagaggccagaataggctaa', 
  d='atggtgatagctaacgtatgctag', 
  e='atgccatgcgaggagccggctgccattgactag')

starts<-lapply(gregexpr("atg", x), function(x) x[ (x-x[1] %% 3) == 0])
stops <- mapply(function(strg, starts) {poss <- gregexpr("taa|tga|tag", strg) ; poss[[1]][ ( (poss[[1]]-starts) )%% 3 == 0]}, x, starts=unlist(starts))
 stops
#--------------
$a
[1] 22

$b
[1] 25

$c
[1]  7 31

$d
[1] 22

$e
[1] 31

您可以通过检查距离是否可被 3 整除来检查终止密码子是否在“帧内”读取：

> (25-1)%%3
[1] 0
于 2013-07-29T05:51:11.660 回答
2

'subject' 参数matchPattern是一个特殊的对象（例如 XString）。您可以通过使用 paste 折叠序列并使用?BString.

因此，使用您的数据：

file = read.fasta(file = "mydata.txt")

# find 'atg' locations
atg <- lapply(file, function(x) {
  string <- BString(paste(x, collapse = ""))
  matchPattern("atg", string)
})

atg[1:2]
# $a
#   Views on a 18-letter BString subject
# subject: atgacccccaccgagtaa
# views:
#     start end width
# [1]     1   3     3 [atg]
#
# $b
#  Views on a 21-letter BString subject
# subject: atgcccactgtcatcacctaa
# views:
#     start end width
# [1]     1   3     3 [atg]

举个简单的例子，在一个序列中查找 'atg' 的数量和位置：

sequence <- BString("atgatgccatgcccccatgcatgatatg")
result <- matchPattern("atg", sequence)
#   Views on a 28-letter BString subject
# subject: atgatgccatgcccccatgcatgatatg
# views:
#     start end width
# [1]     1   3     3 [atg]
# [2]     4   6     3 [atg]
# [3]     9  11     3 [atg]
# [4]    17  19     3 [atg]
# [5]    21  23     3 [atg]
# [6]    26  28     3 [atg]

# Find out how many 'atg's were found
length(result)
# [1] 6

# Get the start site of each 'atg'
result@ranges@start
# [1]  1  4  9 17 21 26

此外，检查?DNAString和?RNAString。它们与BString仅限于核苷酸字符相似，并允许在 DNA 和 RNA 序列之间进行快速比较。

编辑以解决评论中提到的帧移位问题：您可以使用@DWin 提到的模技巧对结果进行子集化以获得帧中的那些'atg'。

# assuming the first 'atg' sets the frame
in.frame.result <- result[(result@ranges@start - result@ranges@start[1]) %% 3 == 0]
# Views on a 28-letter DNAString subject
# subject: ATGATGCCATGCCCCCATGCATGATATG
# views:
#     start end width
# [1]     1   3     3 [ATG]
# [2]     4   6     3 [ATG]

# There are two 'atg's in frame in this result
length(in.frame.result)
# [1] 2

# With your data:
file = read.fasta(file = "mydata.txt")
atg <- lapply(file, function(x) {
  string <- BString(paste(x, collapse = ""))
  result <- matchPattern("atg", string)
  result[(result@ranges@start - result@ranges@start[1]) %% 3 == 0]
})
于 2013-07-29T06:01:39.440 回答