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 </