No.1 R语言在生物信息中的应用——序列读取及格式化输出

目的：读入序列文件（fasta格式），返回一个数据框，内容包括——存储ID、注释行(anno)、长度(len)、序列内容(content)

一、问题思考：

　　1. 如何识别注释行和序列内容行

　　2. 如何快速定位序列内容所在位置

二、你可能需要的知识——基本的R语言基础

　　1. R语言基本数据类型

　　2. 会使用帮助（help,?）及网络资源

　　3. 其他的部分可能需要你针对自己看到的问题自己想办法解决或者留言

##--构建函数--##
seq_import <- function( file ){
  seq <- readLines(file) # 读入序列，每个元素存入一行
  seq <- seq[seq != ""] # 去除空行
  is.anno <- regexpr("^>", seq, perl=T) # 正则匹配（regular expression）注释行,是注释行为1，否则为-1
  seq.anno <- seq[ which(is.anno == 1) ] # 注释内容
  seq.content <- seq[ which(is.anno == -1) ] # 序列内容
  ##--计算每条序列内容所占的行数，便于后来拼接--##
  start <- which(is.anno == 1) # 注释行行号
  end <- start[ 2:length(start) ]-1 # 第二条记录注释行到最后一条记录注释行行号减一，即为每条记录结束行号，这里会统计少一行——最后一行的结束未统计
  end <- c(end, length(seq) ) # 末尾添加一行：所有序列结束行
  distance <- end - start # 每条记录所占行号
  index <- 1:length(start) # 生成一个一到记录总个数的向量
  index <- rep(index, distance) # 分组标签
  seqs <- tapply(seq.content, index, paste, collapse="") # 拼接每条序列内容，返回一个列表，列表每个元素为一条序列的内容
  seq.content<-as.character( seqs ) # 将列表转换为向量，向量每个元素为一条序列的内容
  seq.len <- nchar(seq.content) # 获得序列长度
  seq.ID <- gsub("^>(\\w+\\|){3}([A-Za-z0-9.]+)\\|.*", "\\2", seq.anno, perl = T) # 获取序列的ID
  result <- data.frame( seq.ID, seq.anno, seq.len, seq.content ) # 组件结果：ID，长度，注释行，序列内容
  result # 最后一行作为返回值
}

三、文件内容：（复制时最后一行需要换行符，否则最后一行读取不到）

>gi|10579650|gb|AAG18645.1| hypothetical protein VNG_0001H [Halobacterium sp. NRC-1]
MTRRSRVGAGLAAIVLALAAVSAAAPIAGAQSAGSGAVSVTIGDVDVSPANPTTGTQVLITPSINNSGSA
SGSARVNEVTLRGDGLLATEDSLGRLGAGDSIEHTTHHVPLSSTFTEPGDHQLSVHVRGLNPDGSVFYVQRSVYV
TVDDRTSDVGVSARTTATNGSTDIQATITQYGTIPIKSGEHTTHLQVVSDGRIVERAPVANVSESDSANVTFDG
ASIPSGELVIRGEYTLDDEHSTHTTNTTLTYHHYHQHPQRSADVALTGVEASGGGTTYTISGDAANLGSADAASV
RVNAVGDGLSANGGYFVGKIETSEFATFDMTVQADSAVDEIPITVNYSADGQRYSDVVTVDVSGASSGSA
TSPERAPGQQQKRAPSPSNGASGGGLPLFKIGGAVAVIAIVVVVVRRWRNP
>gi|10579651|gb|AAG18646.1| amino acid ABC transporter, ATP-binding protein [Halobacterium sp. NRC-1]
MSIIELEGVVKRYETGAETVEALKGVDFSAARGEMVTVVGPSGSGKSTMLNMIGLLDSPTAGSVTLDGQD
VTGFSEDERTEERRAELGFVFQSFHLLPMLTAVENVELPSMWDTSVDRHDRAVDLLERVGLGDRLTHTPG
ELSGGQQQRVAIARSLINEPEILLADEPTGNLDQEHTTHTGGTILTEMQRLHTKHTEEENIAVVAITHDTQLEEFSDR
AVNLVDGVLHTTHH

 

四、调用函数，查看结果：

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

 

五、结果截图：

六、问题解决

　　1. 如何识别注释行和序列内容行（正则匹配）

　　2. 如何快速定位序列内容所在位置（注释行之间的就是系列，可利用此特性求出序列行所占行数，然后利用tapply+paste分组连接）

七、函数注释

1. grep(pattern, x, ignore.case = FALSE, perl = FALSE, value = FALSE,
     fixed = FALSE, useBytes = FALSE, invert = FALSE)
作用：从x搜索某种模式，找到返回下标，没找到返回数字0
pattern-匹配规则
x-待匹配对象
ignore.case-是否区分大小写
perl-
value-“否”返回下表，“是”返回匹配到的子串本身
fixed-“否”使用正则表达式匹配，“是”则不使用
useBytes-
invert-“是”反向匹配，即返回未匹配的向量的下表或本身，“否”正向匹配

> grep("w", c("a", "b", "ww", "w"), fixed = T)
[1] 3 4

2. grepl(pattern, x, ignore.case = FALSE, perl = FALSE,
      fixed = FALSE, useBytes = FALSE)
作用：与grep相似，但是返回逻辑值

> grepl("w", c("a", "b", "ww", "w"))
[1] FALSE FALSE  TRUE  TRUE

3. sub(pattern, replacement, x, ignore.case = FALSE, perl = FALSE,
    fixed = FALSE, useBytes = FALSE)
作用：在x中搜索pattern，并且以replacement替换，但是只匹配、替换一次，返回替换后的字串，没有匹配上的则不替换

> sub("(\\w)\\d",replacement = "\\1", c("a1x1", "b1x2", "c1x3") )
[1] "ax1" "bx2" "cx3"

> sub("(\\d)",replacement = "@", c("a1x1", "b1x2", "ww") )
[1] "a@x1" "b@x2" "ww" 

4. gsub(pattern, replacement, x, ignore.case = FALSE, perl = FALSE,

     fixed = FALSE, useBytes = FALSE)

作用：在x中搜索pattern，并且以replacement替换，替换所有符合要求的模式

> gsub("(\\w)\\d",replacement = "\\1", c("a1x1", "b1x2", "c1x3") )
[1] "ax" "bx" "cx"

5. regexpr(pattern, text, ignore.case = FALSE, perl = FALSE,

        fixed = FALSE, useBytes = FALSE)
作用：返回匹配到的字串的起始位置，以及匹配长度，每个元素匹配一次。未匹配到返回-1

> regexpr("t", c("temp", "tmp",  "tmppt" ), useByte = T)
[1] 1 1 1 # 匹配起始位置
attr(,"match.length")
[1] 1 1 1 # 属性值，可用att(对象, 属性)查看
attr(,"useBytes")
[1] TRUE

> tt<-regexpr("t", c("temp", "tmp",  "tmppt" ), useByte = T)
> attr(tt, "match.length")
[1] 1 1 1

6. gregexpr(pattern, text, ignore.case = FALSE, perl = FALSE,

         fixed = FALSE, useBytes = FALSE)

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

> gregexpr("t", c("temp", "tmp",  "tmppt" ), useByte = T)
[[1]]
[1] 1
attr(,"match.length")
[1] 1
attr(,"useBytes")
[1] TRUE

[[2]]
[1] 1
attr(,"match.length")
[1] 1
attr(,"useBytes")
[1] TRUE

[[3]]
[1] 1 5
attr(,"match.length")
[1] 1 1
attr(,"useBytes")
[1] TRUE

7. regexec(pattern, text, ignore.case = FALSE,

        fixed = FALSE, useBytes = FALSE)

作用：返回匹配到的字串的起始位置，以及匹配长度，每个元素匹配一次。未匹配到返回-1

> regexec("t", c("temp", "tmp",  "tmppt" ), useByte = T)
[[1]]
[1] 1
attr(,"match.length")
[1] 1
attr(,"useBytes")
[1] TRUE

[[2]]
[1] 1
attr(,"match.length")
[1] 1
attr(,"useBytes")
[1] TRUE

[[3]]
[1] 1
attr(,"match.length")
[1] 1
attr(,"useBytes")
[1] TRUE

7.tapply(X, INDEX, FUN = NULL, ..., simplify = TRUE)

作用：对向量中的数据进行分组处理（可以把每组当成一个新的向量）

X-待处理对象，通常为向量

INDEX-与X相同长度的列表，会被强制性转换为因子变量，作为分组的一句

FUN-函数

...-函数参数

simplify-“FALSE”返回结果为列表，“TRUE”返回结果为与X相同

 

> x<-c("abc", "de", "f", "gh")
> index<-c(1, 1, 2, 2)
> tapply(x, index, paste, collapse="+", simplify=F)
$`1`
[1] "abc+de"

$`2`
[1] "f+gh"

> tapply(x, index, paste, collapse="+", simplify=T)
       1        2 
"abc+de"   "f+gh" 

> paste(c("abc","de"),collapse ="+")
[1] "abc+de"