目的：

　　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. 输出文件如何构建