114-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
