http://mp.weixin.qq.com/s/Rw6ivF_qBPCWiiGL7lxc_g
文件内容操作(二)
文件排序
seq:产生一系列的数字;man seq查看其具体使用。我们这使用seq产生下游分析所用到的输入文件。
#产生从1到10的数，步长为1
ct@ehbio:~$ seq 1 10
1
2
3
4
5
6
7
8
9
10
#产生从1到10的数，步长为1，用空格分割
ct@ehbio:~$ seq -s ' ' 1 10
1 2 3 4 5 6 7 8 9 10
#产生从1到10的数，步长为2
#如果有3个数，中间的数为步长，最后一个始终为最大值
ct@ehbio:~$ seq -s ' ' 1 2 10
1 3 5 7 9
#还记得前面提到的标准输入和标准输出吧
#后台回复 标准输入 查看
ct@ehbio:~$ cat <(seq 0 3 17) <(seq 36 18) >test
ct@ehbio:~$ cat test
0
3
6
9
12
15
3
9
15
sort:排序，默认按字符编码排序。如果想按数字大小排序，需添加-n参数。
#可能不符合预期的排序，系统首先排0，然后排1, 3, 6, 9
ct@ehbio:~$ sort test
0
12
15
15
3
3
6
9
9
#按数字大小排序
ct@ehbio:~$ sort -n test
0
3
3
6
9
9
12
15
15
sort -u:去除重复的行，等同于sort | uniq。
ct@ehbio:~$ sort -nu test
0
3
6
9
12
15
sort file | uniq -d:获得重复的行。(d=duplication)
ct@ehbio:~$ sort -n test | uniq -d
3
9
15
sort file | uniq -c:获得每行重复的次数。
#第一列为每行出现的次数，第二列为原始的行
ct@ehbio:~$ sort -n test | uniq -c
1 0
2 3
1 6
2 9
112
215
#换一个文件看的更清楚
ct@ehbio:~$ cat <test2
> a
> b
> c
> b
> a
> e
> d
> a
> END
#第一列为每行出现的次数，第二列为原始的行
ct@ehbio:~$ sort test2 | uniq -c
3 a
2 b
1 c
1 d
1 e
#在执行uniq操作前，文件要先排序，不然结果很诡异
ct@ehbio:~$ 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列
ct@ehbio:~$ sort test2 | uniq -c | awk'BEGIN{OFS="\t";}{print $2, $1}'
a3
b2
c1
d1
e1
对两列文件，安照第二列进行排序, sort -k2,2n。
#第二列按数值大小排序
ct@ehbio:~$ sort test2 | uniq -c | awk'BEGIN{OFS="\t";}{print $2, $1}' | sort -k2, 2n
c1
d1
e1
b2
a3
#第二列按数值大小排序
#第二列相同的再按第一列的字母顺序的逆序排序(-r)
#注意看前3行的顺序与上一步结果的差异
ct@ehbio:~$ sort test2 | uniq -c | awk'BEGIN{OFS="\t";}{print $2,$1}' | sort -k2,2n -k1,1r
e1
d1
c1
b2
a3
FASTA序列提取
生成单行序列FASTA文件，提取特定基因的序列，最简单的是使用grep命令。grep在前面也提到过，以后还会经常提到，主要用途是匹配文件中的字符串，以此为基础，进行一系列的操作。如果会使用正则表达式，将会非常强大。正则表达式版本很多，几乎每种语言都有自己的规则，本文档不会展开，用到哪个提哪个。
#生成单行序列FASTA文件
ct@ehbio:~$ cat <test.fasta
> >SOX2
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> >POU5F1
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> >NANOG
> CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
> END
ct@ehbio:~$ cat test.fasta
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
>POU5F1
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
>NANOG
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
# grep匹配含有SOX2的行
# -A 1表示输出的行中，包含匹配行的下一行(A: after)
ct@ehbio:~$ grep -A 1 'SOX2' test.fasta
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
#也可以使用AWK
#先判断当前行是不是>开头，如果是，表示是序列名字行，替换掉大于号，取出名字。
# sub替换, sub(被替换的部分，要替换成的，待替换字符串)
#如果不以大于号开头，则为序列行，存储起来。
# seq[name]:相当于建一个字典，name为key，序列为值。然后就可以使用name调取序列。
ct@ehbio:~$ awk'BEGIN{OFS=FS="\t"}{if($0~/>/) {name=$0; sub(">","", name);} else seq[name]=$0;}END{print ">SOX2"; printseq["SOX2"]}' test.fasta
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
多行FASTA序列提取要麻烦些，一个办法就是转成单行序列，用上面的方式处理。sed和tr都为最常用的字符替换工具。
ct@ehbio:~$ cat <test.fasta
> >SOX2
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> ACGAGGGACGCATCGGACGACTGCAGGAC
> >POU5F1
> CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
> CGGAAGGTAGTCGTCAGTGCAGCGAGTCC
> >NANOG
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> ACGAGGGACGCATCGGACGACTGCAGG
> ACGAGGGACGCATCGGACGACTGCAGGACTGTC
> ACGAGGGACGCATCGGACGACTGCAGGACTGT
> END
#给>号开头的行的行尾加个TAB键，以便隔开名字和序列
# TAB键不可见，直接看看不大
# \(\)表示记录匹配的内容，\1则表示()中记录的匹配的内容
#后面我们专门讲sed
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGAC
>POU5F1
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT
CGGAAGGTAGTCGTCAGTGCAGCGAGTCC
>NANOG
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGG
ACGAGGGACGCATCGGACGACTGCAGGACTGTC
ACGAGGGACGCATCGGACGACTGCAGGACTGT
#使用cat -A可以显示文件中所有的符号
# ^I表示tab键
# $表示行尾
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta | cat -A
>SOX2^I$
ACGAGGGACGCATCGGACGACTGCAGGACTGTC$
ACGAGGGACGCATCGGACGACTGCAGGACTGTC$
ACGAGGGACGCATCGGACGACTGCAGGAC$
>POU5F1^I$
CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT$
CGGAAGGTAGTCGTCAGTGCAGCGAGTCC$
>NANOG^I$
ACGAGGGACGCATCGGACGACTGCAGGACTGTC$
ACGAGGGACGCATCGGACGACTGCAGG$
ACGAGGGACGCATCGGACGACTGCAGGACTGTC$
ACGAGGGACGCATCGGACGACTGCAGGACTGT$
#把所有的换行符替换为空格
# tr这个命令，前面提到过，若想不起来`man tr`查看
#主意第二个参数，引号内为空格。
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta | tr '\n' ' '
>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC >POU5F1CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC >NANOGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT
#把最后一个空格替换为换行符
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta | tr '\n' ' ' | sed -e 's/ $/\n/'
>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC >POU5F1CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC >NANOGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT
#把' >'替换为换行符 注意被替换的是 空格+大于号
#当连用多个替换命令时，使用-e隔开
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta | tr '\n' ' ' | sed -e 's/ $/\n/' -e 's/ >/\n>/g'
>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC
>POU5F1CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC
>NANOGACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGT
#把所有的空格替换掉
ct@ehbio:~$ sed 's/^\(>.*\)/\1\t/'test.fasta | tr '\n' ' ' | sed -e 's/ $/\n/' -e 's/ >/\n>/g' -e 's/ //g'
>SOX2ACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGAC
>POU5F1CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGTCGGAAGGTAGTCGTCAGTGCAGCGAGTCC
>NANOGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACGAGGGACGCATCGGACGACTGCAGGACTGTCACGAGGGACGCATCGGACGACTGCAGGACTGT
#把TAB键转换为换行符
ct@ehbio:~$ 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
ct@ehbio:~$ 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
