Linux檔案排序和FASTA檔案操作

標籤：換行   產生   rto   工具   文本處理工具   agg   log   包含   數字   

檔案排序seq: 產生一系列的數字; man seq查看其具體使用。我們這使用seq產生下遊分析所用到的輸入檔案。

# 產生從1到10的數，步長為1$ seq 1 1012345678910# 產生從1到10的數，步長為1，用空格分割$ seq -s ‘ ‘ 1 101 2 3 4 5 6 7 8 9 10# 產生從1到10的數，步長為2# 如果有3個數，中間的數為步長，最後一個始終為最大值$ seq -s ‘ ‘ 1 2 101 3 5 7 9$ cat <(seq 0 3 17) <(seq 3 6 18) >test$ cat test 036912153915

sort: 排序，預設按字元編碼排序。如果想按數字大小排序，需添加-n參數。

# 可能不符合預期的排序，系統首先排0，然後排1, 3, 6, 9$ sort test012151533699# 按數字大小排序$ sort -n test033699121515

sort -u: 去除重複的行，等同於sort | uniq

$ sort -nu test03691215

sort file | uniq -d: 獲得重複的行(d = duplication)

$ sort -n test | uniq -d3915

sort file | uniq -c: 獲得每行重複的次數。

# 第一列為每行出現的次數，第二列為原始的行$ sort -n test | uniq -c  1 0  2 3  1 6  2 9  1 12  2 15 # 換一個檔案看的更清楚$ cat <<END >test2> a> b> c> b> a> e> d> a> END # 第一列為每行出現的次數，第二列為原始的行$ sort test2 | uniq -c      3 a      2 b      1 c      1 d      1 e # 在執行uniq操作前，檔案要先排序，不然結果很詭異$ cat test2 | uniq -c      1 a      1 b      1 c      1 b      1 a      1 e      1 d      1 a

整理下uniq -c的結果，使得原始行在前，每行的計數在後。 awk是一個強大的文本處理工具，其處理資料模式為按行處理。每次讀入一行，進行操作。OFS: 輸出檔案的資料行分隔符號 (output file column separtor)；FS為輸入檔案的資料行分隔符號 (預設為空白字元)。awk中的列從第1到n列，分別記錄為$1, $2 … $n。BEGIN表示在檔案讀取前先設定基本參數；與之相對應的是END，只檔案讀取完成之後進行操作。不以BEGIN, END開頭的{}就是檔案讀取、處理的部分。

# awk的操作就是鍍金上一步的結果，去除多餘的空白，然後調換2列$ sort test2 | uniq -c | awk ‘BEGIN{OFS="\t";}{print $2, $1}‘a    3b    2c    1d    1e    1

對兩列檔案，安照第二列進行排序, sort -k2,2n。

# 第二列按數值大小排序$ sort test2 | uniq -c | awk ‘BEGIN{OFS="\t";}{print $2, $1}‘ | sort -k2, 2nc    1d    1e    1b    2a    3 # 第二列按數值大小排序# 第二列相同的再按第一列的字母順序的逆序排序 (-r)# 注意看前3行的順序與上一步結果的差異$ sort test2 | uniq -c | awk ‘BEGIN{OFS="\t";}{print $2,$1}‘ | sort -k2,2n -k1,1re    1d    1c    1b    2a    3

FASTA序列提取產生單行序列FASTA檔案，提取特定基因的序列，最簡單的是使用grep命令。主要用途是匹配檔案中的字串，以此為基礎，進行一系列的操作。如果會使用Regex，將會非常強大。Regex版本很多，幾乎每種語言都有自己的規則。

# 產生單行序列FASTA檔案$ cat <<END >test.fasta> >SOX2> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> >POU5F1> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> >NANOG> CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT> END$ cat test.fasta >SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC>POU5F1ACGAGGGACGCATCGGACGACTGCAGGACTGTC>NANOGCGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT# grep匹配含有SOX2的行# -A 1 表示輸出的行中，包含匹配行的下一行 (A: after)$ grep -A 1 ‘SOX2‘ test.fasta >SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC# 先判斷當前行是不是 > 開頭，如果是，表示是序列名字行，替換掉大於符號，取出名字。# sub 替換, sub(被替換的部分，要替換成的，待替換字串)# 如果不以大於符號開頭，則為序列行，儲存起來。# seq[name]: 相當於建一個字典，name為key，序列為值。然後就可以使用name調取序列。$ awk ‘BEGIN{OFS=FS="\t"}{if($0~/>/) {name=$0; sub(">", "", name);} else seq[name]=$0;}END{print ">SOX2"; print seq["SOX2"]}‘ test.fasta>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC

多行FASTA序列提取要麻煩些，一個辦法就是轉成單行序列，用上面的方式處理。 sed和tr都為最常用的字元替換工具。

$ cat <<END >test.fasta> >SOX2> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> ACGAGGGACGCATCGGACGACTGCAGGAC> >POU5F1> CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT> CGGAAGGTAGTCGTCAGTGCAGCGAGTCC> >NANOG> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> ACGAGGGACGCATCGGACGACTGCAGG> ACGAGGGACGCATCGGACGACTGCAGGACTGTC> ACGAGGGACGCATCGGACGACTGCAGGACTGT> END # 給>號開頭的行的行尾加個TAB鍵，以便隔開名字和序列# TAB鍵不可見，直接看看不大# \(\)表示記錄匹配的內容，\1則表示()中記錄的匹配的內容# 後面我們專門講sed$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta >SOX2    ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC>POU5F1    CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGTCGGAAGGTAGTCGTCAGTGCAGCGAGTCC>NANOG    ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT #使用cat -A 可以顯示檔案中所有的符號# ^I 表示tab鍵# $表示行尾$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | cat -A>SOX2^I$ACGAGGGACGCATCGGACGACTGCAGGACTGTC$ACGAGGGACGCATCGGACGACTGCAGGACTGTC$ACGAGGGACGCATCGGACGACTGCAGGAC$>POU5F1^I$CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT$CGGAAGGTAGTCGTCAGTGCAGCGAGTCC$>NANOG^I$ACGAGGGACGCATCGGACGACTGCAGGACTGTC$ACGAGGGACGCATCGGACGACTGCAGG$ACGAGGGACGCATCGGACGACTGCAGGACTGTC$ACGAGGGACGCATCGGACGACTGCAGGACTGT$ # 把所有的分行符號替換為空白格# 主意第二個參數，引號內為空白格$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | tr ‘\n‘ ‘ ‘>SOX2     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGAC >POU5F1     CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC >NANOG     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGT  # 把最後一個空格替換為分行符號$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | tr ‘\n‘ ‘ ‘ | sed -e ‘s/ $/\n/‘>SOX2     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGAC >POU5F1     CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC >NANOG     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGT # 把  ‘ >‘替換為分行符號 注意被替換的是 空格+大於符號# 當連用多個替換命令時，使用-e 隔開$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | tr ‘\n‘ ‘ ‘ | sed -e ‘s/ $/\n/‘ -e ‘s/ >/\n>/g‘>SOX2     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGAC>POU5F1     CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC>NANOG     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGT # 把所有的空格替換掉$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | tr ‘\n‘ ‘ ‘ | sed -e ‘s/ $/\n/‘ -e ‘s/ >/\n>/g‘ -e ‘s/ //g‘>SOX2    ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC>POU5F1    CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGTCGGAAGGTAGTCGTCAGTGCAGCGAGTCC>NANOG    ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT # 把TAB鍵轉換為分行符號$ sed ‘s/^\(>.*\)/\1\t/‘ test.fasta | tr ‘\n‘ ‘ ‘ | sed -e ‘s/ $/\n/‘ -e ‘s/ >/\n>/g‘ -e ‘s/ //g‘ -e ‘s/\t/\n/g‘ >SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC>POU5F1CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGTCGGAAGGTAGTCGTCAGTGCAGCGAGTCC>NANOGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT

或者簡單點，直接用前面的awk略微做下修改。

# 差別只在一點# 對於單行fasta檔案，只需要記錄一行，seq[name]=$0# 對於多好fasta檔案，需要把每一行序列都加到前面的序列上，seq[name]=seq[name]$0$ awk ‘BEGIN{OFS=FS="\t"}{if($0~/>/) {name=$0; sub(">", "", name);} else seq[name]=seq[name]$0;}END{print ">SOX2"; print seq["SOX2"]}‘ test.fasta>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC

Linux檔案排序和FASTA檔案操作