Linux学习-文件排序和FASTA文件操作

环境变量的补充

PATH只是众多环境变量中的一个变量,用于存储可执行文件所在的目录,以便在用户输入命令时可以查询的到。尤其是自己写的脚本或安装的程序,系统不会知道它们在哪个路径下,需要我们去提供给系统这些新的路径,学名叫设置环境变量。

此外常用到的环境变量还有LD_LIBARY_PATH: 指定动态链接库 (so文件)的位置,一般在安装软件出错时会用到;PYTHONPATH: 指定Python的安装包的路径;PERL5LIB: 指定perl的安装包的路径。

设置环境变量要注意2点:1. 设置新的环境变量时一般要包含原始的环境变量,不能覆盖;2. 注意自己的目录和系统环境变量的目录的顺序,想让哪个先被找到,就先放哪个。

文件排序

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 3 6 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
  1 12
  2 15

# 换一个文件看的更清楚
ct@ehbio:~$ cat <<END >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$nBEGIN表示在文件读取前先设置基本参数;与之相对应的是END,只文件读取完成之后进行操作。不以BEGIN, END开头的{}就是文件读取、处理的部分。

# 管道符还记得吧,后台回复 管道 可查看
# awk的操作就是镀金上一步的结果,去除多余的空白,然后调换2列
ct@ehbio:~$ sort test2 | uniq -c | awk 'BEGIN{OFS="\t";}{print $2, $1}'
a    3
b    2
c    1
d    1
e    1

对两列文件,安照第二列进行排序, sort -k2,2n

# 第二列按数值大小排序
ct@ehbio:~$ 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行的顺序与上一步结果的差异
ct@ehbio:~$ 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命令。

grep在前面也提到过,以后还会经常提到,主要用途是匹配文件中的字符串,以此为基础,进行一系列的操作。如果会使用正则表达式,将会非常强大。正则表达式版本很多,几乎每种语言都有自己的规则,本文档不会展开,用到哪个提哪个。

# 生成单行序列FASTA文件
ct@ehbio:~$ cat <<END >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"; print seq["SOX2"]}' test.fasta
>SOX2
ACGAGGGACGCATCGGACGACTGCAGGACTGTC

多行FASTA序列提取要麻烦些,一个办法就是转成单行序列,用上面的方式处理。

sedtr都为最常用的字符替换工具。

ct@ehbio:~$ cat <<END >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' ' '
>SOX2     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGAC >POU5F1     CGGAAGGTAGTCGTCAGTGCAGCGAGTCCGT CGGAAGGTAGTCGTCAGTGCAGCGAGTCC >NANOG     ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGG ACGAGGGACGCATCGGACGACTGCAGGACTGTC ACGAGGGACGCATCGGACGACTGCAGGACTGT 

# 把最后一个空格替换为换行符
ct@ehbio:~$ 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 隔开
ct@ehbio:~$ 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

# 把所有的空格替换掉
ct@ehbio:~$ 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键转换为换行符
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

原文发布于微信公众号 - 生信宝典(Bio_data)

原文发表时间:2017-06-16

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏潇涧技术专栏

Builtin Lint Detectors (1)

本文主要介绍的是Lint工具中自带的与Android开发相关的lint检查项,通过查看lint检查项的描述及其代码实现,我发现这里面存在不少应用开发编码的Bes...

731
来自专栏Python研发

Django之Model世界

django为使用一种新的方式,即:关系对象映射(Object Relational Mapping,简称ORM)

1242
来自专栏大内老A

ASP.NET MVC以ModelValidator为核心的Model验证体系: ModelValidatorProviders

前面篇文章我们分别介绍用真正用于实施Model验证的ModelValidator(《ASP.NET MVC以ModelValidator为核心的Model验证体...

1805
来自专栏开发技术

Cassandra-java操作——基本操作

  接着上篇博客,我们来谈谈java操作cassandra; 上篇博客的环境:jdk1.7 + python2.7.10 + cassandra2.2.8; 由...

1132
来自专栏逍遥剑客的游戏开发

VS2005中Nebula3数据类型的调试信息显示

1717
来自专栏hotqin888的专栏

用golang递归构建无限级树状目录json数据和数据库

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/hotqin888/article/det...

3032
来自专栏技术博客

Entity Framework 系统约定配置

Code First之所以能够让开发人员以一种更加高效、灵活的方式进行数据操作有一个重要的原因在于它的约定配置。现在软件开发越来越复杂,大家都试图将软件设计的越...

962
来自专栏SDNLAB

Open vSwitch系列之openflow版本兼容

众所周知Open vSwitch支持的openflow版本从1.0到1.5版本(当前Open vSwitch版本是2.3.2)通过阅读代码,处理openflow...

56413
来自专栏逆向技术

PE格式第六讲,导出表

                PE格式第六讲,导出表 请注意,下方字数比较多,其实结构挺简单,但是你如果把博客内容弄明白了,对你受益匪浅,千万不要看到字...

1936
来自专栏生信宝典

Bookdown文档生成教程

bookdown是一款及其方便的编写技术文档或教材的工具,语法简洁,数据处理灵活。支持Rmarkdown或普通markdown通过pandoc软件转换为HTML...

5705

扫码关注云+社区

领取腾讯云代金券