Linux文件排序和FASTA文件操做

文件排序

seq: 產生一系列的數字; man seq查看其具體使用。咱們這使用seq產生下游分析所用到的輸入文件。

# 產生從1到10的數，步長爲1
$ seq 1 10
1
2
3
4
5
6
7
8
9
10
# 產生從1到10的數，步長爲1，用空格分割
$ seq -s ' ' 1 10
1 2 3 4 5 6 7 8 9 10
# 產生從1到10的數，步長爲2
# 若是有3個數，中間的數爲步長，最後一個始終爲最大值
$ seq -s ' ' 1 2 10
1 3 5 7 9
$ cat <(seq 0 3 17) <(seq 3 6 18) >test
$ cat test 
0
3
6
9
12
15
3
9
15

sort: 排序，默認按字符編碼排序。若是想按數字大小排序，需添加-n參數。

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

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

$ sort -nu test
0
3
6
9
12
15

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

$ sort -n test | uniq -d
3
9
15

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    3
b    2
c    1
d    1
e    1

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

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

FASTA序列提取

生成單行序列FASTA文件，提取特定基因的序列，最簡單的是使用grep命令。主要用途是匹配文件中的字符串，以此爲基礎，進行一系列的操做。若是會使用正則表達式，將會很是強大。正則表達式版本不少，幾乎每種語言都有本身的規則。

# 生成單行序列FASTA文件
$ cat <<END >test.fasta
> >SOX2
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> >POU5F1
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> >NANOG
> CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
> END
$ cat test.fasta 
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
>POU5F1
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
>NANOG
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT

# grep匹配含有SOX2的行
# -A 1 表示輸出的行中，包含匹配行的下一行 (A: after)
$ grep -A 1 'SOX2' test.fasta 
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC

# 先判斷當前行是否是 > 開頭，若是是，表示是序列名字行，替換掉大於號，取出名字。
# 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
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC

多行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    
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGAC
>POU5F1    
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
CGGAAGGTAGTCGTCAGTGCAGCGAGTCC
>NANOG    
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGG
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGACTGT
 
#使用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' 
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC
>POU5F1
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGTCGGAAGGTAGTCGTCAGTGCAGCGAGTCC
>NANOG
ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT

或者簡單點，直接用前面的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
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC