linux极简统计分析工具datamash教程

linux极简统计分析工具datamash教程

18 min read

引子

之前写 awk 教程的时候,曾经提到过一些对文本中行列进行某些计算统计的需求,例如使用数组分类求和。一些基本需求 awk 都可以实现,但是写起来稍显复杂。在 R 中使用 dplyr 或者基础函数 aggregate() 可以方便的进行分组操作,如果能在 linux 中使用更加简洁的单行命令针对数值和字符进行一些基本运算就省去了在 R 终端操作的时间。这篇文章介绍一个 linux 中能满足这类需求的工具 GNU datamash

两个使用 awk 的例子如下:

cat  awk2.txt
USSR    8649    275     Asia
Canada  3852    25      North_America
China   3705    1032    Asia
USA     3615    237     North_America
Brazil  3286    134     South_America
India   1267    746     Asia
awk '{ pop[$4] += $3 }END{ for (name in pop) print name,pop[name]}' awk2.txt

#North_America 262
#Asia 2053
#South_America 134

例如利用getline通过第一列信息判断,将第一列相同的第二列内容合并。

cat awk3.txt

a       qw
a       we
a       wet
b       wer
b       klj
b       piu
c       eie
c       tmp
c       ike

awk 'BEGIN{getline;a=$1;printf ("%s\t%s",$1,$2)}{if(a==$1){printf ","$2}else{printf "\n%s\t%s",$1,$2;a=$1}}END{printf "\n"}' awk3.txt

a       qw,we,wet
b       wer,klj,piu
c       eie,tmp,ike

datamash 是什么

在 gnu 官网 中,用了下面一句话来介绍 datamash。

GNU datamash is a command-line program which performs basic numeric, textual and statistical operations on input textual data files.

datamash 作为一个命令行程序可以对文本文件进行数字和文本相关的基本统计操作(虽说基本但是所有的操作都足够常用高频)。

使用前首先进行安装,可以通过sudo apt-get install datamash来进行安装或者通过源码来安装最新版本,这里采用第二种方式。

wget https://ftp.gnu.org/gnu/datamash/datamash-1.4.tar.gz
tar -xzf datamash-1.4.tar.gz
cd datamash-1.4
./configure
make
make check
sudo make install

调用格式与参数

datamash 的基本调用格式如下:

datamash [option]… op1 column1  [op2 column2 …]

上面的内容转换为描述语言就是:在option的参数下,对column1列进行op1操作,对column2列进行op2操作。如果使用--group 参数,所有的operation都将会分组进行;如果没有使用--group,所有的operation会对所有的值进行。需要说明的是这里的column1 可以是表示第几列的数字,当使用-H或者--header-in时可以是所选字段的名称,可以使用列名,当 operation 要求输入成对数据的时候使用: 连接,比如spearson 5:6

主要 operation

不同版本支持的参数不同,以下参数适用于 v1.4 版本。在分组统计相关的参数中,p/s 前缀分别代表 population 或者 sample。一般而言,sample 对应的计算等同于 R 中对应函数,例如sstdevsd()是一致的。

Primary operations:

primary operations 决定了文件将被如何处理,如果 primary operatino 没有设置,整个文件将会被逐行处理(针对 per-line operation)或者将所有行按照一组进行处理。

Line-Filtering operations:

Per-Line operations:

Group-by Numeric operations:

Group-by Textual/Numeric operations:

Group-by Statistical operations:

主要参数

操作示例

使用软件自带的测试数据集 scores_h.txt,有三列,分别是学生姓名科目和成绩。

cat /usr/local/share/datamash/examples/scores_h.txt |head
Name    Major   Score
Shawn   Arts    65
Marques Arts    58
Fernando        Arts    78
Paul    Arts    63
Walter  Arts    75
Derek   Arts    60
Nathaniel       Arts    88
Tyreque Arts    74
Trevon  Arts    74

文本基本操作

#检查行列完整性
$ cat scores_h.txt |datamash check

84 lines, 3 fields

# 交换列顺序
$ cat scores_h.txt |datamash reverse |head -n5

Score   Major   Name
65      Arts    Shawn
58      Arts    Marques
78      Arts    Fernando
63      Arts    Paul

# 根据前两列构建透视表,且显示分数
$ cat scores_h.txt |datamash --header-in crosstab 1,2 unique 3 |head

# 因为这个文件有 header 所以需要使用 --header-in 参数
    Arts    Business        Engineering     Health-Medicine Life-Sciences   Social-Sciences
Aaron   N/A     83      N/A     N/A     58      N/A
Allen   N/A     N/A     N/A     N/A     50      N/A
Andre   N/A     N/A     N/A     72      N/A     N/A
Angel   N/A     N/A     N/A     100     N/A     N/A
Anthony N/A     N/A     N/A     N/A     32      N/A
Antonio N/A     N/A     88      N/A     N/A     N/A
Austin  N/A     N/A     N/A     N/A     91      N/A
Avery   N/A     N/A     51      N/A     N/A     N/A
Brandon N/A     N/A     N/A     N/A     72      N/A

分组统计操作

统计每门课的选课人数

$ datamash --header-in --sort groupby 2 count 2 < scores_h.txt
Arts    19
Business        11
Engineering     13
Health-Medicine 13
Life-Sciences   12
Social-Sciences 15

查看每门课程的最高分和最低分,在这里使用--header-out 输出时会自动添加列名

$ datamash --header-in --header-out --sort groupby 2  min 3 max 3 < scores_h.txt
GroupBy(Major)  min(Score)      max(Score)
Arts    46      88
Business        79      94
Engineering     39      99
Health-Medicine 72      100
Life-Sciences   14      91
Social-Sciences 27      90

因为这是有列名的输入文件,使用--header-in参数后可以使用列名来进行计算。统计每门课程的平局分和标准差,在显示的结果中要求只保留两位小数即可。

$ datamash --header-in --header-out --sort -R 2 \
groupby 2  mean Score pstdev Score < scores_h.txt

GroupBy(Major)  mean(Score)     pstdev(Score)
Arts    68.95   10.14
Business        87.36   4.94
Engineering     66.54   19.10
Health-Medicine 90.62   8.86
Life-Sciences   55.33   19.73
Social-Sciences 60.27   16.64

如果想要对文本中的数字进行简化,可以使用如下几个命令:

$ (echo num; seq -1.25 0.25 1.25) \
|datamash --full -H round 1 ceil 1 floor 1 trunc 1 frac 1

num     round(num)      ceil(num)       floor(num)      trunc(num)      frac(num)
-1.25   -1      -1      -2      -1      -0.25
-1.00   -1      -1      -1      -1      0
-0.75   -1      0       -1      0       -0.75
-0.50   -1      0       -1      0       -0.5
-0.25   0       0       -1      0       -0.25
0.00    0       0       0       0       0
0.25    0       1       0       0       0.25
0.50    1       1       0       0       0.5
0.75    1       1       0       0       0.75
1.00    1       1       1       1       0
1.25    1       2       1       1       0.25

在实际的服务器管理中,可能需要统计一些系统用户信息。这些信息存放在/etc/passwd中,:分割。这里使用login shells 进行分组,然后统计用户数。

$ cat /etc/passwd | datamash -t : --output-delimiter $'\t' --sort groupby 7 count 7

/bin/bash       10
/bin/false      22
/bin/sh 3
/bin/sync       1
/usr/sbin/nologin       17

生信相关文本

基因表达矩阵

假设有一个行为基因列为组织的表达矩阵如下所示

$ cut -f1-5 input.matrix |head
geneid  Flag    Leaf    Node    Sheath
G1011   0       0.11069839039242035     4.578044680376127       2.202828823551847
G1031   0.09338340723174347     0.11069839039242035     0.46883590100237443     0
G5133   0       0       0.05515716482380876     0
G3987   2.5213519952570738      3.09955493098777        583.17670368213 16.598059042576708
G5604   0.09338340723174347     0       2.923329735661864       0.051228577291903415
G5664   0       0.11069839039242035     0.6067288130618963      0
G4107   0       0.11069839039242035     4.21952310902137        0.15368573187571025
G4147   0.09338340723174347     0       0.19305007688333065     0.051228577291903415
G5786   0       0.33209517117726106     5.1296163286142145      0.9221143912542615

我想要计算每个基因在所有组织中的表达情况(均值标准差),可以进行如下操作

$ cat input.matrix |datamash transpose \
|datamash --header-in --header-out  mean 2-10 sstdev 2-10 |datamash transpose

mean(G1011)     3.4218560658562
mean(G1031)     0.13477809345566
mean(G5133)     0.026272278033758
mean(G3987)     237.63686611126
mean(G5604)     1.6544141646051
mean(G5664)     0.25398335468086
mean(G4107)     0.92018619120978
mean(G4147)     0.074272982619755
mean(G5786)     1.8763718749021
sstdev(G1011)   3.6565336756012
sstdev(G1031)   0.15591791251859
sstdev(G5133)   0.025561380705389
sstdev(G3987)   276.34845044642
sstdev(G5604)   1.9673226520581
sstdev(G5664)   0.23076182316703
sstdev(G4107)   1.5120321659008
sstdev(G4147)   0.063632028432619
sstdev(G5786)   2.1600115292083

但是上面的格式看起来不是很美观,可以利用 datamash 自身的命令稍加调整

$ cat input.matrix |datamash transpose \
|datamash --header-in --header-out  sstdev 2-10 mean 2-10 \
|datamash transpose |sed 's/(/\t/;s/)//' \
|datamash crosstab 1,2 unique 3 |datamash transpose

    mean    sstdev
G1011   3.4218560658562 3.6565336756012
G1031   0.13477809345566        0.15591791251859
G3987   237.63686611126 276.34845044642
G4107   0.92018619120978        1.5120321659008
G4147   0.074272982619755       0.063632028432619
G5133   0.026272278033758       0.025561380705389
G5604   1.6544141646051 1.9673226520581
G5664   0.25398335468086        0.23076182316703
G5786   1.8763718749021 2.1600115292083

GenePred 文件操作

genepred 和 gff3 文件格式信息一致,但是更加便于我们进行操作。genepred 文件有两种格式,其中一种是标准格式,一种是 extended 格式,具体解释如下。

table genePred
"A gene prediction."
    (
    string  name;               "Name of gene"
    string  chrom;              "Chromosome name"
    char[1] strand;             "+ or - for strand"
    uint    txStart;            "Transcription start position"
    uint    txEnd;              "Transcription end position"
    uint    cdsStart;           "Coding region start"
    uint    cdsEnd;             "Coding region end"
    uint    exonCount;          "Number of exons"
    uint[exonCount] exonStarts; "Exon start positions"
    uint[exonCount] exonEnds;   "Exon end positions"
    )

table genePredExt
"A gene prediction with some additional info."
    (
    string name;        	"Name of gene (usually transcript_id from GTF)"
    string chrom;       	"Chromosome name"
    char[1] strand;     	"+ or - for strand"
    uint txStart;       	"Transcription start position"
    uint txEnd;         	"Transcription end position"
    uint cdsStart;      	"Coding region start"
    uint cdsEnd;        	"Coding region end"
    uint exonCount;     	"Number of exons"
    uint[exonCount] exonStarts; "Exon start positions"
    uint[exonCount] exonEnds;   "Exon end positions"
    int score;            	"Score"
    string name2;       	"Alternate name (e.g. gene_id from GTF)"
    string cdsStartStat; 	"Status of CDS start annotation (none, unknown, incomplete, or complete)"
    string cdsEndStat;   	"Status of CDS end annotation (none, unknown, incomplete, or complete)"
    lstring exonFrames; 	"Exon frame offsets {0,1,2}"
    )

首先检查一下我们的数据情况,可以看出是 extended 格式

$ datamash check < test.genepred
133744 lines, 15 fields

统计每个基因的转录本个数,同时输出这些转录本的 id

cat test.genepred|datamash -s -g 12 count 1 collapse 1

看看每个染色体有多少基因

$ cat test.genepred|datamash -s -g 2 countunique 12  |head
chr1A_part1	2694
chr1A_part2	1650
chr1B_part1	2428
chr1B_part2	2286
chr1D_part1	3698
chr1D_part2	786
chr2A_part1	2712
chr2A_part2	3085
chr2B_part1	3076
chr2B_part2	3033

如果结合awk还可以进行很多筛选。例如统计每个基因转录本的个数,这些转录本中外显子的最大值最小值,然后挑选出多余 3 个转录本的基因。

$ cat test.genepred|datamash -s -g 12 count 8 min 8 max 8 |awk '$2>2' |head
test1A02G001900	10	9	10
test1A02G005600	6	12	13
test1A02G009000	3	11	13
test1A02G012500	4	13	14
test1A02G013800	4	12	13
test1A02G016900	3	8	12
test1A02G030900	3	6	8
test1A02G035500	4	3	5
test1A02G039100	5	5	9
test1A02G047200	3	10	12

除此以外,针对 bed 文件也有非常高频的使用场景这里就不再展示了。

其它 one-liner

说到使用单行命令解决问题,通常都是 perl 的骄傲,当然 awk 也绝对是一把好手。如果特别喜欢用 R,其实命令行里也可以通过Rscript -e来强行单行。以下简单举一些例子进行比较。

awk 和 datamash

awk 计算 mean min 和 max 都没有压力,使用 END 可以解决问题。

$ seq 10 | datamash sum 1
55
$ seq 10 | awk '{sum+=$1} END {print sum}'
55

$ seq -5 1 7 | datamash min 1
-5
$ seq -5 1 7 | awk 'NR==1 {min=$1} NR>1 && $1<min { min=$1 } END {print min}'
-5

$ seq 10 | datamash mean 1
5.5
$ seq 10 | awk '{sum+=$1} END {print sum/NR}'
5.5

如果要进行分组操作,awk 就需要数组和 getline 来搞定。有时候需要写很长一行命令,并不是很划算。

首先构造一个两个的矩阵:

printf "%s\t%d\n" a 1 b 2 a 3 b 4 a 3 a 6 > test.txt

awk 根据第一列的 key 进行分组统计个数以及求和,可以看出awk 需要写的各种括号实在是不少。

# 统计个数
$ cat test.txt |datamash -s -g 1 count 2
a       6
b       4
$ cat test.txt |awk '{a[$1]++} END {for(i in a) { print i, a[i] }}'
a 6
b 4

#求和
$ cat test.txt| datamash -s -g 1 sum 2
a	13
b	6

$ cat test.txt |awk '{a[$1]+=$2} END {for(i in a) { print i, a[i] }}'
a 13
b 6

R 和 datamash

尽管可能 90% 以上的人说到单行命令都不会想到 R,但 R 确实可以强行单行。不过要单行操作用到的更多是 R 基础函数,很多喜闻乐见的包用起来就不方便了。

当统计基本的描述统计量时,以下为两者的对比,此时感觉 R 的 summary() 更方便。

$ cat test.txt| datamash --header-out  min 2 q1 2 median 2 mean 2 q3 2 max 2
min(field-2)	q1(field-2)	median(field-2)	mean(field-2)	q3(field-2)	max(field-2)
1	2.25	3	3.1666666666667	3.75	6

$ cat test.txt| Rscript -e 'summary(read.table("stdin"))'
 V1          V2
 a:4   Min.   :1.000
 b:2   1st Qu.:2.250
       Median :3.000
       Mean   :3.167
       3rd Qu.:3.750
       Max.   :6.000

如果要分组计算,R 可以使用aggregate

$ cat test.txt|datamash -s --header-out -g 1 min 2 q1 2 median 2 mean 2 q3 2 max 2
GroupBy(field-1)	min(field-2)	q1(field-2)	median(field-2)	mean(field-2)	q3(field-2)	max(field-2)
a	1	2.5	3	3.25	3.75	6
b	2	2.5	3	3	3.5	4

$ cat test.txt|Rscript -e 'a=read.table("stdin")' -e 'aggregate(a$V2,by=list(a$V1),summary)'
  Group.1 x.Min. x.1st Qu. x.Median x.Mean x.3rd Qu. x.Max.
1       a   1.00      2.50     3.00   3.25      3.75   6.00
2       b   2.00      2.50     3.00   3.00      3.50   4.00

总之,在 Linux 下有这么好用的一个工具,命令简单且可以轻松通过管道符和 awk 以及 sed 等连用,完成平时分析数据时的很多基本需求。推荐学习和使用。

参考网站:


本文作者:思考问题的熊

版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

扫码即刻沟通

前一篇

每周文献-190622-ATAC-seq和转录本从头组装新方法

后一篇

极速安装软件的升级版 conda