用Sed进行流编辑
sed命令从文本或者标准输入中每次读入一行数据。
我们先从简单的实例出发，看下该命令怎么将一列中的chrm12,chrom2等转换成chr12，chr2的格式。
wangsx@SC-201708020022:~/tmp$ cat chrms.txt
chrom1 3214482 3216968
chrom1 3214234 3216968
chrom1 3213425 3210653
wangsx@SC-201708020022:~/tmp$ sed 's/chrom/chr/' chrms.txt
chr1 3214482 3216968
chr1 3214234 3216968
chr1 3213425 3210653
虽然示例文件处理仅仅只有三行，但我们可以将这种处理方式运用到上G甚至更大的数据文件中，而不用打开整个文件进行处理。并且，可以借助重导向实现对数据处理结果的输出。
sed替换命令采用的格式是
s/pattern/replacement/
sed会自动搜索符合pattern的字符串，然后修改为replacement（我们想要修改后的样子）。一般默认sed只替换第一个匹配的pattern，我们可以通过添加全局标识g将其应用到数据的所有行中。
s/pattern/replacement/g
如果我们想要忽略匹配的大小写，使用i标识
s/pattern/replacement/i
默认sed命令支持基本的POSIX正则表达式（BRE），可以通过-E选项进行拓展（ERE）。很多的Linux命令都这种方式，像常用的grep命令。
再看一个实例，如果我们想把chr1:28647389-28659480这样格式的文字转换为三列，可以使用：
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | \
> sed -E 's/^(chr[^:]+):([0-9]+)-([0-9]+)/\1\t\2\t\3/'
chr1 28647389 28659480
我们聚焦在第二个命令sed上。初看杂乱无章，但是从最大的结构看依旧是
s/pattern/replacement/
先看pattern部分，这是由几个简单正则表达式组成的复合体，几个()括起来的字符串可以单独看。第一个匹配chr加上一个非冒号的字符，第二个和第三个都是匹配多个数字。最开始的^表示以chr起始（前面没有字符），各个括号中间的是对应的字符。整体的pattern的目的就是为了找到文本中符合这种模式的字符串，如果只是想把这个模式找出来的话，几个括号可以不用加。显然这几个括号的作用就是将它们划分成多个域，帮助sed进行处理。可以看到replacement部分存在\1,\2,\3，它恰好对应()的顺序。这样我们在中间插入\t制表符，就可以完成我们想要的功能：将原字符串转换为三列。
我本身对字符串并不是非常熟悉，懂一些元字符，可能讲解的不是很到位。不熟悉正则表达式的朋友，可以学习和参考下学习正则表达式，是我从Github上Copy到的非常好的学习资料，有兴趣也可以Fork学习。
上山的路总是有很多条，我们下面看下其他实现该功能的办法：
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | sed 's/[:-]/\t/g'
chr1 28647389 28659480
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | sed 's/:/\t/' | sed 's/-/\t/'
chr1 28647389 28659480
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | tr ':-' '\t'
chr1 28647389 28659480
这三种方式看起来都非常简单有效。它处理字符串的思路不是从匹配pattern然后替换入手，不对，应该说是不是从匹配所有pattern然后替换入手。处理的关键是只处理字符串中看似无用的连字符:与-，将其替换成制表符从而轻松完成分割。
sed 's/:/\t/' | sed 's/-/\t/'可以通过-e选项写为sed -e 's/:/\t/' -e 's/-/\t/'，效果等价。
默认sed命令支持基本的POSIX正则表达式（BRE），可以通过-E选项进行拓展（ERE）。很多的Linux命令都这种方式，像常用的grep命令。
再看一个实例，如果我们想把chr1:28647389-28659480这样格式的文字转换为三列，可以使用：
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | \
> sed -E 's/^(chr[^:]+):([0-9]+)-([0-9]+)/\1\t\2\t\3/'
chr1 28647389 28659480
我们聚焦在第二个命令sed上。初看杂乱无章，但是从最大的结构看依旧是
s/pattern/replacement/
先看pattern部分，这是由几个简单正则表达式组成的复合体，几个()括起来的字符串可以单独看。第一个匹配chr加上一个非冒号的字符，第二个和第三个都是匹配多个数字。最开始的^表示以chr起始（前面没有字符），各个括号中间的是对应的字符。整体的pattern的目的就是为了找到文本中符合这种模式的字符串，如果只是想把这个模式找出来的话，几个括号可以不用加。显然这几个括号的作用就是将它们划分成多个域，帮助sed进行处理。可以看到replacement部分存在\1,\2,\3，它恰好对应()的顺序。这样我们在中间插入\t制表符，就可以完成我们想要的功能：将原字符串转换为三列。
我本身对字符串并不是非常熟悉，懂一些元字符，可能讲解的不是很到位。不熟悉正则表达式的朋友，可以学习和参考下学习正则表达式，是我从Github上Copy到的非常好的学习资料，有兴趣也可以Fork学习。
上山的路总是有很多条，我们下面看下其他实现该功能的办法：
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | sed 's/[:-]/\t/g'
chr1 28647389 28659480
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | sed 's/:/\t/' | sed 's/-/\t/'
chr1 28647389 28659480
wangsx@SC-201708020022:~/tmp$ echo "chr1:28647389-28659480" | tr ':-' '\t'
chr1 28647389 28659480
这三种方式看起来都非常简单有效。它处理字符串的思路不是从匹配pattern然后替换入手，不对，应该说是不是从匹配所有pattern然后替换入手。处理的关键是只处理字符串中看似无用的连字符:与-，将其替换成制表符从而轻松完成分割。
sed 's/:/\t/' | sed 's/-/\t/'可以通过-e选项写为sed -e 's/:/\t/' -e 's/-/\t/'，效果等价。
默认，sed会输出每一行的结果，用replacement替换pattern，但实际中我们可能会因此得到不想要的结果。比如下面的这个例子。
如果我们想要抓出gtf文件第九列的转录名，可能会使用以下命令
wangsx@SC-201708020022:~/database$ zgrep -v "^#" gencode.v27lift37.annotation.gtf.gz | head -n 3 | \
> sed -E 's/.*transcript_id "([^"]+)".*/\1/'
chr1 HAVANA gene 11869 14409 . + . gene_id "ENSG00000223972.5_2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2_2"; remap_status "full_contig"; remap_num_mappings 1; remap_target_status "overlap";
ENST00000456328.2_1
ENST00000456328.2_1
我们可以发现一些没有转录名行的结果是输出整行，这可不是我们想要的。一种解决办法是在使用sed之前先抓出有transcript_id的行。其实sed命令本身也可以通过选项和参数设定解决这个问题，这里我们可以用-n选项关闭sed输出所有行，在最末的/后加p只输出匹配项。
wangsx@SC-201708020022:~/database$ zgrep -v "^#" gencode.v27lift37.annotation.gtf.gz | head -n 3 | sed -E -n 's/.*transc
ript_id "([^"]+)".*/\1/p'
ENST00000456328.2_1
ENST00000456328.2_1
注意方括号内^是非（取反）的意思。
解释如下：
- 首先，匹配字符串"transcript_id"之前0或多个任意字符（
.表示除换行键的任意字符）。 - 然后，匹配和捕获一个或多个不是引号的字符，用的是
[^"]
+号的使用是一种非贪婪的方法。很多新手会用*，这是贪婪操作，往往会得不偿失，需要注意喔。
wangsx@SC-201708020022:~/database$ sed -E 's/transcript_id "([^"]+)".*/\1/' greedy_example.txt
ENSMUST00000160944
wangsx@SC-201708020022:~/database$ sed -E 's/transcript_id "(.*)".*/\1/' greedy_example.txt
ENSMUST00000160944"; gene_name "Gm16088
使用*时它会尽量多地去匹配符合要求的模式。
我们也可以用sed命令来获取特定范围的行，比如说我要取出头10行，可以使用
sed -n '1,10p' filename
20到50行
sed -n '20,50p' filename
当然sed的功能特性远远不止这些，有待于大家更多地挖掘。不过需要注意的是，尽量让工具干它最擅长的事情。如果是复杂地大规模计算，还是最好写个Python脚本。
KISS原则:
Keep Incredible Sed Simple
高级Shell用法
子shell
首先需要记住连续命令和管道命令的区别：前者是简单地一个一个按顺序运行程序（一般用&&或者;）；后者前一个程序的输出结果会直接传到下一个命令程序的输入中（这不就是流程化操作么，用|分隔）。
子shell可以让我们在一个独立的shell进程中执行连续命令。
首先看个例子
wangsx@SC-201708020022:~/database$ echo "this command"; echo "that command" | sed 's/command/step/'
this command
that step
wangsx@SC-201708020022:~/database$ (echo "this command"; echo "that command") | sed 's/command/step/'
this step
that step
发现仅仅加了个括号，结果就不同了。第二个命令就用了子shell，它把两个echo命令放进单独的空间执行后将结果传给下游。
子shell在对gtf文件进行操作时有个非常有意思有用的用处。我们如果想对gtf文件排序，但是又想要保留文件头部注释信息，我们就能够用两次grep操作分别抓出注释和非注释信息，然后又把它结合在一起。下面看看效果，用less进行检查：
wangsx@SC-201708020022:~/database$ (zgrep "^#" gencode.v27lift37.annotation.gtf.gz; \
> zgrep -v "^#" gencode.v27lift37.annotation.gtf.gz | \
> sort -k1,1 -k4,4n) | less
可以看到，子shell确实能够给我们提供非常有用的操作去组合命令实现想要的功能。
命令管道及进程替换
很多生信命令行工具需要提供多个输入和输出参数，这用在管道命令里可能会导致非常低效的情形（管道只接受一个标准输入和输出）。幸好，我们可以使用命令管道来解决此类问题。
命名管道，也成为FIFO（先入先出，额，这不是队列么:smile:）。它是一个特殊的排序文件，命名管道有点像文件，它可以永久保留在你的文件系统上（估计本质就是文件吧~）。
我们用mkfifo来生成它
wangsx@SC-201708020022:~/tmp$ ls -l fqin
prw-rw-rw- 1 wangsx wangsx 0 9月 3 20:45 fqin
可以它看它权限的第一个字符是p，指代是pipe。说明是个特殊文件。
我们像文件一样对它进行一些操作
wangsx@SC-201708020022:~/tmp$ cat fqin
hello, named pipes
[1]+ 已完成 echo "hello, named pipes" > fqin
wangsx@SC-201708020022:~/tmp$ rm fqin
比如当使用一个生信命令行工具
processing_tool --in1 in1.fq --in2 in2.fq --out1 out1.fq --out2 out2.fq
in1.fq in2.fq就可以上游输出数据到processing_tool的命名管道；同理out1.fq out2.fq可以是命名管道用来写进输出数据。
但这样我们每次都得不停地创建和删除这些文件，解决办法是使用匿名管道，也叫进程替换。
不能光说，看看例子就知道和理解了。
wangsx@SC-201708020022:~/tmp$ cat <(echo "hello, process substitution")
hello, process substitution
echo命令运行后使用了进程替换，产生匿名文件，然后匿名文件被重导向cat命令。
把它用到工具上，就变成了(假定上游zcat下游执行grep命令)
processing_tool --in1 < (zcat file1) --in2 < (zcat file2) --out1 (gzip > outfile1) --out2 (gzip > outfile2)
关于Linux数据处理工具内容全部整理发布在我的博客上。详情点击
网友评论