之前写 awk 教程的时候，曾经提到过一些对文本中行列进行某些计算统计的需求，例如使用数组分类求和。一些基本需求 awk 都可以实现，但是写起来稍显复杂。在 R 中使用 dplyr 或者基础函数 两个使用 awk 的例子如下： 例如利用 在 gnu 官网 中，用了下面一句话来介绍 datamash。 GNU datamash is a command-line program which performs basic numeric, textual and statistical operations on input textual data files. datamash 作为一个命令行程序可以对文本文件进行数字和文本相关的基本统计操作（虽说基本但是所有的操作都足够常用高频）。 使用前首先进行安装，可以通过 datamash 的基本调用格式如下： 上面的内容转换为描述语言就是：在 不同版本支持的参数不同，以下参数适用于 v1.4 版本。在分组统计相关的参数中，p/s 前缀分别代表 population 或者 sample。一般而言，sample 对应的计算等同于 R 中对应函数，例如 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: 关于 sskew 值对应区间的说明 使用软件自带的测试数据集 统计每门课的选课人数 查看每门课程的最高分和最低分，在这里使用 因为这是有列名的输入文件，使用 如果想要对文本中的数字进行简化，可以使用如下几个命令： 在实际的服务器管理中，可能需要统计一些系统用户信息。这些信息存放在 假设有一个行为基因列为组织的表达矩阵如下所示 我想要计算每个基因在所有组织中的表达情况（均值标准差），可以进行如下操作 但是上面的格式看起来不是很美观，可以利用 datamash 自身的命令稍加调整 genepred 和 gff3 文件格式信息一致，但是更加便于我们进行操作。genepred 文件有两种格式，其中一种是标准格式，一种是 extended 格式，具体解释如下。 首先检查一下我们的数据情况，可以看出是 extended 格式 统计每个基因的转录本个数，同时输出这些转录本的 id 看看每个染色体有多少基因 如果结合 除此以外，针对 bed 文件也有非常高频的使用场景这里就不再展示了。 说到使用单行命令解决问题，通常都是 perl 的骄傲，当然 awk 也绝对是一把好手。如果特别喜欢用 R，其实命令行里也可以通过 awk 计算 mean min 和 max 都没有压力，使用 END 可以解决问题。 如果要进行分组操作，awk 就需要数组和 首先构造一个两个的矩阵： awk 根据第一列的 key 进行分组统计个数以及求和，可以看出 awk 需要写的各种括号实在是不少。 尽管可能 90% 以上的人说到单行命令都不会想到 R，但 R 确实可以强行单行。不过要单行操作用到的更多是 R 基础函数，很多喜闻乐见的包用起来就不方便了。 当统计基本的描述统计量时，以下为两者的对比，此时感觉 R 的 如果要分组计算，R 可以使用 总之，在 Linux 下有这么好用的一个工具，命令简单且可以轻松通过管道符和 awk 以及 sed 等连用，完成平时分析数据时的很多基本需求。推荐学习和使用。 参考网站： 本文作者：思考问题的熊 版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。 如果你对这篇文章感兴趣，欢迎通过邮箱或者微信订阅我的 「熊言熊语」会员通讯，我将第一时间与你分享肿瘤生物医药领域最新行业研究进展和我的所思所学所想，点此链接即可进行免费订阅。引子
aggregate() 可以方便的进行分组操作，如果能在 linux 中使用更加简洁的单行命令针对数值和字符进行一些基本运算就省去了在 R 终端操作的时间。这篇文章介绍一个 linux 中能满足这类需求的工具 GNU datamash。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 是什么
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 [option]… op1 column1 [op2 column2 …]
option的参数下，对column1列进行op1操作，对column2列进行op2操作。如果使用--group 参数，所有的operation都将会分组进行；如果没有使用--group，所有的operation会对所有的值进行。需要说明的是这里的column1 可以是表示第几列的数字，当使用-H或者--header-in时可以是所选字段的名称，可以使用列名，当 operation 要求输入成对数据的时候使用: 连接，比如spearson 5:6。主要 operation
sstdev 和 sd()是一致的。
groupby：分组等同于 --group或者-g参数。后面指定用于分组的列crosstab：类似于 excel 中的数据透视表 Pivot Table，可以按照两列来处理矩阵，默认是计算在 a,b 中出现的次数transpose：交换行列，转置，等价于 R 中的t()reverse：反转字段顺序，交换列check：验证数据结构，保证每行字段相同
rmdup：删除有重复键值的行
round：四舍五入整数值floor：不小于输入的最小整数值，向上取ceil：不大于输入的最大整数值，向下取trunc：取整数部分frac：取小数部分
sum：求和min：最小值max：最大值absmin：最小绝对值absmax：最大的绝对值range：最大值-最小值
count ：分组中元素个数first：分组中第一个值last：分组中最后一个值rand：分组中一个随机值unique：逗号分割的所有唯一值collapse：逗号分科的所有值countunique：唯一值的个数
mean：均值mode：众数median：中位数q1：第一四分位数q3：第三四分位数iqr：四分位距perc：百分位值，默认取 95，perc:25 等同于q1antimode：最少出现的数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 -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
$ 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 文件操作
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}"
)
$ datamash check < test.genepred
133744 lines, 15 fields
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
其它 one-liner
Rscript -e来强行单行。以下简单举一些例子进行比较。awk 和 datamash
$ 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
getline 来搞定。有时候需要写很长一行命令，并不是很划算。printf "%s\t%d\n" a 1 b 2 a 3 b 4 a 3 a 6 > test.txt
# 统计个数
$ 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
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
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
· 分享链接 https://kaopubear.top/blog/2019-06-19-datamashmanual/