劉小澤寫于19.3.9
終于搞完了預答辯，想想昨晚熬夜到四點做完了PPT，真是不容易啊，推送繼續(xù)來，未來的時間會越來越充裕，我也期待繼續(xù)和大家一起學習??

今天繼續(xù)挑一些我個人認為比較重要的單行命令，記錄下來

循環(huán)少不了

尤其是在處理多條染色體的時候，比如人類1-22條常染色體+X、Y、M染色體

有幾個注意的問題：

# 只處理數(shù)字的話，例如只處理常染色體
for i in {1..22};do echo $i;done
# 如果再加上其他非數(shù)字，例如這樣
for i in {1..22,X,Y,M};do echo $i;done
# 得到的結果是不對的，并不是輸出1-22以及XYM，而是輸出
1..22
X
Y
M
# 如果想同時得到數(shù)字和非數(shù)字輸出，要把非數(shù)字內(nèi)容放大括號外邊，并且空格分隔
for i in {1..22} X Y M;do echo $i;done
# 如果想逆序輸出，就在大括號中先寫大數(shù)，后寫小數(shù)，例如
for i in {5..2} X Y M;do echo $i;done
5
4
3
2
X
Y
M

用cut提取列

我們都知道的用法是：

# 提取1-5以及10-15列(注意中間的逗號)
cut -f1-5,10-15 data
# 還可以指定分隔符(默認tab分隔)
cut -d' ' -f1-5 data

另外，如果想從第5列輸出到最后一列，其實我們不需要知道一共有多少列，直接這樣：

cut -f5- data
# 注意5后面的-就是表示一直到最后一列

這樣普通的用法cut還是很方便的，但是相比awk，cut還是略顯簡陋，因為不支持自定義列，例如：

# 先提取第5列，然后提取第2列(注意如果操作gtf文件，需要看看有沒有頭注釋信息，有的話需要先去掉)
cat gencode.v28.gtf | grep -v "#" | awk '{print $5"\t"$2}' | head
# 但是用cut就不能實現(xiàn)
cat gencode.v28.gtf | grep -v "#" | cut -f5,2 | head
# 得到的結果還是按照先第2列后第5列排序

多學點awk

假設現(xiàn)在有一個test.txt文件，其中包含了1-22條染色體的基因信息。其中第一列是基因名，第2列是染色體編號。

我們現(xiàn)在想按照test.txt中不同染色體的信息區(qū)分開，1號染色體基因信息存放到chr1.txt中，2號染色體的基因信息存放到chr2.txt中...

可能一開始我們想這么做

for i in {1..22};do awk '$2 == $i' test.txt >chr$i.txt;done

但是要知道awk和shell是不融合的，awk中并不知道在shell中定義的$i

因此我們指定尋找$2 == $i，awk不清楚要找什么，但事實上，我們就是想讓awk去根據(jù)在shell中設定的變量去引用，只需要加一個參數(shù)-v：

for i in {1..22};do awk -v i=$i '$2 == i' test.txt >chr$i.txt;done

-v的作用就是在awk中引入變量variable

例如還有：

# 自定義過濾指標(例如要求第四列值大于10)
for i in {1..22}; do awk -v i=$i -v threshold=10 '$2 == i && $4 > threshold' test.txt > chr$i.txt; done

今天時間有限，先發(fā)這些吧

歡迎關注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!