（三）R语言的生物信息学入门——Function, data.frame, 简单DNA读取与分析

注：本博客旨在分享个人学习心得，有不规范之处请多多包涵！

Function 函数

像其它的编程语言一样，R语言也有内置函数（如前面用到的c()）和自定义函数。函数一般由三个重要的部分组成：输入参数，函数主体，返回参数。R语言的函数也是允许无输入参数或返回参数的的。以下的例子为在R语言中构造与调用函数：

#用R语言内置的function()函数来声明函数，并在括号内声明输入参数。可以用=为参数设置默认值
#getDouble为构造的函数的名称
getDouble <- function(num=1){
     
  result <- num * 2 #函数主体
  return(result) #返回参数用return()
}

myResult <- getDouble(num=20) #调用函数，myResult应为40

return()和print()的区别就在与print()会在Console显示输出而return()不会。

data.frame 数据框

比起list，data.frame在计算生物学中R语言应用更为广泛。就像matrix一样，data.frame也是二维的数据结构，即由行列组成。然而，matrix的所有元素就像vector一样一定要是同一种数据类型，但data.frame只需保证每一列的数据类型相同。这就使得data.frame能够很方便地被用来读取类似于.csv或.xlsx这样的可能有多种数据类型的表格文件。data.frame的构造请见下例：

id <- c("ENO2", "TDH3", "RPL39", "GAL4") #蛋白质的代号

protName <- c("enolase", "glyceraldehyde-3-phosphate dehydrogenase", "60S ribosomal protein L39", "regulatory protein Gal4") #蛋白质的全名

abundance <- c(24563, 22369, 16232, 32.3) #细胞中的含量，单位ppm

length <- c(437, 332, 51, 881) #蛋白质中氨基酸的数量

yeastProt <- data.frame(id, protName, abundance, length, stringsAsFactors = FALSE) 
#构造data.frame；每个输入vector的长度要相同
#vector的变量名会成为表头，vector的值会按列填入data.frame

构造完成的data.frame应如下图：
data.frame的索引访问与matrix类似，具体见下例：

yeastProt[,1] #输出第1列的全部内容
yeastProt[2, ] #输出第2行的全部内容，包括表头
yeastProt[1, 2] #输出第1行第2列的元素，enolase
yeastProt[1:3, c(1, 3)] #输出第1-3行第1和3个元素，包括表头
yeastProt$length[1:2] #输出列标题为length的第1和第2个元素，437 332
yeastProt[1, "abundance"] #输出第1行列标题为abundance的元素，24563
yeastProt[yeastProt[, 4] > 500, 1] #输出yeastProt第4列大于500的行的第1列元素，GAL4

以下表格是vector, matrix, data.frame的互相转换：

	vector转为	matrix转为	data.frame转为
vector	\	as.vector()
matrix	as.matrix()	\	as.matrix()
data.frame	as.data.frame()	as.data.frame()	\

在将其它数据类型转化成data.frame时，表头可能是不正确的。修改表头请见下例：

myMatrix <- matrix(data=1:12, ncol=3, nrow=4) #3x4的matrix
myDF <- as.data.frame(myMatrix) #转为data.frame
#转化完默认的列标题为V1, V2, V3... 默认行标题为1，2，3...
colnames(myDF) <- c("c1", "c2", "c3") #自定义列标题
rownames(myDF) <- c("r1", "r2", "r3", "r4") #自定义行标题

设置完成后，matrix被转变成了如下图的data.frame：

DNA序列读取与分析实例

.fasta文件是常用的读取核酸或蛋白质序列产生的文件格式。文件中开头有一些说明信息，如来源物种和染色体编号；接下来就是以单个字母代表的核酸/蛋白质编码产生的序列。本次请从这个链接下载使用NC_045512.2版本COVID19的DNA测序数据，提取码: vfti。
请看下例用R语言进行对这一序列的分析：

install.packages('seqinr', repos='http://cran.us.r-project.org') #安装第三方库
library("seqinr") #载入第三方库
covid <- read.fasta(file="covid19.fasta") #读取.fasta文件，注意路径最好写完整。读取完是list格式
covid_seq <- covid[[1]] #把完整序列单独存储为一个vector

DNA由4种base组成：A、T、C、G。A和T配对，C和G配对。根据查格夫法则(Chargaff’s law)可推导，DNA中A和T的数量相同，C和G的数量相同。因为C和G配对后由三条氢键链接，CG组合是不容易被拆散的，也就是DNA的高CG含量会减小它变异的可能性。
计算CG含量请见下例：

countTable <- table(covid_seq) #统计每个base出现的次数
cgCount <- countTable[["c"]]+countTable[["g"]] #C和G的数量
totalCount <- length(covid_seq) #序列的总长
cgContent <- cgCount/totalCount #CG在序列中的占比

另一个比较重要的概念是DNA word，即长度大于1个base的DNA组合。通过分析DNA word是否被过表达或低表达，可以在一定程度上推测该物种演化的轨迹。$\rho$用来表示一个DNA word是否被过表达或低表达。对于一个长度为2的DNA word，计算公式如下：
$$\rho (xy)=\frac{f_{xy}}{f_x\cdot f_y}$$
其中$f_{xy}$是DNA word的出现频率除以序列总长，$f_x$和$f_y$是单个base的出现频率除以序列总长。计算某DNA word被如何表达请见下例：

baseCount <- count(covid_seq, 1) #用seqinr这个第三方库中的count()函数找长度为1的base的出现次数，和R内置的table()类似
wordCount <- count(covid_seq, 2) #找长度为2的DNA word的出现次数
totalLength <- length(covid_seq) #序列总长

#比如，我们要找TT这个DNA word的出现次数
rouTT <- (wordCount[["tt"]]/totalLength)/((baseCount[["t"]]/totalLength)*(baseCount[["t"]]/totalLength))

当$\rho$等于1时，这个DNA word为完全随机出现，与单个base出现概率的乘积相同。然而，当$\rho$大于1时，它的出现概率大于随机出现的概率，说明它被过表达，所以可推测这个DNA word在这一物种的进化方面可能有某些有益作用。反之亦然。

结束语

下次会以分子生物学和生态学的数据为例，介绍R语言中的矩阵图、条形统计图、饼图、点阵图与线性回归、带状图的使用方法，敬请期待！有任何问题或想法欢迎留言和评论！