linux统计分析命令datamash

文章目录

引子
datamash 是什么
调用格式与参数
- 主要 operation
- Primary operations:
- Line-Filtering operations:
- Per-Line operations:
- Group-by Numeric operations:
- Group-by Textual/Numeric operations:
- **主要参数**
操作示例
文本基本操作
分组统计操作
生信相关文本
- 基因表达矩阵
- GenePred 文件操作
R 和 datamash
猜你喜欢
写在后面

本文转载自“生信菜鸟团”

引子

之前写 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 是什么

https://www.gnu.org/software/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 中对应函数，例如 sstdev 和 sd() 是一致的。

Primary operations:

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

groupby：分组等同于 --group 或者 -g 参数。后面指定用于分组的列

crosstab：类似于 excel 中的数据透视表 Pivot Table，可以按照两列来处理矩阵，默认是计算在 a,b 中出现的次数

transpose：交换行列，转置，等价于 R 中的 t()

reverse：反转字段顺序，交换列

check：验证数据结构，保证每行字段相同

Line-Filtering operations:

rmdup：删除有重复键值的行

Per-Line operations:

round：四舍五入整数值

floor：不小于输入的最小整数值，向上取

ceil：不大于输入的最大整数值，向下取

trunc：取整数部分

frac：取小数部分

Group-by Numeric operations:

sum：求和

min：最小值

max：最大值

absmin：最小绝对值

absmax：最大的绝对值

range：最大值 - 最小值

Group-by Textual/Numeric operations:

count ：分组中元素个数

first：分组中第一个值

last：分组中最后一个值

rand：分组中一个随机值

unique：逗号分割的所有唯一值

collapse：逗号分科的所有值

countunique：唯一值的个数

Group-by Statistical operations:

mean：均值

mode：众数

median：中位数

q1：第一四分位数

q3：第三四分位数

iqr：四分位距

perc：百分位值，默认取 95，perc:25 等同于 q1

antimode：最少出现的数

pstdev：总体标准差

sstdev：样本标准差

pvar：总体方差

svar ：样本方差

mad：scaled 绝对中位差

madraw：原始绝对中位差

sskew：样本偏度

pskew：总体偏度

skurt：样本超值峰度

pkurt：总体超值峰度

jarque：Jarque-Beta test for normality P 值

dpo：D’Agostino-Pearson Omnibus test for normality P 值

scov：样本协方差（要求 pair 数据）

pcov：总体协方差（要求 pair 数据）

spearson：样本 pearson 相关系数（要求 pair 数据）

ppearson：总体 pearson 相关系数（要求 pair 数据）

主要参数

-f, --full: 进行计算前打印所有行

-g, --group：按照列分组

–header-in：输入文件第一行为列名，不进行计算

–header-out：输出第一行为列名

-i：对文本操作是忽略大小写

-s：分组前先排序，等同于 sort

-t：指定输入文件分隔符

–output-delimiter：指定输出文件分隔符

–narm：跳过 NA

-R, --round：输出结果保留几位小数

-W, --whitespace: 使用一个或者多个空格或者 tab 作为分隔符

操作示例

使用软件自带的测试数据集 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 -n5Score   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 transposemean    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.250Median :3.000Mean   :3.1673rd Qu.:3.750Max.   :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 等连用，完成平时分析数据时的很多基本需求。推荐学习和使用。

参考网站：

main page

datamash manual

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA