No.2R语言在生物信息中的应用—模式匹配

OStack程序员社区-中国程序员成长平台 › 门户 › 编程› R语言›R语言教程

原作者: [db:作者] 来自: [db:来源] 收藏邀请

目的：

　　1. 计算自定义模序在所有蛋白质的匹配位点和次数

　　2. 输出超过阈值的蛋白质序列到Hit_sequences.fasta

　　3. Hit_sequences.fasta中序列用小写字母，匹配用大写字母

　　4. 返回一个数据框，内容包存储ID、注释行(anno)括——、长度(len)、匹配位置（Position），匹配次数(Hits)，相应序列（tag）

一、问题思考：

1. 如何快速计算匹配位点

2. 输出文件如何构建

>序列ID（ACCESSION）

序列内容

二、流程图

三、代码详解

 1 pattern_match<-function(pattern, sequences, hit_num){
 2   # 1. 因为字符型在数据框中被设置为因子型所以需要转换；
 3   # 2. 返回匹配起始位置，以及匹配长度（属性值：match.length）,返回值为列表
 4   pos<-gregexpr(pattern, as.character(sequences[, 4]), perl=T) 
 5   posv<-unlist(lapply(pos, paste, collapse=",")) # 把每条序列的匹配起始位置用“，”连接
 6   posv[posv == -1]<-0 
 7   fun<-function(x){ 
 8     if(x[1] == -1)
 9       0
10     else
11       length(x)
12   }
13   hitsv<-unlist(lapply(pos, fun)) # 获取每条序列匹配次数
14   sequences<-data.frame(sequences[, 1:3], Position = as.vector(posv), 
15                         Hits = hitsv, sequences[, 4]) # 构建数据框：序列id，注释，长度，Position（匹配位置），Hits(匹配次数)，序列内容
16   tag<-gsub("([A-Z])", "\\L\\1", as.character(sequences[sequences[, 5]>hit_num,6]), 
17             perl=T, ignore.case = T) # 把蛋白质序列中匹配次数大于阈值的序列转换成小写字母，这里的perl = T 为必须
18   pattern2<-paste("(", pattern, ")", sep="" ) # 重新构建模式，这样做是因为没法在模式中引入变量，变通之后就可以
19   tag<-gsub(pattern2, "\\U\\1", tag, perl=T, ignore.case=T ) # 把匹配到的模式转换为大写
20   export<-data.frame(sequences[sequences[, 5]>hit_num, -6], tag) # 构建输出数据框，大于阈值的蛋白质序列所有信息
21   selected<-export
22   # 构建写入文件的数据框格式，包括“>序列号”和序列内容
23   export<-data.frame(Acc = paste(">",export[, 1], sep = ""), seq = export[, 6]) 
24   # 先转置->转换为字符型->转换为向量（按列合并）   
25   # e.g:x<-matrix(1:4,nrow = 2, byrow = T); as.vector(x); 结果为1 3 2 4
26   write.table(as.vector(as.character(t(export))), file="Hit_sequences.fasta", quote = F, 
27               row.names = F, col.names = F)
28   cat("含有模序\"", pattern, "\"超过", hit_num, 
29       "个的所有蛋白序列已写入当前工作目录下的文件‘Hit_sequences.fasta’", "\n", seq = "")
30   cat("极度嗜盐古菌蛋白组中以下序列含有模序\"", pattern, "\"的数量超过", hit_num, "个：", "\n", seq = "")
31   print(selected[, 1:5])
32   selected
33 }

四、调用函数，查看结果（这里需要用到No.1 R语言在生物信息中的应用——序列读取及格式化输出的结果）

setwd("E:/bioinfor/bioBook/") # 设定工作目录
rm(list = ls()) # 清空变量
my_file<-"seq.txt" # 指定序列文件
source("./seq_import.R") # 载入函数
my_sequences<-seq_import(file = my_file) # 调用函数

source("./pattern_match.R") # 载入函数
hit_sequences<-pattern_match(pattern = "H..H{1,2}", sequences = my_sequences, 
                             hit_num = 2) # 调用函数

五、结果截图：

六、问题解决

　　1. 如何快速计算匹配位点

gregexpr(pattern, text, ignore.case = FALSE, perl = FALSE,
                 fixed = FALSE, useBytes = FALSE)

作用：返回匹配到的字串的起始位置，以及匹配长度（属性值），匹配所有元素的所有位置.未匹配到返回-1

　　　-> grepexpr函数可以返回匹配位点的起始位置，计算起始位置个数就可以快速计算匹配位点

　　2. 输出文件如何构建