Linux生信小工具–sort、uniq

豆豆写于19.5.9
整理一下linux的一些小工具，这次写一写sort与uniq这对CP

Sort 准备工作

就像cut一样，它也是对列进行处理，它默认对空字符(例如：空格、tab)进行操作，如果要指定其他的分隔符，需要用到-t

还是用例子说明问题：

新建一个文件，记录了染色体编号、起始与终止位点

$ cat >id.txt
chr1 26 39 
chr1 32 47 
chr3 11 28
chr1 40 49 
chr3 16 27 
chr1 9 28 
chr2 35 54 
chr1 10 19

Sort  实际练习

1 先按照染色体编号排，再按照起始位点排

$  sort -sk1,1 -k2,2n id.txt >id.sorted.txt
$ cat id.sorted.txt
chr1 9 28
chr1 10 19
chr1 26 39
chr1 32 47
chr1 40 49
chr2 35 54
chr3 11 28
chr3 16 27

这里用-k指定要排序的列，它表示sorting keys，每个-k 选项有两个参数代表列的范围start,end(指定起始与终止)。

如果只是排序某一列，那么只写start,start就好，就像-k1,1就表示对第一列排序；

对第一列排序后，其实依然无法确定后面的顺序，结果第二、三列还是无规则的，因此同样使用-k2,2对第二列进行排序，但是要注意，第二列是数字，因此需要增加一个-n的参数，表示numerical

会了这个，就可以很方便地对gtf文件进行排序了

2 检查文件是否为排序后的

因为sort和一般的linux工具还不同，它是同时处理多行数据进行比较，因此拿到一个标记为sorted文件，不需要重新sort，我们只需要检查一下就好

# 对一个已经排序过的文件
$ sort -k1,1 -k2,2n -c id.sorted.txt
$ echo $?
0 # 返回true的状态信息
# 如果文件没有排序过
$ sort -k1,1 -k2,2n -c id.txt
sort: id.txt:4: disorder: chr1 40 49
$ echo $?
1

3 逆向排序

使用-r参数

$ sort -k1,1 -k2,2n -r id.txt
chr3 11 28
chr3 16 27
chr2 35 54
chr1 9 28
chr1 10 19
chr1 26 39
chr1 32 47
chr1 40 49
# 可以看到，它只对第一列进行了逆向排序

# 如果要针对其他列进行逆向排序，就要把-r放在那一列的参数中
$ sort -k1,1 -k2,2nr  id.txt
chr1 40 49
chr1 32 47
chr1 26 39
chr1 10 19
chr1 9 28
chr2 35 54
chr3 16 27
chr3 11 28
# 注意这里并不是完全将第二列的数字逆向输出，因为我们还指定了-k1,1这个限制条件，因为单看chr1的第二列是逆向的，chr3亦如此
# 因此最后实现了：染色体正向输出，起始位点逆向输出

4 如何避免按照ASCII码输出

什么叫ASCII码输出结果？

简单讲，比如现在有个文件如下：

$ cat>id2.txt
chr2 15 19
chr22 32 46
chr10 31 47
chr1 34 49
chr11 6 16
chr2 17 22
chr10 30 42
$ sort -k1,1 -k2,2n id2.txt
chr1 34 49
chr10 30 42
chr10 31 47
chr11 6 16
chr2 15 19
chr2 17 22
chr22 32 46

结果发现竟然不是按照常规的chr1、chr2..这样的顺序输出，而是先chr1，然后再chr10、chr11，最后才是chr2

因为字符串中的数字，linux默认先识别第一个数字，再识别第二个数字，因此chr10、chr11要在chr2的前面出现

为了可以智能识别字符串中的数字，并按照它们排序，可以使用-V参数

$ sort -k1,1V -k2,2n id2.txt
chr1 34 49
chr2 15 19
chr2 17 22
chr10 30 42
chr10 31 47
chr11 6 16
chr22 32 46

Uniq 准备工作

uniq这个命令一般用到管道中，它接受上游的输入，然后输出的下游文件是连续的去重复的

简单的小例子

$ cat >letters.txt
A
A
B
C
B
C
C
C

Uniq 实际练习

1 直接上手

不管什么命令，拿到手第一感觉就是不去管参数，直接用试试

$ uniq  letters.txt
A
B
C
B
C

发现并没有想像中那样，直接输出A、B、C

这是因为，uniq只移除连续的重复行，比如文件中的前2行就是连续重复的，这样的只输出一次，但是像文件的3-6行的B、C、B、C，它就认为这是无重复的

但是，有解决方案：先使用sort将文件排序，再进行uniq去重复(也就是开篇说的sort与uniq是一对CP)

$ sort letters.txt | uniq
A
B
C

另外如果存在大小写区分，那么我们是可以使用-i忽略的

2 统计非重复项的次数

$  uniq -c letters.txt
2 A
1 B
1 C
1 B
3 C
$ sort letters.txt | uniq -c
2 A
2 B
4 C

3 只输出重复的行

使用-d参数即可

$ uniq -d letters.txt
A
C

4 与其他命令搭配使用

4.1 例如，要利用统计基因组各种区域的统计信息

$ grep -v "^#" gencode.v28.annotation.gff3 | cut -f3 | sort | uniq -c
735618 CDS
1237914 exon
 148007 five_prime_UTR
  58381 gene
  85439 start_codon
  77451 stop_codon
    119 stop_codon_redefined_as_selenocysteine
 144591 three_prime_UTR
 203835 transcript

如果看着数字有点混乱，还可以降序排列

$ grep -v "^#" gencode.v28.annotation.gff3 | cut -f3 | sort | uniq -c | sort -k1,1nr
1237914 exon
 735618 CDS
 203835 transcript
 148007 five_prime_UTR
 144591 three_prime_UTR
  85439 start_codon
  77451 stop_codon
  58381 gene
    119 stop_codon_redefined_as_selenocysteine

4.2 统计每条链上有多少gene feature 信息

因为sort和uniq都是基于行的，因此可以将多列组成行，然后统计(比如这里将第3列—gene feature和第7列–strand信息提取了出来)

$ grep -v "^#" gencode.v28.annotation.gff3 | cut -f3,7 | sort | uniq -c
372649 CDS    +
 362969 CDS    -
 627928 exon    +
 609986 exon    -
  74043 five_prime_UTR    +
  73964 five_prime_UTR    -
  29575 gene    +
  28806 gene    -
  43212 start_codon    +
  42227 start_codon    -
  39425 stop_codon    +
  38026 stop_codon    -
     47 stop_codon_redefined_as_selenocysteine    +
     72 stop_codon_redefined_as_selenocysteine    -
  72862 three_prime_UTR    +
  71729 three_prime_UTR    -
 103624 transcript    +
 100211 transcript    -

4.3 除此以外，还能看看某一个基因id的具体gene feature信息

$  grep ENSG00000187634.11 gencode.v28.annotation.gff3 | cut -f3 | sort | uniq -c
123 CDS
    153 exon
     21 five_prime_UTR
      1 gene
     11 start_codon
     11 stop_codon
     20 three_prime_UTR
     17 transcript

这样不用将数据导入R或者Excel，就能方便对某个或某些基因(一个流程)进行统计得到summary stats