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