全文6,271字,閱讀16分鐘。

awk 是處理文本文件的一個應用程序,幾乎所有的Linux以及MacOS都自帶這個程序。

在這篇文章中,我想給大家介紹如何用這個程序來解決一些基本的生物信息數據處理和文本處理的問題,特別適合對此不熟悉的同學和讀者朋友。

簡述

我們日常項目中很多的數據分析和處理工作並非都需要編寫複雜的程序才能完成,很多小修小改、查找、替換、簡單的數據計算等工作,其實可以用一些Linux/MacOS中自帶的命令行工具來完成。awk 就是這一類工具中的一個,它依次處理文件中的每一行,並讀取裡面的每一個欄位,對於我們在生信中很多每行格式都相同的文本文件來說,awk 可能是最方便的一個工具,不但可以省去很多不必要的腳本和程序,還可以通過對它的靈活應用,快速理解手中的數據。

其實,把 awk 說成是一個程序工具並不十分準確。實際上,它還是一種解釋型的編程語言(類似於Perl),即寫即用,響應快,錯了重改也方便,也有人習慣稱這一類編程語言為腳本語言。不過在這裡我只介紹它的命令行用法,對於很多生物信息的數據分析場景,應該是足夠的,與之類似的還有 sed。

事先說明一下,awk 畢竟是命令行工具,所以我在這篇文章中所用到的例子都只能在Linux或者MacOS命令行中才能執行。如果你用的是windows電腦,那麼需要安裝一個命令行工具纔可以(比如git bash——segmentfault.com/a/1190,而不能是win自帶的命令行)。

基本用法

接下來就是正文了,awk 其實十分簡單,它在命令行中的基本用法就是下面這個形式:

# 標準用法的形式
$ awk 處理動作 文件名
# 例子
$ awk {print $0} demo.vcf

在這個例子中,demo.vcf 是 awk 要處理的文本文件——注意我這裡反覆強調必須是文本文件,而不是BAM或者.gz這一類非文本文件,如果想用 awk 處理這類文件,那麼需要先轉換為文本文件纔行,假如文件不大,那麼可以不做單獨轉換,直接用管道操作來完成即可。回到剛剛的例子,demo.vcf 前面的單引號內有一個大括弧(注意,這個單引號是必須的,而在包含判斷、輸出等複雜語句的時候大括弧也是必須的),裡面是對文件中每一行內容的處理動作,比如這裡是:print $0,其中 print 是列印命令,而 $0 代表當前完整的一行,所以上面這個命令的執行結果就是把 demo.vcf 每一行都原樣列印出來

由於文章這裡不便打開 demo.vcf 文件,我們就用標準輸入來演示這個例子吧。

$ echo "this is a variant in vcf file" | awk {print $0}
this is a variant in vcf file

上面代碼中,echo 也是 linux 一個命令行工具,作用是輸出標準輸入的文本內容,這裡不展開,print $0 就是把標準輸入的 "this is a variant in vcf file" 這一句話,重新輸出到屏幕。

大家應該也注意到上面的命令裏有一個 「|」 豎線,這個就是 Linux/Unix/MacOS 的管道操作符,一個非常有用的符號。它可以把前一個命令的結果作為標準輸入傳輸到後一個命令中去處理,而且還可以多重串聯下去,就像成語接龍一樣,前一個管道處理完再傳給下一個管道去處理,然後再下一個,如果你願意的話,甚至可以一直接下去,這樣做的好處是減少系統 IO,加快處理速度。我前面說到 awk 只能處理文本文件,那當我們的文件不是文本格式時,比如是 gz 壓縮文件或者BAM文件的時候,要用 awk 處理的話,就需要先做轉換然後通過管道把數據傳過給 awk 來分析,比如:

$ samtools view demo.bam | awk {print $0}

這裡就是先通過 samtools view 將 demo.bam 轉為可讀的文本,然後用管道("|")把數據傳到後面的 awk 中去處理——我們這裡為例子的方便只是做原樣輸出。不過,通過這種形式進行數據分析的時候,應該注意的地方是,被處理的 demo.bam 文件不能太大,否則,管道前一個命令(samtools view)轉換出來的文本信息會一直累積到計算機內存中,最後很可能把機器內存撐爆。

默認情況下,awk 將根據空格和製表符(tab),把每一行自動切分成若干個欄位,並在系統裏依次用 $1,$2,$3,... 代表第一個欄位、第二個欄位、第三個欄位等等。

$ echo "this is a variant in vcf file" | awk {print $4}
variant

上述代碼中,$4 代表了 "this is a variant in vcf file" 這一句話中的第四個欄位 "variant"。如果把這一段話換為一份文件,那麼這個命令就會把那份文件中各行的第四列都列印輸出出來。

除此之外,對於某些不是以空格和tab作為分隔符存儲的文件,或者在文件中的某一列的信息中是以其它分隔符串接起來的,比如 VCF 的 INFO 那一列,它是 VCF 的第八列,該列中的信息往往比較豐富,並且各個欄位之間是通過逗號(,)連接在一起的,如下:

CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525

此時,如果我們想獲得該 INFO 中的某一個欄位信息,那麼這個時候除了要提取出這一列之外,還需要通過自定義輸入分隔符才能將其進行切割。自定義輸入分隔符,在 awk 中用的是 -F 參數,例子:

$ awk {if($1!~"^#"){print $8}} demo.vcf | awk -F , {print $2}
CMDB_AC=420
CMDB_AC=441
CMDB_AC=441
CMDB_AC=842
CMDB_AC=842
CMDB_AC=534

這個命令完成了 demo.vcf 中 INFO 這一列信息中第2個欄位信息的提取。其中 通過 -F 參數重新設置了輸入分隔符為逗號,從而完成了對INFO的切分,然後再提取出欄位。該操作命令中前半部分的語句 "if($1!~"^#"){print $8}" 是把VCF 的header信息過濾掉,由於 VCF 的 Header 中每一行都是以 # 開頭的,所以 $1!~"^#" 就可以忽略掉這些 # 開頭的行。

BEGIN 語句

另外在上面的例子中,除了使用 -F 參數之外,還有另一個方法也可以完成這個操作,就是通過 BEGIN 語句,在執行實際命令之前初始化輸入分隔符:

$ awk {if($1!~"^#"){print $8}} demo.vcf | awk BEGIN{FS=","}{print $2}

上面這個命令也能夠實現同樣的功能,只不過,在 BEGIN 中就不是通過 -F 參數了,代替之的是一個內置變數 FS。同時,如果需要的話,我們還可以在其中設置多重分隔符,如 FS="[:,]"(或者 -F [:,]),代表同時用冒號和逗號作為輸入分隔符切分數據,這種方式在比較複雜的文本環境中應用起來會更加方便。

此外,既然可以設置輸入分隔符,自然也可以定義輸出分隔符。我這裡還是用 BEGIN 定義作為例子:

$ awk {if($1!~"^#"){print $8}} demo.vcf | awk BEGIN{FS=","; OFS="###"}{print $1,$2,$3}
CMDB_AF=0.030044###CMDB_AC=420###CMDB_AN=13442
CMDB_AF=0.031047###CMDB_AC=441###CMDB_AN=13553
CMDB_AF=0.031047###CMDB_AC=441###CMDB_AN=13553
CMDB_AF=0.050419###CMDB_AC=842###CMDB_AN=16135
CMDB_AF=0.050419###CMDB_AC=842###CMDB_AN=16135
CMDB_AF=0.053564###CMDB_AC=534###CMDB_AN=9525

awk 默認的輸出分隔符是空格,這個例子在 BEGIN 語句中則通過 OFS 參數將輸出分隔符修改為 "###",當然,最後想用什麼輸出分隔符,完全取決於我們的實際需要。

有BEGIN就有END

與 BEGIN 語句對應的是 END 語句。awk 在默認情況,是每處理完一行數據,就可以輸出一次。但是有時候,我們不希望那麼快地把信息輸出,特別是當我們進行累加統計的時候,比如,計算基因組的平均覆蓋深度,這個時候,我們希望能夠在讀取完整份文件之後,再統一輸出結果。這個時候 END 就派上用場了,顧名思義,只有直到文件結束了,它才把最後結果輸出出來。

$ awk {total_depth += ($3-$2)*$4; total_length += ($3-$2);}END{print total_depth/total_length} seq_depth.bed
53.4

上面的代碼裡面 seq_depth.bed 是通過 bedtools2 生成的一份基因組覆蓋深度文件,為 bed 格式,第一列是染色體ID,第二列是起始位置,第三列是終止位置,第四列是該區域各個位點的覆蓋深度,其中每一個bed區域裏各個位點的深度都是一樣的,所以只留下一個值,這也是為什麼我在上面累加深度的時候需要用 ($3-$2)*$4 的原因。在整個命令中,直到最後讀完整份 seq_depth.bed 才print 出最終的平均深度,比如這裡的 53.4。

內置變數

其實,除了上述通過 $+數字 的形式表示某個欄位之外,awk 本身還有一些默認變數。其中包括,變數 NF 表示當前行按照輸入分隔符切分之後一共有多少列(或者說多少欄位),所以 $NF就表示最後一個欄位,在一些列數非常多的文件中 NF 是很有用的,我們不用數數 數到眼花,也能立刻獲得最後一個欄位,或者立刻知道每一行都有多少欄位。如:

$ echo "this is a variant in vcf file" | awk {print $NF}
file

$(NF-1)則代表倒數第二個欄位,$(NF-2)表示倒數第三,依此類推。

有表示列數,自然也就有表示行數的。awk 中的變數 NR 就是表示當前所處理的是第幾行。

$ awk {if($1!~"^#"){print $8}} demo.vcf | awk -F , {print NR")", $2}
1) CMDB_AC=420
2) CMDB_AC=441
3) CMDB_AC=441
4) CMDB_AC=842
5) CMDB_AC=842
6) CMDB_AC=534

在這個例子中唯一需要注意的是,print 輸出的欄位中,如果各個欄位之間沒通過逗號隔開,那麼輸出時,中間也不會加入任何分隔符,比如這裡 NR 後面直接跟了 ")",輸出的時候 ")" 就緊貼著行數出來。

awk 內置的變數還有這些,其實有不少我們在上面已經用過了,這裡再做匯總:

FILENAME:當前文件名
FS:欄位分隔符,默認是空格和製表符
RS:行分隔符,用於分割每一行,默認是換行符
OFS:輸出欄位的分隔符,用於列印時分隔欄位,默認為空格
ORS:輸出記錄的分隔符,用於列印時分隔記錄,默認為換行符
OFMT:數字輸出的格式,默認為%.6g

內置函數

awk 除了有好用的內置變數之外,也提供了不少好用的內置函數。這些函數可以讓我們很方便地對原始數據進行一些基本的處理。比如,tolower() 用於將字元轉換為小寫。

$ awk {if($1!~"^#"){print $8}} demo.vcf | awk -F , tolower($2)}
cmdb_ac=420
cmdb_ac=441
cmdb_ac=441
cmdb_ac=842
cmdb_ac=842
cmdb_ac=534

上述代碼中,tolower把對應欄位都轉換為了小寫。其他的常用函數還有如下這些:

tolower():字元轉為小寫。
length():返回字元串長度。
substr():返回子字元串。
sin():正弦。
cos():餘弦。
sqrt():平方根。
rand():隨機數。

但是真正完整的內置函數列表,需要查看awk手冊

條件判斷

awk還可以自定條件判斷語句,只把符合條件預期的結果輸出。命令模式如:

$ awk 條件 動作 文件名

需要注意的是,條件判斷要寫在動作之前。請看下面一個例子:

$ awk $6 > 40 demo.vcf
這裡只把 demo.vcf 中第六列大於40(也就是質量值>40)的行輸出出來。

我們也可以寫一個正則表達式,把符合匹配條件的行輸出,比如上述例子也出現過,把VCF的Header過濾掉:

$ awk $1!~/^#/ demmo.vcf

條件判斷是很自由的,我們可以依據自己的需要任意設計條件,包括大於、小於、等於、匹配、與或非、異或等等邏輯判斷條件都可以設置。

if 語句

除了上面提到的條件判斷之外,awk也有 if 語句,可以用來編寫更加靈活複雜的條件,上述例子也已經有所涉及。

$ awk {if($1!~"^#" && $6>40){print $8}} demo.vcf
CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525

這個代碼的意思就是,只把VCF文件中,質量值大於40的變異的INFO信息 print 出來。

有 if 結構自然也會有與它配對的 else 部分。

$ awk {if(/^#/){print $0}else{if($6>40){print $0}}} demo.vcf

小結

有關 awk 的基本用法,這篇文章就介紹到這裡吧,下一篇是 awk 的進階(進階篇已經優先在我的知識星球中給所有星友分享了)。雖然還不是十分全面,但我覺得能夠掌握好上面的使用方法,並靈活應用,其實已經可以用一行命令處理很多基本的分析需求了,不必為了一個小功能費勁去寫一個程序。更重要的是,像 awk 這一類的處理程序,可以很方便且快速地幫我們瞭解一個剛接觸到的數據,對於加深對數據的理解都是很有好處的,並不用總是要等到你學會寫程序之後纔行。更多的 awk 教程推薦這本書《Effective awk programming》,貌似只有英文版。

參考鏈接

ruanyifeng.com/blog/201

runoob.com/linux/linux-

如果喜歡更多的生物信息和組學文章,歡迎搜索並關注我的微信公眾號: helixminer

你還可以讀

  • zhuanlan.zhihu.com/p/59
  • 兩種不同的科研模式
  • GATK4全基因組數據分析最佳實踐 ,我以這篇文章為標誌,終結當前WGS系列數據分析的流程主體問題
  • 人類基因組的Phasing原理是什麼?

推薦閱讀:

相關文章