用Sed進行流編輯
sed命令從文本或者標準輸入中每次讀入一行數(shù)據(jù)舶吗。
我們先從簡單的實例出發(fā)征冷，看下該命令怎么將一列中的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甚至更大的數(shù)據(jù)文件中肴捉，而不用打開整個文件進行處理。并且叔收，可以借助重導向實現(xiàn)對數(shù)據(jù)處理結果的輸出齿穗。
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這樣格式的文字轉換為三列洛二，可以使用：
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加上一個非冒號的字符垒迂，第二個和第三個都是匹配多個數(shù)字械姻。最開始的^表示以chr起始（前面沒有字符），各個括號中間的是對應的字符机断。整體的pattern的目的就是為了找到文本中符合這種模式的字符串楷拳，如果只是想把這個模式找出來的話，幾個括號可以不用加吏奸。顯然這幾個括號的作用就是將它們劃分成多個域欢揖，幫助sed進行處理》芪担可以看到replacement部分存在\1,\2,\3她混，它恰好對應()的順序。這樣我們在中間插入\t制表符泊碑，就可以完成我們想要的功能：將原字符串轉換為三列坤按。
我本身對字符串并不是非常熟悉，懂一些元字符馒过，可能講解的不是很到位臭脓。不熟悉正則表達式的朋友，可以學習和參考下學習正則表達式腹忽，是我從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這樣格式的文字轉換為三列，可以使用：
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加上一個非冒號的字符，第二個和第三個都是匹配多個數(shù)字冗恨。最開始的^表示以chr起始（前面沒有字符）答憔，各個括號中間的是對應的字符。整體的pattern的目的就是為了找到文本中符合這種模式的字符串掀抹，如果只是想把這個模式找出來的話虐拓，幾個括號可以不用加。顯然這幾個括號的作用就是將它們劃分成多個域傲武，幫助sed進行處理蓉驹。可以看到replacement部分存在\1,\2,\3谱轨，它恰好對應()的順序戒幔。這樣我們在中間插入\t制表符吠谢，就可以完成我們想要的功能：將原字符串轉換為三列土童。
我本身對字符串并不是非常熟悉，懂一些元字符工坊，可能講解的不是很到位献汗。不熟悉正則表達式的朋友敢订，可以學習和參考下學習正則表達式，是我從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會輸出每一行的結果虹统，用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
我們可以發(fā)現(xiàn)一些沒有轉錄名行的結果是輸出整行，這可不是我們想要的夸赫。一種解決辦法是在使用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
注意方括號內^是非（取反）的意思。
解釋如下：
- 首先切平，匹配字符串"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ū)別：前者是簡單地一個一個按順序運行程序（一般用&&或者;）；后者前一個程序的輸出結果會直接傳到下一個命令程序的輸入中（這不就是流程化操作么福扬，用|分隔）腕铸。
子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)僅僅加了個括號铛碑，結果就不同了恬惯。第二個命令就用了子shell，它把兩個echo命令放進單獨的空間執(zhí)行后將結果傳給下游亚茬。
子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確實能夠給我們提供非常有用的操作去組合命令實現(xiàn)想要的功能。
命令管道及進程替換
很多生信命令行工具需要提供多個輸入和輸出參數(shù)颂砸，這用在管道命令里可能會導致非常低效的情形（管道只接受一個標準輸入和輸出）噪奄。幸好，我們可以使用命令管道來解決此類問題人乓。
命名管道勤篮，也成為FIFO（先入先出，額色罚，這不是隊列么:smile:）碰缔。它是一個特殊的排序文件，命名管道有點像文件戳护，它可以永久保留在你的文件系統(tǒng)上（估計本質就是文件吧~）金抡。
我們用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就可以上游輸出數(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ù)處理工具內容全部整理發(fā)布在我的博客上愕把。詳情點擊
