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

linux极简统计分析工具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 是什么

在 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)或者将所有行按照一组进行处理。

  • 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:样本偏度

    关于 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 -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) 进行许可。

扫码即刻沟通

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×