第一，这一节讲什么？
这一节讲述一些关于R语言处理数据框的基本操作，会详细演示上一节常用函数的实际使用方法
根据我平时工作的经验，只将一些最常用的操作介绍给大家
第二，这一节如何讲？
我会提供一些示例数据，然后大家根据我的讲述，一一朝下操作
在操作过程中，我会强调自己认为一些比较重要的点
第三，这一节的目的？
演示，当获得最基本的差异分析结果后，如何快速筛选出自己的目标基因
数据格式说明
生信分析的大部分结果是表格形式，这个表格文件的格式一般为txt文档。
这些txt文档的列通常为tab制表符分割，如果使用Notepad++打开文档的话，你可能会看到类似下面的结果
不同列之间的
→即上面所说tab制表符，纯文本形式表示为\t制表符一般是不可见，同样不可见的常用字符还有每行结束时的换行符
\n, 以及空格Notepad++是一个常用的文本编辑器，在视图设置里可以将所有字符显示出来，例如在上图中，空格和制表符分别被显示为
·和→
稍微啰嗦一下，为什么生信的结果多以txt文档和csv文档(以,分割的表格)？
原因有两个：
excel文档有行数上限，对2003版最大行数是65536行，对2007以上版本，最大行数是1048676行；生信分析结果数据量较大，有时候会有限制
-
excel不是纯文本格式（主要原因），生信分析基本都在Linux操作系统下完成，而excel文件在Linux系统下是无法直接处理的
写到这里，我实在忍不住想要吐槽一下微软，这真是一个神奇的公司……生信人员必须掌握的技能表中，
dos2unix一定是名列前茅的……感觉有点莫名其妙的同学请跳过这段吐槽哈~~~
对于大多数非生信专业同学，可以鼠标右键选择excel方式打开；或者直接将文件后缀修改为xls，然后就可以直接用excel打开了。但是请大家务必注意，excel会自动将时间类的字符进行转化，例如在excel中输入Sep 2字符会自动变为2-Sep。
测试数据说明
数据存放在百度网盘test_data目录下，链接: https://pan.baidu.com/s/1iNo6_ni9mDaNzadE__ZNMg, 提取码: y3bm（上次准备的lncRNA相关数据，有的同学反映对lncRNA不是很熟悉，所以我更改为gene的数据类型）
Case-vs-Normal.xls是人细胞的差异表格示例，Case和Normal都是三个生物学重复
表格里面提供了差异比较的倍数和P值，以及每一个基因的GO和KEGG注释信息
这种表格应该是非常常见的，一般在生信公司做完基本的RNAseq分析后，都可以拿到类似的数据
那么接下来，我便给大家演示一下，如何通过R快速筛选出自己感兴趣的基因
细心的同学可能会发现，测试数据里面有一些比较奇怪的数据
有的单元格是空的，有的单元格内容为NA，这些都表示此处无内容
还有的单元格数据为Inf，这代表无穷大
数据框读取
数据框读取的时候，大家可能会比较熟悉read.table函数。但是我比较推荐read.delim函数，读取我上面说的生信格式文件会更加方便。
在下面的案例中，可以看到，read.table函数会报错。
另外一点，虽然read.delim函数默认参数sep='\t',header = T（列和列之间的分隔符，第一行是否为表头），但是我每次都会自己手动写出来。这样会增加自己对于即将处理的数据的熟悉度，算是个人一个小的工作习惯。
其实还有一个参数是
check.names=F，后面会讲
> df <- read.table('Case-vs-Normal.xls',sep='\t',header = T)
Warning messages:
1: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
EOF within quoted string
2: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
number of items read is not a multiple of the number of columns
line 10 did not have 10 elements
> df <- read.delim('Case-vs-Normal.xls',sep='\t',header = T)
为什么此处
read.table使用会出错，我暂时不做过多解释。如果有人比较纠结这个事情的话，可以尝试一下下面的代码。df <- read.table('Case-vs-Normal.xls',sep='\t',header = T, quote = '')
点击数据框预览，发现怎么有的列名变了？比如原来是FPKM_Case-1，现在变为FPKM_Case.1。大家注意一下，R在读取数据框是，会自动将中横杠-转换为点.，添加check.names=F，可以避免这种情况。
所以，在数据框读取时，我自己一般习惯的读取方式为
df <- read.delim('Case-vs-Normal.xls',sep='\t',header = T, check.names=F)
你可以再预览一下，看看是不是正常了？
还有啰嗦的一个点是，大家仔细看看一下数据框预览的结果就会发现，不同列中空值结果的表示是不一样的
在foldChange和pval两列中，无论原始数据是NA还是空的，此处都是NA
GO和GO_description两列中，空值没有任何元素
但是在KEGG列第一行中，空值表示为NA
大家记住
1.纯数字列，无论原始数据是什么形式，读取的结果，空值都是NA；字符串列，空值是和其原始数据形式统一
2.GO和GO_description两列中的空值，其实本质上空字符串，代码表示为""
3.""和NA在进行空值剔除时会有影响，这个后面碰到了再讲
数据框基本操作
快速浏览
读取完表格之后，我们通常会查看一下数据框的基本信息
1.通过summary函数查看整体情况
summary(df)
仔细看一下结果，可以发现summary返回的结果分为两种
一种是针对纯数字的列，例如foldChange列，返回值分别是该列最小值、第1四分位数、中位值、均值、第3四分位数、最大值等信息，同时还统计了空值NA的个数
在这个案例中，由于
foldChange列中存在Inf，所以最大值和均值均为无穷大
另一种是针对含有字符串的列，例如GO列，返回值为各个字符串的频数统计，并且是由大到小依次展示
在这个案例中，出现次数最多的空值，有4727次；其次是
GO:0016021，有377次
summary只显示6行结果，多余的内容都用other来代替了
2.如果想要快速浏览数据，可以通过head(df)查看数据库前6行，或者head(df, 10)查看前10行
由于结果过多，就不展示图片了
3.查看数据框维度
> dim(df) # 查看行列数
[1] 20030 13
> nrow(df) # 查看行数
[1] 20030
> ncol(df) # 查看列数
[1] 13
我一般就使用
dim（），dim(df)[1]是行数，dim(df)[2]是列数，不同记那么多函数名称了代码后面的
#表示注释内容，注释内容是不会被运行的
4.行列名查看
rownames(df) # 查看行名
colnames(df) # 查看列名
刚刚我们读取表格的时候，只用
header参数定义的列名，所以行名默认是从1开始的数字
切片、提取子集
在R里面，数据框操作的基本原则为：（我试着写着顺口溜哈）
逗号(,)间隔行和列，行的操作必有逗(,)
单个数字处理列，知道名称也能做
冒号:对应连续值，间断数据要c括(c())
列的提取样式多，数据格式不一样
1.提取某个指定元素
> df[5,6] # 提取第5行第6列元素
[1] 0
> df[8,9] # 提取第8行第9列元素
[1] 0.9832643
2.提取整行数据
#### 提取单行数据
df[2,] # 提取第2行
##### 提取多行数据
df[1:4,] # 提取1-4行
df[c(2,4,5),] # 提取2、4、5行
3.提取整列数据
对于数据框而言，提取列的方式有多种
#### 按照位置信息提取列
df[,1] # 提取第1列
df[1] # 提取第1列
#### 按照列名提取列
df['gene_id'] #提取gene_id列
df$gene_id #提取gene_id列
#### 提取多列
df[2:6] # 提取2-6列
df[c(5,8,10)] # 提取5、8、10列
但是上面不同方式提取出来的结果，其数据类型有所区别，可以通过class函数来查看相应数据类型
> class(df[1])
[1] "data.frame"
> class(df[,1])
[1] "factor"
> class(df['gene_id'])
[1] "data.frame"
> class(df$gene_id)
[1] "factor"
对于普通用户，仅进行简单的数据处理或绘图，上面几种方法可以随便使用，
如果以后有高级分析需求，例如做差异分析等复杂处理，需要注意一下数据结构的问题
这个我们以后碰到了再说
4.提取子集
df[1:4,2:6] #提取第1-4行,2-6列
df[c(1,5,10),c(2,9,10)] # 提取第1、5、10行,第2、9、10列
df[1:4,c('gene_id','pval')] # 提取'gene_id'和’pval'两列的第1-4行
数据处理
数据框基本操作还有很多，就不过多讲解了，
下面直接进入实际操作——数据筛选，
将其他常见的数据框操作融合进实际案例处理，
这样能够让大家加深理解
去除NA
前面提到过，foldChange和pval两列中有NA值
这说明差异比较的时候，Normal的表达值为0(分母为0，没有意义)
所以我们首先要将这些基因剔除了，可以通过na.omit()函数处理
> dim(df)
[1] 20030 13
> dim(na.omit(df)) # 剔除含有NA的整行内容
[1] 17161 13
注意，na.omit（）只剔除含有NA的行
上面提到那些GO和GO_description等列，虽然有空值，但是na.omit()无法处理这些数据
df_delNA <- na.omit(df)
head(df_delNA) # 可以看到，GO和GO_description等列还是有空值
如果执行
df <- na.omit(df)的话，会直接替换df数据本身通常，我会将处理过后的数据保存进新的变量，如
df_delNA，这样的好处是不会改变原始数据，如果那里有问题了，便于重新处理数据
条件筛选
直接筛选
剔除完空行之后，我们希望获得差异倍数较大（foldChange绝对值大于等于2），且极显著(pval小于等于0.01)的基因
首先，进行显著性筛选
> head(df_delNA$pval <= 0.01)
[1] FALSE FALSE FALSE FALSE FALSE FALSE
> dim(df_delNA[df_delNA$pval <= 0.01,])
[1] 145 13
df_delNA$pval <= 0.01会将该列的数值同0.01一一比较，符合条件的话返回TRUE，否则返回FALSE
这些TRUE或FALSE组成的列表对应不同行，可以将其作为行的筛选条件
最后，我们还是将筛选出来的结果保存为新的变量
df_delNA_0.01 <- df_delNA[df_delNA$pval <= 0.01,]
然后，对差异倍数进行筛选
筛选条件foldChange绝对值大于等于2，可以拆分为foldChange大于等于2或foldChange小于等于-2
用代码表示就是
df_delNA_0.01$foldChange >=2 | df_delNA_0.01$foldChange <=-2
在R里面，多个条件通过
&(并且)和|(或)来表示
另外一种表示方式，可以直接用绝对值函数abs()来处理
abs(df_delNA_0.01$foldChange) >=2
这两种的结果是一样的
> dim(df_delNA_0.01[df_delNA_0.01$foldChange >= 2 |
+ df_delNA_0.01$foldChange <= -2, ])
[1] 58 13
>
> dim(df_delNA_0.01[abs(df_delNA_0.01$foldChange) >= 2, ])
[1] 58 13
1.在R里面撰写代码的时候，如果一条命令太长了，可以将其分为多行来撰写
虽然
df_delNA_0.01$foldChange >= 2 | df_delNA_0.01$foldChange <= -2分为两行撰写了，但是中间
|会将其连接为一个整体分行撰写会增加代码的可读性
2.通过
dim(df_delNA_0.01[abs(df_delNA_0.01$foldChange) >= 2, ])可以看出在R里面，函数是可以直接叠加使用的
上面的代码逻辑是
先
abs()求该列的绝对值，然后判断绝对值是否大于2其次将判断结果作为行筛选条件赋给数据框，
最后
dim()计算最终筛选结果的行列数
先处理，再筛选
上面展示的是，我们根据表格中的原始数据，设置筛选条件，缩减表格
但是也有很多情况是，我们需要先对原始数据进行一定的理，然后根据处理后的结果，再进行筛选
比如我们希望拿到实验组中低丰度(FPKM小于1)，但是差异却极显著(pval小于0.01)的基因
那么首先，我们需要计算实验组的FPKM均值
实验组的数据在哪里呢？
> colnames(df_delNA)
[1] "gene_id" "FPKM_Case-1" "FPKM_Case-2"
[4] "FPKM_Case-3" "FPKM_Normal-1" "FPKM_Normal-2"
[7] "FPKM_Normal-3" "foldChange" "pval"
[10] "GO" "GO_description" "KEGG"
[13] "KEGG_description" "caseMean"
数一下，实验组(Case)位于2-4列，所以计算他们的均值就是
df_delNA$caseMean <- rowMeans(df_delNA[2:4])
1.
df_delNA[2:4]直接提取数据框df_delNA第2-4列但是我一般会写成
df_delNA[,2:4]这样，我会很清楚的知道，我在操作列
2.
rowMeans函数可以计算数据框的行均值，对应的colMeans可以计算列均值类似的函数还有
rowSums、colSums等3.
df_delNA$caseMean会在数据框df_delNA中直接添加新列caseMean> colnames(df_delNA) [1] "gene_id" "FPKM_Case-1" "FPKM_Case-2" [4] "FPKM_Case-3" "FPKM_Normal-1" "FPKM_Normal-2" [7] "FPKM_Normal-3" "foldChange" "pval" [10] "GO" "GO_description" "KEGG" [13] "KEGG_description" "caseMean"
聪明的同学可能已经发现，既然实验组的名称里面有Case关键字，那么是否有快捷的方法？
上一节常用函数里面，我们提过一个grep函数，大家可以试一下
df_delNA$caseMean <- rowMeans(df_delNA[,grep('Case',colnames(df_delNA))])
grep可以进行关键词搜索，返回匹配的位置信息如果忘了的同学可以回顾一下上一节的内容
> grep('Case',colnames(df_delNA)) [1] 2 3 4
接下来筛选就和上面讲的内容差不多了
df_delNA_lowEx <- df_delNA[df_delNA$caseMean <=1, ]
df_delNA_lowEx_sig <- df_delNA_lowEx[df_delNA_lowEx$pval <= 0.01,]
最终，我们筛选到63个在实验组中表达丰度低，但却具有显著调控作用的基因
> dim(df_delNA_lowEx_sig)
[1] 63 14
