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"