No.2 R语言在生物信息中的应用—模式匹配
目的：
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. 输出文件如何构建
