用Sed進行流編輯
sed命令從文本或者標準輸入中每次讀入一行數(shù)據(jù)。
我們先從簡單的實例出發(fā)，看下該命令怎么將一列中的chrm12,chrom2等轉(zhuǎn)換成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甚至更大的數(shù)據(jù)文件中，而不用打開整個文件進行處理。并且，可以借助重導向?qū)崿F(xiàn)對數(shù)據(jù)處理結(jié)果的輸出。
sed替換命令采用的格式是
s/pattern/replacement/
sed會自動搜索符合pattern的字符串，然后修改為replacement（我們想要修改后的樣子）。一般默認sed只替換第一個匹配的pattern，我們可以通過添加全局標識g將其應用到數(shù)據(jù)的所有行中。
s/pattern/replacement/g
如果我們想要忽略匹配的大小寫，使用i標識
s/pattern/replacement/i
默認sed命令支持基本的POSIX正則表達式（BRE），可以通過-E選項進行拓展（ERE）。很多的Linux命令都這種方式，像常用的grep命令。
再看一個實例，如果我們想把chr1:28647389-28659480這樣格式的文字轉(zhuǎn)換為三列，可以使用：
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上。初看雜亂無章，但是從最大的結(jié)構(gòu)看依舊是
s/pattern/replacement/
先看pattern部分，這是由幾個簡單正則表達式組成的復合體，幾個()括起來的字符串可以單獨看。第一個匹配chr加上一個非冒號的字符，第二個和第三個都是匹配多個數(shù)字。最開始的^表示以chr起始（前面沒有字符），各個括號中間的是對應的字符。整體的pattern的目的就是為了找到文本中符合這種模式的字符串，如果只是想把這個模式找出來的話，幾個括號可以不用加。顯然這幾個括號的作用就是將它們劃分成多個域，幫助sed進行處理?？梢钥吹?code>replacement部分存在\1,\2,\3，它恰好對應()的順序。這樣我們在中間插入\t制表符，就可以完成我們想要的功能：將原字符串轉(zhuǎn)換為三列。
我本身對字符串并不是非常熟悉，懂一些元字符，可能講解的不是很到位。不熟悉正則表達式的朋友，可以學習和參考下學習正則表達式，是我從Github上Copy到的非常好的學習資料，有興趣也可以Fork學習。
上山的路總是有很多條，我們下面看下其他實現(xiàn)該功能的辦法：
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這樣格式的文字轉(zhuǎn)換為三列，可以使用：
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上。初看雜亂無章，但是從最大的結(jié)構(gòu)看依舊是
s/pattern/replacement/
先看pattern部分，這是由幾個簡單正則表達式組成的復合體，幾個()括起來的字符串可以單獨看。第一個匹配chr加上一個非冒號的字符，第二個和第三個都是匹配多個數(shù)字。最開始的^表示以chr起始（前面沒有字符），各個括號中間的是對應的字符。整體的pattern的目的就是為了找到文本中符合這種模式的字符串，如果只是想把這個模式找出來的話，幾個括號可以不用加。顯然這幾個括號的作用就是將它們劃分成多個域，幫助sed進行處理?？梢钥吹?code>replacement部分存在\1,\2,\3，它恰好對應()的順序。這樣我們在中間插入\t制表符，就可以完成我們想要的功能：將原字符串轉(zhuǎn)換為三列。
我本身對字符串并不是非常熟悉，懂一些元字符，可能講解的不是很到位。不熟悉正則表達式的朋友，可以學習和參考下學習正則表達式，是我從Github上Copy到的非常好的學習資料，有興趣也可以Fork學習。
上山的路總是有很多條，我們下面看下其他實現(xiàn)該功能的辦法：
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會輸出每一行的結(jié)果，用replacement替換pattern，但實際中我們可能會因此得到不想要的結(jié)果。比如下面的這個例子。
如果我們想要抓出gtf文件第九列的轉(zhuǎn)錄名，可能會使用以下命令
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
我們可以發(fā)現(xiàn)一些沒有轉(zhuǎn)錄名行的結(jié)果是輸出整行，這可不是我們想要的。一種解決辦法是在使用sed之前先抓出有transcript_id的行。其實sed命令本身也可以通過選項和參數(shù)設定解決這個問題，這里我們可以用-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
注意方括號內(nèi)^是非（取反）的意思。
解釋如下：
- 首先，匹配字符串"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的功能特性遠遠不止這些，有待于大家更多地挖掘。不過需要注意的是，盡量讓工具干它最擅長的事情。如果是復雜地大規(guī)模計算，還是最好寫個Python腳本。
KISS原則:
Keep Incredible Sed Simple
高級Shell用法
子shell
首先需要記住連續(xù)命令和管道命令的區(qū)別：前者是簡單地一個一個按順序運行程序（一般用&&或者;）；后者前一個程序的輸出結(jié)果會直接傳到下一個命令程序的輸入中（這不就是流程化操作么，用|分隔）。
子shell可以讓我們在一個獨立的shell進程中執(zhí)行連續(xù)命令。
首先看個例子
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
發(fā)現(xiàn)僅僅加了個括號，結(jié)果就不同了。第二個命令就用了子shell，它把兩個echo命令放進單獨的空間執(zhí)行后將結(jié)果傳給下游。
子shell在對gtf文件進行操作時有個非常有意思有用的用處。我們?nèi)绻雽?code>gtf文件排序，但是又想要保留文件頭部注釋信息，我們就能夠用兩次grep操作分別抓出注釋和非注釋信息，然后又把它結(jié)合在一起。下面看看效果，用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確實能夠給我們提供非常有用的操作去組合命令實現(xiàn)想要的功能。
命令管道及進程替換
很多生信命令行工具需要提供多個輸入和輸出參數(shù)，這用在管道命令里可能會導致非常低效的情形（管道只接受一個標準輸入和輸出）。幸好，我們可以使用命令管道來解決此類問題。
命名管道，也成為FIFO（先入先出，額，這不是隊列么:smile:）。它是一個特殊的排序文件，命名管道有點像文件，它可以永久保留在你的文件系統(tǒng)上（估計本質(zhì)就是文件吧~）。
我們用mkfifo來生成它
wangsx@SC-201708020022:~/tmp$ ls -l fqin
prw-rw-rw- 1 wangsx wangsx 0 9月 3 20:45 fqin
可以它看它權(quán)限的第一個字符是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就可以上游輸出數(shù)據(jù)到processing_tool的命名管道；同理out1.fq out2.fq可以是命名管道用來寫進輸出數(shù)據(jù)。
但這樣我們每次都得不停地創(chuàng)建和刪除這些文件，解決辦法是使用匿名管道，也叫進程替換。
不能光說，看看例子就知道和理解了。
wangsx@SC-201708020022:~/tmp$ cat <(echo "hello, process substitution")
hello, process substitution
echo命令運行后使用了進程替換，產(chǎn)生匿名文件，然后匿名文件被重導向cat命令。
把它用到工具上，就變成了(假定上游zcat下游執(zhí)行grep命令)
processing_tool --in1 < (zcat file1) --in2 < (zcat file2) --out1 (gzip > outfile1) --out2 (gzip > outfile2)
關于Linux數(shù)據(jù)處理工具內(nèi)容全部整理發(fā)布在我的博客上。詳情點擊
