原创:熊+jimmy+SX 生信技能树
主要是
awk/grep/sed
这三驾马车,加上vi
这个神器,最后辅助一些小工具,包括wc,cat,diff,join,paste,cut,uniq
这里简要地整理下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
这样就可以一步步的查看中间结果啦~
用man可以查看less的详细手册,注意两个参数即可,分别是 -S(可以比较规整地查看内容)和 -N(查看内容每行加相应地行数)。
同理还有 more 命令, cat 命令也可以查看文本。
纯文本信息汇总
文件行数
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(基因组注释文件)为例
(zzz) [zzz@gpu01 mm10]$ head -8 Mus_musculus.GRCm38.84.gtf
#!genome-build GRCm38.p4
#!genome-version GRCm38
#!genome-date 2012-01
#!genome-build-accession NCBI:GCA_000001635.6
#!genebuild-last-updated 2016-01
1 havana gene 3073253 3074322 . + . gene_id "ENSMUSG00000102693"; gene_version "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotype "TEC"; havana_gene "OTTMUSG00000049935"; havana_gene_version "1";
1 havana transcript 3073253 3074322 . + . gene_id "ENSMUSG00000102693"; gene_version "1"; transcript_id "ENSMUST00000193812"; transcript_version "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotype "TEC"; havana_gene "OTTMUSG00000049935"; havana_gene_version "1"; transcript_name "4933401J01Rik-001"; transcript_source "havana"; transcript_biotype "TEC"; havana_transcript "OTTMUST00000127109"; havana_transcript_version "1"; tag "basic"; transcript_support_level "NA";
1 havana exon 3073253 3074322 . + . gene_id "ENSMUSG00000102693"; gene_version "1"; transcript_id "ENSMUST00000193812"; transcript_version "1"; exon_number "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotype "TEC"; havana_gene "OTTMUSG00000049935"; havana_gene_version "1"; transcript_name "4933401J01Rik-001"; transcript_source "havana"; transcript_biotype "TEC"; havana_transcript "OTTMUST00000127109"; havana_transcript_version "1"; exon_id "ENSMUSE00001343744"; exon_version "1"; tag "basic"; transcript_support_level "NA";
可以看到注释行是5行,我们利用tail
试试去掉它
# 注意此处 -n后接的"+"号
(zzz) [zzz@gpu01 mm10]$ tail -n +6 Mus_musculus.GRCm38.84.gtf | head -n 1
1 havana gene 3073253 3074322 . + . gene_id "ENSMUSG00000102693"; gene_version "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotype "TEC"; havana_gene "OTTMUSG00000049935"; havana_gene_version "1";
成功搞定,然后结合前面提到的 awk命令即可计算行数。
上面方法robust不够(人为地确定行数),更为通用的方法是 grep结合 awk命令
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 1
1 havana gene 3073253 3074322 . + . gene_id "ENSMUSG00000102693"; gene_version "1"; gene_name "4933401J01Rik"; gene_source "havana"; gene_biotype "TEC"; havana_gene "OTTMUSG00000049935"; havana_gene_version "1";
cut详解
cut
可以处理列数据, -f
选项指定列,可以是一个范围(比如2-8),注意不能用它给列排序
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 10 | cut -f 3
gene
transcript
exon
gene
transcript
exon
gene
transcript
exon
exon
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 10 | cut -f 3-5
gene 3073253 3074322
transcript 3073253 3074322
exon 3073253 3074322
gene 3102016 3102125
transcript 3102016 3102125
exon 3102016 3102125
gene 3205901 3671498
transcript 3205901 3216344
exon 3213609 3216344
exon 3205901 3207317
-d
选项可以指定分隔符,比如 -d
,指定 ,
为分隔符。
使用 column
命令来格式化输出,上次的命令结果输出明显没对齐,我们把它对齐看看:
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 10 | cut -f 3-5 | column -t
gene 3073253 3074322
transcript 3073253 3074322
exon 3073253 3074322
gene 3102016 3102125
transcript 3102016 3102125
exon 3102016 3102125
gene 3205901 3671498
transcript 3205901 3216344
exon 3213609 3216344
exon 3205901 3207317
使用这个命令是为了好观察,不要把用它处理然后把结果传入文本(会导致程序处理文件效率降低,因为文本解析速度会下降)。
cut
和 column
默认以 \t
为分隔符,这里也能够用 -s
选项指定
先把之前的tab分隔文件弄成逗号分隔文件,然后使用 -s选项看看:
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 10 | cut -f 3-5 | awk '{FS="\t";OFS=",";}{print $1,$2,$3}'
gene,3073253,3074322
transcript,3073253,3074322
exon,3073253,3074322
gene,3102016,3102125
transcript,3102016,3102125
exon,3102016,3102125
gene,3205901,3671498
transcript,3205901,3216344
exon,3213609,3216344
exon,3205901,3207317
(zzz) [zzz@gpu01 mm10]$ grep -v "^#" Mus_musculus.GRCm38.84.gtf | head -n 10 | cut -f 3-5 | awk '{FS="\t";OFS=",";}{print $1,$2,$3}'| column -s "," -t
gene 3073253 3074322
transcript 3073253 3074322
exon 3073253 3074322
gene 3102016 3102125
transcript 3102016 3102125
exon 3102016 3102125
gene 3205901 3671498
transcript 3205901 3216344
exon 3213609 3216344
exon 3205901 3207317
简单了解grep
grep
处理速度非常之快,能用它尽量用它。 --color=auto
可以激活颜色(标记匹配文字),更方便查看。
-v
选项排除匹配到的, -w
进行完全匹配。这样可以防止,你想排除 abc结果把 abc1, abcd也排除掉了。
-B
指定输出包括匹配到的前多少行,比如 -B1
就是前一行;
-A
指定输出包括匹配到的后多少行,比如 -A2
就是包括了后两行;
-C
指定输出包括匹配到的前后多少行;
grep
支持基本正则表达式;
-E
指定支持扩展表达式,或者用 egrep
命令;
-c
选项对匹配的行计数;
-o
选项只抽离输出匹配的部分。
(zzz) [zzz@gpu01 mm10]$ grep -E -o 'gene_id "\w+"' Mus_musculus.GRCm38.84.gtf | head -n 5
gene_id "ENSMUSG00000102693"
gene_id "ENSMUSG00000102693"
gene_id "ENSMUSG00000102693"
gene_id "ENSMUSG00000064842"
gene_id "ENSMUSG00000064842"
发现冗余项非常多,如果我们只要唯一的呢,怎么办?
[zzz@gpu01 mm10]$ grep -E -o 'gene_id "\w+"' Mus_musculus.GRCm38.84.gtf| cut -d " " -f 2| sed 's/"//g' | sort | uniq| head -n 10
ENSMUSG00000000001
ENSMUSG00000000003
ENSMUSG00000000028
ENSMUSG00000000031
ENSMUSG00000000037
ENSMUSG00000000049
ENSMUSG00000000056
ENSMUSG00000000058
ENSMUSG00000000078
ENSMUSG00000000085
(zzz) [zzz@gpu01 mm10]$ grep -E -o 'gene_id "\w+"' Mus_musculus.GRCm38.84.gtf| cut -d " " -f 2| sed 's/"//g' | sort -u| head -n 10
ENSMUSG00000000001
ENSMUSG00000000003
ENSMUSG00000000028
ENSMUSG00000000031
ENSMUSG00000000037
ENSMUSG00000000049
ENSMUSG00000000056
ENSMUSG00000000058
ENSMUSG00000000078
ENSMUSG00000000085
#这俩命令等价
file查看文件编码
(zzz) [zzz@gpu01 mm10]$ file GRCm38.84.id2biotype.txt
GRCm38.84.id2biotype.txt: ASCII text
常用的大型数据文件一般存为ASCII
码形式(像几大基因bank的数据文件),而我们自己认为创建的常为UTF-8
,所以有时候认为处理文件需要会碰到把UTF-8
编码的字符插入到ASCII
码文件里去了。遇到这种问题,我们可以用 hexdump -c
命令查看出错的地方
sort对文本排序
先创建一个bed格式文件:
(zzz) [zzz@gpu01 mm10]$ cat > test.bed
chr1 26 39
chr3 32 47
chr1 40 50
chr1 9 28
chr2 35 54
chr1 10 19
^Z
[5]+ Stopped cat > test.bed
(zzz) [zzz@gpu01 mm10]$ sort test.bed
chr1 10 19
chr1 26 39
chr1 40 50
chr1 9 28
chr2 35 54
chr3 32 47
可以明显看到文本按照第一列进行了排序。 默认, sort用空格或tab键作为域(列)分隔符。如果我们用其他形式的分隔符,需要用 -t选项指定。
下面是对 bed文件最通用的排序命令:
(zzz) [zzz@gpu01 mm10]$ sort -k1,1 -k2,2n test.bed
chr1 9 28
chr1 10 19
chr1 26 39
chr1 40 50
chr2 35 54
chr3 32 47
bedtools
软件都会先用这个命令对bedtools
文件排序。
-k
选项指定某列的排序方式。而每次使用 -k选项都要带上指定列的范围(start, end)。如果只指定一列,就为(start,start)了,像上面命令的 -k1,1
就是。也许你会觉得 -k2,2n
很奇怪,这里的n指定程序把第二列当做数值对待。如果不做设定,都是当做字符对待(shell都是这么对待数值数据的)。所以总结其他这一行命令就是对第一列按照字符排序,第二列按照数值排序。
我们可以用 -c
选项检查一个文件是不是已经按照过某种方式排过序了,一般shell
返回0
表示成功执行。
反向排序用-r
选项。如果你只想反转一列,可以把它加在 -k
选项后。
sort -k1,1 -k2,2nr test.bed
chr1 40 50
chr1 26 39
chr1 10 19
chr1 9 28
chr2 35 54
chr3 32 47
现在我给 test.bed加一行:
cat test.bed
chr1 26 39
chr3 32 47
chr1 40 50
chr1 9 28
chr2 35 54
chr1 10 19
chr11 22 56
你会发现有点奇怪
sort -k1,1 -k2,2n test.bed
chr1 9 28
chr1 10 19
chr1 26 39
chr1 40 50
chr11 22 56
chr2 35 54
chr3 32 47
怎么 chr11在 chr2前面?其实 sort
排序的方式有点像查字典。例子中,命令先比较 c,然后比较 h,然后比较 r,接着比较 1,自然 11会在 2前面了。这里可以添加 -V
选项修改。
sort -k1,1V -k2,2n test.bed
chr1 9 28
chr1 10 19
chr1 26 39
chr1 40 50
chr2 35 54
chr3 32 47
chr11 22 56
是不是觉得这样更可观一些?不过通常在处理数据时不做此处理,符合规范的数据可以让后续处理程序效率更高。
用uniq寻找唯一值
cat test.letter
A
A
B
C
B
C
C
C
D
F
D
uniq test.letter
A
B
C
B
C
D
F
D
#只把连续的相同的字符去掉了,`sort`命令可以把相同字符搞到一起
sort test.letter | uniq
A
B
C
D
F
sort test.letter | uniq -c
2 A
2 B
4 C
2 D
1 F
#加 -c选项计数
sort test.letter | uniq -c | sort -rn
4 C
2 D
2 B
2 A
1 F
Join 命令
用来连接文件。 假设现在我们有两个文件:
cat example.bed
chr1 26 39
chr1 32 47
chr3 11 28
chr1 40 49
chr3 16 27
chr1 9 28
chr2 35 53
cat example_length.txt
chr1 53453
chr2 34356
chr3 24356
我想把第二个文件说明染色体长度添加到第一个文件对应染色体的第三列。 我们首先要给文件排序(使用 join前必须做),然后使用 join命令。
sort -k1,1 example.bed > example_sorted.bed
sort -c -k1,1 example_length.txt
join -1 1 -2 1 example_sorted.bed example_length.txt > example_with_length.txt
cat example_with_length.txt
chr1 26 39 53453
chr1 32 47 53453
chr1 40 49 53453
chr1 9 28 53453
chr2 35 53 34356
chr3 11 28 24356
chr3 16 27 24356
命令基本语法是
join -1 -2
既然名字叫join
,就是两者必须有共同之处,通过共同的支点将两者连为一体。 -1和 -2选项后接参数分别指定了这个支点,也就是连接的域(列)。比如例子中,都是两个文件的第一列。
两个文件中,第一列都共有 chr1(2)(3)。如果第二个文件没有 chr3, join之后也没了!!
我们可以通过 -a选项指定哪一个文件可以不遵循配对
join -1 1 -2 1 -a 1 example_sorted.bed example_length_alt.txt
chr1 26 39 53453
chr1 32 47 53453
chr1 40 49 53453
chr1 9 28 53453
chr2 35 53 34356
chr3 11 28
chr3 16 27
awk 可以说是一门语言了
awk是文本处理的一把好手,虽然它不能像 python, R干一些高级复杂的主题工作,但是它具备完整的命令操作和编程体系,手册参考:http://man.linuxde.net/awk。
首先要明白的是,awk
按行处理数据。在shell知识里,如果把一个文档看做一张表。那么一行就是一个记录,一列就是一个域。可以看出, awk就是按记录处理文本的。
其次是awk
的程序结构是
pattern {action}
pattern
可以是表达式或者正则表达式。pattern
有点像 if
语句,当它满足时就会执行相应的动作。
另一个 awk
核心是它用$0
表示所有列,$1
,$2
...等等表示对应的列。我们可以很方便地用它进行操作。
awk '{print $0}' example.bed
chr1 26 39
chr1 32 47
chr3 11 28
chr1 40 49
chr3 16 27
chr1 9 28
chr2 35 53
awk '{print $1}' example.bed
chr1
chr1
chr3
chr1
chr3
chr1
chr2
#print语句就像动作一样输出你操作的结果。
awk '{ print $2 "\t" $3}' example.bed
26 39
32 47
11 28
40 49
16 27
9 28
35 53
awk '{ print $2 $3}' example.bed
2639
3247
1128
4049
1627
928
3553
wsx@wsx-ubuntu:/tmp$ awk '{ print $2 , $3}' example.bed
26 39
32 47
11 28
40 49
16 27
9 28
35 53
了解上述几个语句的不同。
表示染色体名一般用带 chr或者不带 chr标志两种方式。当我们要用到这两种时,肯定要让它们能够对应起来,也就是转换。 awk命令可以非常方便地添加 chr标记。
下面我先把例子文件的 chr去掉,然后加上试试。
awk '{ print $1}' example.bed
chr1
chr1
chr3
chr1
chr3
chr1
chr2
awk '{ print $1}' example.bed | cut -c 4
1
1
3
1
3
1
2
awk '{ print $1}' example.bed | cut -c4 | awk '{print "chr"$1}'
chr1
chr1
chr3
chr1
chr3
chr1
chr2
awk作为一门编程语言,它支持各种操作符(运算,逻辑,判断)喔。
awk '$3 - $2 >18' example.bed
chr1 9 28
awk '$1 ~/chr1/ && $3 - $2 > 10' example.bed
chr1 26 39
chr1 32 47
chr1 9 28
这里 ~ 符号用来匹配正则表达式
还有 awk
存在一些变量,像NR
表示行号, OFS
表示输出分隔符等。
awk 'NR >= 3 && NR <= 5' example.bed
chr3 11 28
chr1 40 49
chr3 16 27
如果我们想把 gtf文件转换成为 bed格式,可以使用
head -n 1000 Homo_sapiens.GRCh37.75.gtf | awk '!/^#/{ print $1 "\t" $4-1 "\t" $5} ' | head -n 3
1 11868 14412
1 11868 14409
1 11868 12227
用sed进行流编辑
sed
工作流:读取$\to$
执行命令$\to$
显示。默认情况,所有的命令都会一个叫做在模式空间(pattern buffer)的缓冲区进行。因此不会改变原始输入文件的内容。
语法规则
sed
即支持在命令行中用单引号输入执行命令,也支持执行含有sed命令的文件。使用方式如下:
sed [‐n] [‐e] 'command(s)' files
sed [‐n] ‐f script files
选项
sed
后面首先需要跟参数,支持的参数有:
-n
:禁止显示所有输入内容,只显示经过sed
处理的行(常用)
-e
:直接在命令列模式上进行 sed
的动作编辑,接要执行的一个或者多个命令
-f
:执行含有sed
动作的文件-r
:sed
的动作支持的扩展正则(默认基础正则)
-i
:直接修改读取的文件内容,不输出。
作用区域
默认情况下,sed
命令会作用于文本数据的所有行。如果只想作用于某些行时,则需要使用在命令通过行号或者文本过滤的方式前指明作用区域。
行号
使用数字行号时,类似于R中的向量子集提取。单独数字表示某一行,逗号分割指示行范围,另外 m,+n
表示从m行开始向下n行, m~n
表示从m行开始的每n行。
例如sed ‐n '2~2 p' test.txt
输出偶数行, sed ‐n '1~2 p'
test.txt输出奇数行内容。
文本过滤
'/pattern/ command'
可以只在包含pattern 的行中执行命令。如 sed ‐n '/hello/ p' test.txt
只会打印出包含hello的行。sed ‐n '/hello/, /world/ p' test.txt
打印两者之间的所有行。
特殊情况下也可以将文本过滤和行号结合使用,sed ‐n '/hello/,+5 p' test.txt
打印第一次出现hello的下面5行
命令
p 复制
复制模式空间中的内容,如果不和 -n参数连用,每一行都会在屏幕输出两次,一行正常输出一行复制,结合 -n 参数后就可以打印需要的内容。
d 删除
支持按照行号或者匹配来删除。
i 插入
有的时候一个结果文件没有essay-header,使用sed 可以轻松完成。在匹配位置之前插入内容。
sed '1i name\tlength\foldchange' test.txt
a 追加
和插入命令的区别在于在匹配位置后一行插入内容,,如果想在末尾插入一行信息时将a auther:zhaofei' test.txt`
c 行替换
有了行的删除插入和追加自然也就会有行替换。
sed '$c auther:zhaofei' test.txt
y 字符转换
sed中的y命令可以实现一一映射的字符替换(注意和s命令的区别)。
[address]y/inchars/outchars/
#inchars中的第一个字符会被转换为outchars中的第一个字符,第二个字符会被转换成outchars中的第二个字符,直持续到处理完指定字符。如果inchars和outchars的长度不同会报错
l 输出隐藏字符
类似与cat -A,但是显示隐藏字符形式不同。
w 写入新文件
增强版的cp,只复制自己想要的东西,也可以将一个文件按不同的筛选条件分开保存。
sed -n -e '/a/ w a.txt' -e '1,10 w b.txt' test.txt
`sed '1,10d;w new.txt' test.txt
r 读取文件
如果现在a文件的某个地方插入b文件,如在第三行插入
sed '3r b.txt' a.txt
e 执行外部命令
首先要区别于参数-e,这个e是在''里面的命令
sed 'e echo "hello"' test.txt
i 反向执行
在命令前加!反向执行命令
sed '/hello/!d' test.txt
#效果和下面的命令一致
sed -n '/hello/p' test.txt
n 匹配行下移一行操作
提前读取当前行的下一行内容,并且覆盖当前模式空间中的行
# 首先匹配到第三行,然后移动到第四行进行删除
# 所及结果是
1
2
3
5
seq 5 |sed 'n;d' #效果类似与输出奇数行
seq 5 |sed -n 'n;p' #效果类似与输出偶数行
= 打印行号
sed '/hello/!d;=' test.txt
s 替换
通用写法 [address1[,address2]]s/pattern/replacement/[flags]
这里pattern部分支持正则匹配,flags包括n(替换第n个匹配项),g(全局替换),p(输出改变的行,结合-n),i(忽略大小写匹配),w(保存改变的行到新文件)。
如果要替换的的内容包括了 /,第一种方式是使用 /进行转义,第二种方法是使用 @ | ! ^作为分隔符。
有时候我们会对文件中的目录进行替换,可以下面的写法
pwd |sed 's@/home/zhaofei@/home/feizhao@'
pwd |sed 's^/home/zhaofei^/home/feizhao^'
pwd |sed 's|/home/zhaofei|/home/feizhao|
pattern支持各种正则表示法,例如
^ 行首
$ 行尾
. 换行符之外的任意单个字符
? 匹配之前项0次或者一次
+ 匹配1次或者多次
* 匹配0次或者多次
{n} 匹配n次
{n,} 匹配至少n次
{m,n} 至少m,最多n
[] 匹配任意一个
[-] 匹配范围中的一个
[^] 排除字符
| 或者
# posix字符类
[:alnum:] 字母数字
[:alpha:] 字母
[:blank:] 空格制表符
[:digit:] 数字
[:lower:] 小写字母
[:upper:] 大写字母
[:punct:] 标点
[:space:] 所有空白符(换行符和回车)
#元字符
\s 单个空白
\w 单词
在进行匹配替换时,我们有时候并不想删除匹配的内容,只是希望其以另一种形式和替换内容一起出现。在sed中,特殊字符 &用来存储匹配模式中的内容。
例如
sed 's/[[:digit:]]/number = &/ test.txt
sed 单行命令 可以完全替代上面所有的命令。
# 删除空行
sed '/^$/d'
#每行后增加空行
sed G
# 在每5行后增加一空白行
gsed '0~5G'
# 在匹配式样“regex”的行之后插入一空行
sed '/regex/G'
# 在匹配式样“regex”的行之前和之后各插入一空行
sed '/regex/{x;p;x;G;}'
# 过滤所有的html标签
sed 's/<[^>]*>//g ; /^$/d' html.txt
# 代替 cat 功能
sed ''
# 代替 head 功能
sed '10 q'
# Dos2unix
sed 's/^M$//'
# Unix2dos
sed 's/$/\r/'
# nl(添加行号)
# sed行号会独占一行
sed = input.file | sed 'N;s/\n/\t/'
# tee
sed ‐n 'p; w new.txt'
# 代替 grep 功能
sed ‐n '/pattern/p'
# grep -v
sed ‐n '/pattern/p!'
# 计算行数
sed -n '$='
# 多个内容同时替换
sed 's/a\|b\|c/d/' tmp.txt
# 将每两行连接成一行
sed '$!N;s/\n/ /'
# 如果当前行以等号开头,将当前行并到上一行末尾
# 并以单个空格代替原来行头的“=”
sed -e :a -e '$!N;s/\n=/ /;ta' -e 'P;D'
# 显示包含“AAA”、“BBB”或“CCC”的行(任意次序)
sed '/AAA/!d; /BBB/!d; /CCC/!d'