学习R语言里的排序函数

我们的生信入门班和数据挖掘线上直播课程已经有了三年多的历史，培养了一波又一波优秀的生信人才。新一期的生信入门班学员热情高涨，他们富有探索精神，又乐于助人，这是一个在帮助别人的过程中自己也获得进步的例子，令我们非常开心。

学习R语言里的排序函数

(生信技能树学员徐谦)

正在上5月生信入门课程的小伙伴们应该初步掌握了一些R语言的基本函数，其中有一类函数可以称为排序函数，例如周二细讲的sort，order函数，以及不常使用的rank函数，这里说的只是R基础包中的排序函数，另外还有其他包中的一些更方便的函数，例如dplyr包中的arrange，这些后续我们再共同学习。先总结一下这几个函数的使用：

1.sort、order、rank

都是接受一个对象，通常为向量，运行之后根据这个向量中的元素位置，返回不同的值：

例如我们将x赋值如下，分别运行三个函数：

> x <- c(3,6,4,8,2)
> order(x)
[1] 5 1 3 2 4
> sort(x)
[1] 2 3 4 6 8
> rank(x)
[1] 2 4 3 5 1

这里很容易理解，order返回的是x当中元素的位置即下标，默认从小到大。这里可以补充一点，在R语言当中，我们看到的向量实际上是有两组信息的，一组是向量当中的元素，我们可以理解为萝卜，另外一组就是他的下标，我们可以理解为萝卜对应的坑。这里即x当中最小的萝卜是2，对应第5个坑。所以结果是51324。

sort呢，sort即将萝卜都拔出来，重新排序，默认从小到大，所以返回了23468。

那rank呢，rank是挨个比萝卜，返回值是它应该种的坑，默认还是从小到大，比如x里面第一个萝卜是3，它在x所有萝卜中排老2，所以先返回2，第二个萝卜6排老4，所以返回4。

总结：order对萝卜进行排序，但返回的是坑的位置；sort是将萝卜拔出来重新排，直接了当；rank是按现有顺序挨个比萝卜大小，返回它该种的坑的位置。

这里说的都是数字，那么其他类型的向量呢？

2. 其他类型向量的排序：

> x <- c("B","A","D","F","C")
> order(x)
[1] 2 1 5 3 4
> sort(x
[1] "A" "B" "C" "D" "F"
> rank(x)
[1] 2 1 4 5 3
> x <- c("张三","李四","王五","赵六","朱七")
> order(x)
[1] 2 3 1 4 5
> sort(x)
[1] "李四" "王五" "张三" "赵六" "朱七"
> rank(x)
[1] 3 1 2 4 5

同样是可以排序的，只不过字符串会按照字母或者汉字拼音首字母进行排序，三个函数返回的值就不需要再赘述了。这里需要强调一点，如果一个数字向量，排序的时候会按数字大小排序，但是当数字和字母在一起组成一起，就不一定会按照数字大小排序了。举个例子：

> x <- c("chr1","chr3","chr2","chr23","chr21")
> order(x)
[1] 1 3 5 4 2
> sort(x)
[1] "chr1"  "chr2"  "chr21" "chr23" "chr3" 
> rank(x)
[1] 1 5 2 4 3

可以理解了，字符排序，是有优先级的，不同染色体的名称，chr2后面不管有没有其他东西，有多少，它一定是在chr3前面的。

3. 参数的更改

R语言里所有的函数都是有参数的，我们可以根据函数作者的设定，赋予不同的参数，例如查阅帮助文档，可以看到sort，order都可以设定decreasing = T或者F来控制顺序，这个大家都知道了，就不重复了。

4. 元素重复了怎么排

如果向量x当中有重复元素呢？再举个例子试一下就可以了：

> x <- c(3,6,2,4,8,2);x
[1] 3 6 2 4 8 2
> order(x,decreasing = F)
[1] 3 6 1 4 2 5
> order(x,decreasing = T)
[1] 5 2 4 1 3 6

两个一样的萝卜，如果是从小到大，先返回先出现的萝卜的坑，如果从大到小，就先返回后出现萝卜的坑。sort返回的是萝卜，两个一样的萝卜谁先谁后没区别。那rank呢，严格来说rank返回的是秩次，这个学过秩和检验的都知道，一样的元素，实际上秩次也是一样的，具体算法大家一看便知：

> x <- c(3,6,2,4,8,2);x
[1] 3 6 2 4 8 2
> rank(x)
[1] 3.0 5.0 1.5 4.0 6.0 1.5
> x <- c(3,6,2,4,8,2,2);x
[1] 3 6 2 4 8 2 2
> rank(x)
[1] 4 6 2 5 7 2 2

对了，返回的是坑的平均值。前两个元素2占坑1和2，平均1.5。如果前三个元素都是2，占坑123，平均还是2。

事情到这里应该就结束了，直到有个小伙伴在群里发了个截图：

我的第一反应是order只会接受一个向量，如果给他2个向量，他会选择性忽略第二个（我把order(x)和order(x,y)的结果看成一样的了）。但是后来一想不太对，R语言中几乎所有的函数都是有严格的对象和参数要求的，如果给了它函数里没写的东西，那大部分时候就会报错，如果没报错，那就是函数接受了，当然也有其他特殊例外的情况。随后我就查了一下帮助文档，来了这么一段示例：

突然有点懵逼，试着自己运行了一下，才恍然大悟，实际上它是可以种好几排萝卜的，理论上可以无限种，但是有一个前提，那就是——坑得一样多，也就是说不管种几排，坑的顺序是不变的。那order(x,y)会返回啥呢？？？我们赋值试一下，只种两排萝卜：

> x <- c(1,3:1);x
[1] 1 3 2 1
> y <- c(9:6);y
[1] 9 8 7 6
> order(x)
[1] 1 4 3 2
> order(x,y)
[1] 4 1 3 2

x里有两个1，为什么把order(x)把坑1排坑4的前面，而order(x,y)，把坑4排的坑1前面，聪明的你应该已经看出来了，刚才我们说过order(x)中如果有一样的萝卜，默认是先出现的萝卜对应的坑排前面，所以1在4前面。而有了另外一排萝卜y后，如果x中有一样的萝卜，它会比y里萝卜的大小，x里有两个萝卜1，对应的y里的萝卜是9和6，6比9小，所以坑4排在坑1前面了。

5. 排序有什么用？

以上是R语言中基础函数中几个排序函数的用法，那排序到底有什么用呢？实际上在R语言中我个人觉得order比sort用的多，原因就是他会返回坑的位置。我们用坑的位置可以做很多事情，因为我们经常操作的数据框中，每一列都是一个向量，每一列都有一样顺序的坑，有了坑的位置我们就可以按行来提取数据框了，就可以按照某一列萝卜的顺序对行进行排序，类似于Excel中按列排序或者筛选扩展到其他列。同样道理，我们也可以用TRUE或者FALSE对坑进行标记，从而扩展到行。基础函数是可以实现数据的筛选的，但是很多时候会有更好的替代方法，等我们再往后接触到dplyr包，apply家族以及正则表达式的时候，就会发现熟练的数据处理是一件能提高效率提高产出的事情。这里可以举一个很简单的例子，就不举数据框的处理了，我们直接举向量的例子。R语言中向量是可以命名的，例如如下代码：

> x <- c(3,6,4,8,2);x
[1] 3 6 4 8 2
> y <- c("张三","李四","王五","赵六","朱七");y
[1] "张三" "李四" "王五" "赵六" "朱七"
> names(x) <- y;x
  张三 李四 王五 赵六 朱七 
   3    6    4    8    2 

好了现在向量x有了三组信息，萝卜，坑的顺序和萝卜的名字，我们继续运行排序函数试试：

> order(x)
[1] 5 1 3 2 4
> sort(x)
朱七 张三 王五 李四 赵六 
   2    3    4    6    8 
> rank(x)
 张三 李四 王五 赵六 朱七 
   2    4    3    5    1 

可以看到order和rank之后把萝卜的名字也给出了，是不是很方便知道同学的名次了。后面还会学到一种富集分析叫GSEA，这是更具体的应用。如果用Y叔的神包clusterProfiler做GSEA的话，需要输入的对象很简单，就是一组排过序的数字向量，这里可以是logFC或者相关系数（具体后面会知道）。

> geneList <- gene_df$logFC
> geneList
   [1] -0.3885455532  0.2467344071  0.0903026505  0.1419874626  0.1666925658
   [6] -0.0332339569 -0.1111574881  0.1920079524  0.2683369453  0.1834692777
  [11]  0.1413467028  0.0139812991  0.2666602717  0.0672510451  0.1784876312
  [16] -0.4048927912  0.0491033960 -0.0391377292 -0.1960619184  0.0295863678
  [21]  0.0668118973 -0.4134487165 -0.0695667235  0.1275862624  0.1677585023
  [26] -0.0628605909  0.0342075099 -0.0218328844 -0.6710892754 -0.1621723730
#后面还有很多

如果我们直接提取logFC，可以看到logFC有高有低有正有负，这时候如果直接做GSEA是不行的，GSEA要求输入对象是排过序的向量，我们就可以用到向量的命名和sort：

> names(geneList) = gene_df$SYMBOL
> geneList = sort(geneList, decreasing = TRUE)
> geneList
         A4GALT            AAAS            AACS           AADAT           AAGAB 
      1.4274086       1.2361003       1.0576695       1.0375003       1.0256938 
           AAK1           AAMDC            AAMP            AAR2           AARS2 
      0.9780116       0.9759169       0.9236183       0.9226613       0.9015509 
         AARSD1           AASDH        AASDHPPT            AASS            AATF 
      0.8799943       0.8419719       0.8249711       0.8224399       0.8160568 
          ABCA1          ABCA10           ABCA2           ABCA5           ABCA7 
      0.8154475       0.8109042       0.8024702       0.7867225       0.7807610
#后面还有很多

向量的名字是基因名，向量的按logFC大小排序，这时候做把geneList这个向量传进去就可以做GESA分析，分析完就可以直接出图了：

> hallmarks <- read.gmt("resource/h.all.v7.5.1.symbols.gmt")
> y <- GSEA(geneList,TERM2GENE =hallmarks)
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done...
There were 11 warnings (use warnings() to see them)
> yd <- as.data.frame(y)
> library(ggplot2)
> dotplot(y,showCategory=6,split=".sign")+facet_grid(~.sign)

是不是挺有用了，这只是一个简单的小例子，在此是想说明对大部分人来说，做生信数据挖掘的过程中，出图从来不是最重要的事，不是GSEA怎么做，dotplot怎么画。最重要的是数据的清洗、想法的获得，是geneList这个向量或者其他要给的对象它怎么来的。我现在就是卡在这才回来补课。tidyverse那些包，apply函数家族的运用，正则表达式这些，都很考验R基础。有了一个好的想法，数据整理好了，出图就是像上面那样分分钟的事儿。

最后，希望各位小伙伴在生信技能树团队的带领下，打好基础，熟练运用各种数据处理函数，早日看到自己的学习成果，完成属于自己的生信分析。与君共勉！