这里简要地整理下Linux用来处理数据文本的工具。具体命令详情请在Linux命令大全中搜索或者查阅其他相关资料。
head
,tail
查看文档头尾。-n
选项可以指定行数。
less
用来查阅文档,q
退出,space bar
翻页,g
第一行,G
最后一行,j
下,k
上,/<pattern>
往下搜索模式,?<pattern>
往上搜索模式,n
前一个匹配字符,N
后一个匹配字符。
less
可以用于debug,查看中间输出结果。比如
step1 input.txt | step2 | step3 > output.txt
# step1,2,3为程序或命令名
可以写为
step1 input.txt | less
step1 input.txt | step2 | less
step1 input.txt | step2 | step3 | less
纯文本信息汇总
wc
命令默认依次输出单词数、行数、总字符数。查看行数使用wc -l
。
如果存在空行,空行会被计数。可以使用grep
命令实现非空行计数
grep -c "[^ \\n\\t]" some_data.bed
ls -lh
以易读形式查看文件大小。
输出文件列数:
# -F指定分隔符,此处假定是table键分隔,默认空格键
awk -F "\t" '{print NF; exit}' some_data.bed
怎么去除注释的元数据行呢?怎么计数非注释行行数呢?
可以使用tail
结合awk
,试试gtf(基因组注释文件)
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ head -n 6 Homo_sapiens.GRCh37.75.gtf
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1 pseudogene gene 11869 14412 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
可以看到注释行是5行,我们利用tail
试试去掉它
# 注意此处 -n后接的"+"号
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ tail -n +5 Homo_sapiens.GRCh37.75.gtf | head -n 1
#!genebuild-last-updated 2013-09
发现还有一行没去掉
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ tail -n +6 Homo_sapiens.GRCh37.75.gtf | head -n 1
1 pseudogene gene 11869 14412 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
成功搞定,然后结合前面提到的awk
命令即可计算行数。
上面方法鲁棒性不够(人为的确定行数),一种更为通用的方法是grep
结合awk
命令
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 1
1 pseudogene gene 11869 14412 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
推荐使用这种。
cut
可以处理列数据,-f
选项指定列,可以是一个范围(比如2-8),注意不能用它给列排序。
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 10 | cut -f 3
gene
transcript
exon
exon
exon
transcript
exon
exon
exon
transcript
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 10 | cut -f 3-5
gene 11869 14412
transcript 11869 14409
exon 11869 12227
exon 12613 12721
exon 13221 14409
transcript 11872 14412
exon 11872 12227
exon 12613 12721
exon 13225 14412
transcript 11874 14409
-d
选项可以指定分隔符,比如-d,
指定,
为分隔符。
使用column
命令来格式化输出,上次的命令结果输出明显没对齐,我们把它对齐看看:
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 10 | cut -f 3-5 | column -t
gene 11869 14412
transcript 11869 14409
exon 11869 12227
exon 12613 12721
exon 13221 14409
transcript 11872 14412
exon 11872 12227
exon 12613 12721
exon 13225 14412
transcript 11874 14409
注意,使用这个命令是为了好观察,不要把用它处理然后把结果传入文本(会导致程序处理文件效率降低,因为文本解析速度会下降)。
cut
和column
默认以\t
为分隔符,这里也能够用-s
选项指定。
先把之前的tab分隔文件弄成逗号分隔文件,然后使用-s
选项看看:
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 10 | cut -f 3-5 | awk '{FS="\t";OFS=",";}{print $1,$2,$3}'
gene,11869,14412
transcript,11869,14409
exon,11869,12227
exon,12613,12721
exon,13221,14409
transcript,11872,14412
exon,11872,12227
exon,12613,12721
exon,13225,14412
transcript,11874,14409
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -v "^#" Homo_sapiens.GRCh37.75.gtf | head -n 10 | cut -f 3-5 | awk '{FS="\t";OFS=",";}{print $1,$2,$3}'| column -s "," -t
gene 11869 14412
transcript 11869 14409
exon 11869 12227
exon 12613 12721
exon 13221 14409
transcript 11872 14412
exon 11872 12227
exon 12613 12721
exon 13225 14412
transcript 11874 14409
grep
处理速度非常之快,能用它尽量用它。--color=auto
可以激活颜色(标记匹配文字),更方便查看。
-v
选项排除匹配到的,-w
进行完全匹配。这样可以防止,你想排除abc
结果把abc1
,abcd
也排除掉了。
-B
指定输出包括匹配到的前多少行,比如-B1
就是前一行;-A
指定输出包括匹配到的后多少行,比如-A2
就是包括了后两行。-C
指定输出包括匹配到的前后多少行。
grep
支持基本正则表达式,-E
指定支持扩展表达式,或者用egrep
命令。
-c
选项对匹配的行计数;-o
选项只抽离输出匹配的部分
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -E -o 'gene_id "\w+"' Homo_sapiens.GRCh37.75.gtf | head -n 5
gene_id "ENSG00000223972"
gene_id "ENSG00000223972"
gene_id "ENSG00000223972"
gene_id "ENSG00000223972"
gene_id "ENSG00000223972"
发现冗余项非常多,如果我们只要唯一的呢,怎么办?
wsx@wsx-ubuntu:~/Work/research/Promoter_Research$ grep -E -o 'gene_id "(\w+)"' Homo_sapiens.GRCh37.75.gtf | cut -f2 -d" "| sed 's/"//g' | sort | uniq | head -n 10
ENSG00000000003
ENSG00000000005
ENSG00000000419
ENSG00000000457
ENSG00000000460
ENSG00000000938
ENSG00000000971
ENSG00000001036
ENSG00000001084
ENSG00000001167
虽然我的笔记本呼啦啦作响,但是还是非常快就跑完了。
后面内容还不少,分两次讲吧~