生信Linux文本处理三剑客

2018-07-13 本文已影响362人刘小泽

慢慢看，憋着急！很有用！

前言：

首先呢，在你的Linux系统中新建一个文件，Thanos.txt(紫薯侠赐予你力量)

@敲入这一段，待会儿肯定给你解释清楚啥意思
echo -e "Hello world\nI love bioinformatics\n\nWelcome to bioinfoplanet\nAnd I will be happy here\n\nWelcome come to our Wechat public account\nbioinfoplanet" > Thanos.txt

然后你会发现多了一个紫薯侠文件，cat Thanos.txt 看是不是输出了几行字？真棒👍你已经开始了第一步

解释一下，echo是打印到终端，你可以理解为打印机，我们赋予它内容，他就给我们打印出来，到Thanos.txt这个本子上。

-e ：内容中有特殊字符就需要加（比如上面👆的换行符\n）

" "中就是要打印的内容

> 输出到哪里

利用这个文件，你可以继续向下学习。当然如果你认为内容太少，你可以vi Thanos.txt 进入vi编辑器输入i编辑

文本查看：

查看命令

head, tail 查看文件头尾

-n 查看指定行 (默认值为10行)
cat 将文件全部打出来
less 逐页显示文本

-S 规则输出

-N显示行号
less界面中可以移动光标，搜索关键字

less可以进行debug，查看中间输出结果
例如有一个管道程序sh1 test.txt | sh2 | sh3 > output.txt
【其中sh1-3都是命令】
在运行之前，我们可以看看每一步是否正确，这样使用：
sh1 test.txt | less
sh1 test.txt | sh2 | less
sh1 test.txt | sh2 | sh3 | less

查看行数：wc -l 输出行数（但是空行也会被算入）=> 解决： grep -c "[^ \\n\\t]" file
去除注释信息（元数据）后查看行数
一般基因组注释GTF文件开始都有几行的注释文件，需要去除以后再查看
grep -v '^@' *.gtf | wc -l
这样做的好处就是接下来还能用awk直接处理列或者计算列数
grep -v '^@' *.gtf | awk -F "\t" '{print NF};exit'
查看列数【以tab分割为例】：
awk -F "\t" '{print NF ; exit}' *.bed

编辑器命令：

q：退出

g：第一行 G：最后一行

/<pattern> 向下👇搜索？<pattern> 向上👆搜索

n向后匹配 N向前匹配

文件输出：

echo -n 不换行输出 -e 处理特殊字符

echo -e "nihao\nworld"                  echo -n "nihao\nworld"                      
=> nihao                                => nihao\nworld
=> world

>覆盖原文件 >>添加到原文件底部

正则表达式regular expression：

它的语法结构有两套系统组成，元字符（metacharacters） + 普通字符

那么问题来了：什么是元字符呢？
A：元字符是表达式的结构，骨架；比如，“我爱生信星球”这个句子中的主谓宾都是固定的，也就是它的元字符是固定的；普通字符可以不同，也就是会有各种语言版本的这句话，来表达我对生信星球的爱

元字符主要由以下字符组成：^ $ . [] {} - ? + () | \

表达式	描述	范例
^	行首标记	`^bioinfoplanet` 匹配以bioinfoplanet起始的行（以下简称bip）
$	行尾标记	`bip$`匹配以bip结尾的行
.	任意字符	`b.p`匹配任意代替.的字母，如bap，bbp...但不能代表两个如baap
[]	其中任意一个	`bi[op]`匹配bio或bip
[^]	除了其中任一个	`bio[^pb]`就是不能匹配biop & biob，其他任意
[a-d]	匹配指定范围内任一个	能匹配a-d任意一个字母
{n}	匹配之前n项	`[0-9]{2}`匹配一个两位数`[0-9][0-9]`
{n, }	至少匹配前面n次	`[0-9]{n, }` 匹配至少是两位数的
{n, m}	最少匹配n次，最多m次	`[0-9]{2,4}` 匹配两位数到四位数
？	匹配之前1个或没有	`bi?p` 匹配bip或bp
*	匹配之前多个或没有	`bi*p`匹配bp或bip/ biip/...
+	匹配之前1个或多个	`bi+p`匹配bip或biip/biiip/...
()	匹配括号中的字符串	(搭配？或或+ 使用 )* `bio(info)？`匹配bio或bioinfo
\|	匹配两侧任一个	`bio\|info` 匹配bio或info
\	转义	`bio\ +info` 匹配bio+info，否则按+格式处理bioinfo

()与[]的差别：()为多选，[]为单选, 对比下就知道

cat SRR519926_1.fastq | egrep 'TA(A+)TA' --color=always | head
@ 匹配了多种模式

re_1.png

cat SRR519926_1.fastq | egrep 'TA[A+]TA' --color=always | head
@ 仅匹配了一种模式

re_2.png

(e)grep：号称“Almost～最快文本搜索”

意思是 global search regular expression(RE) and print out the line
egrep 是grep的拓展模式，支持的元字符较多

生而为搜——它会在每一行搜索匹配内容

-c 显示有多少行被匹配到（count）

-v 过滤掉某些格式的行

-w 完全匹配【比如想删除匹配到abc的，不使用-w可能abc1和abcd都要被删除】

-o 只打印匹配到的内容

-i 忽略匹配字段和匹配内容的大小写

-A/B n：输出匹配内容前/后 n行

--color=always/auto: 始终/自动高亮显示搜索字段

搜索特定信息

@ 搜索特定基因信息(在拟南芥gff文件中查找，下载地址见下方🌰)
grep "AT1G01680" *.gff | head -n5

grep_1.png

排除特定信息

@ 还是搜索这个基因，但排除其中feature项的protein
grep "AT1G01680" *.gff | grep -v "protein" | head -n5

grep_2.png

查找某段序列并输出上下文

@ grep -A n显示后面n行； -B n显示前面n行
@ 下载基因组文件wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas

grep -A 2 "CAAATTGAATTAAG" *.fas
=> 得到下面👇的结果
=> -B的用法你可以试一下

grep_3.png

查找特定序列并计算出现了几次

grep -c "CAAATTGAATTAAG" *.fas 
或者： grep "CAAATTGAATTAAG" *.fas | wc -l
@ 如果单纯输出 就用grep -o

精准匹配某个基因

@ 搭配正则表达式
@ 比如要匹配开头为AT1G250，结尾为3的基因名
grep "AT1G250.*3$" *.gff
=> 得到了两条AT1G25053 & AT1G25083

grep_4.png

并非完美～

虽然grep很强大，但是他也并非十全十美，他的一个不足之处就在于，对于存在换行符的字符串，它会搜不到。例如我们想找TAIR10_chr_all.fas 中95-100之内的‘CCACT’ 碱基

先看一下tail -n100 TAIR10_chr_all.fas | head -n5

grep_5.png

用肉眼就能看到'CCACT'，但是如果用grep搜索是没有结果的
tail -n100 TAIR10_chr_all.fas | head -n5 | egrep 'CCACT'

其实稍加改变，用grep还是能实现，就是不那么优雅：
tail -n100 TAIR10_chr_all.fas | head -n5 | tr -d '\n' | egrep 'CCACT'
tr -d（delete）是删除特定字段

所以，有没有什么更快捷的办法呢？
有的～可以使用emboss套件下的dreg, 它是针对核酸；如果是氨基酸序列和翻译后的蛋白序列，使用preg

tail -n 1000 chr22.fa | head -n 5 | dreg -filter -pattern 'TAATA'

sed:

流编辑器～就像水流一样，按行从数据读取--执行命令--显示结果，一气呵成
整个过程在一个位于内存中叫Pattern Space模式空间中进行，因此不会改变原始文件内容

使用方法：
sed [options] 'Adress Command1;Command2;...' file...

[options]包括：
-n: 只显示经过sed处理的行，保存在模式空间的未处理行不显示【常用】
-i：直接修改读取内容，不输出【慎用，他会修改源文件】
-r：使用拓展正则表达式

Adress工作范围：
可以指定为按行号工作或按过滤条件工作

按行号：使用单独的数字表示某一行；两个数字并逗号分隔表示范围；
n,m+ 表示从n向下m行；n~m表示从n开始的每m行

比如我要输出偶数行: sed -n “2~2 p" Thanos.txt
按匹配条件：'/pattern/command' 在包括pattern的行中执行命令

比如要打印包含bio的行：sed -n '/bio/p' Thanos.txt
比如打印包含bio与包含planet之间的所有行：
sed -n '/bio/, /planet/p' Thanos.txt
还可以与行号连用：sed -n '/bio/, +6 p'Thanos.txt 打印包含bio行以及下面6行

Command命令[多个命令用分号分隔]：
p：复制模式空间的内容，一般与-n连用【否则会一次输出两次】
d：删除【按行号或匹配条件】

删除第5行后面所有的行
sed '5,$d'

i：插入

例如，要在文件开头插入一行，三列name、length、feature
sed '1i name\tlength\tfeature' Thanos.txt

a：追加

例如，要在末尾追加内容【 $表示最后一行】 `sed '$ a bioinfoplanet' Thanos.txt`

c：替换【与a用法相似】
n：匹配的行向下移一行进行操作

例如：seq 6 | sed '2{n;d}' 输出的就是除了3以外的数字
再如：seq 6 |sed 'n;d' 输出奇数
seq 6 | sed -n 'n;p' 输出偶数

!：反向执行

sed '/bio/!d' Thanos.txt 效果等于 sed -n '/bio/p' Thanos.txt

=：打印行号

sed '/bio/!d;=' Thanos.txt

s：替换【支持正则】
s/pattern/replace/flags

其中pattern是支持正则的
flags包括：n（替换第n个匹配）；g（全局匹配）；p（与-n搭配，输出修改的行）；i（忽略大小写）；w（保存修改的行到一个新文件）

一种特殊情况，还比较常用，就是如果替换/怎么办？
方法一：使用转义符\/
方法二：用@ | ! ^ 替换
例如：要更改当前目录中的部分内容
pwd | sed 's@/home/tmp/bio@/home/tmp/bioinfo@'

还有一种特殊情况，就是有时候不想替换掉，只是想把pattern的这部分内容与replace的内容一起输出来：
sed 's/[[:upper:]]/ word = &/' Thanos.txt
意思就是：将Thanos.txt中的大写字母放到&的位置，输出格式就是：word = 大写字母。【特殊字符&用来存储pattern中的内容】

实用的sed单行命令：

@ 删除空行
sed '/^$/d' 

@ 每一行下增加一行空行
sed G

@ 每三行增加一行空白行
sed '0~3G'

@ 在匹配的pattern后面加入一行空白行
sed '/pattern/G'

@ cat的功能实现
sed ''

@ head的功能实现
sed '3 q' =》输出前三行

@ tee功能实现
sed -n 'p; w newfile'

@ grep功能实现
sed -n '/pattern/p'

@ grep -v功能实现
sed -n '/pattern/p!'

@ 计算行数
sed -n '$='

@ 多个内容同时替换(例如将1、2、3替换成4)
sed 's/1\|2\|3/4/'

@ 显示包含“haha”、“xixi”、“yeah”的行
sed '/haha/d!; /xixi/d!; /yeah/d!'

awk：号称“Almost～最强文本操作”

“awk w(o/a)rd”—— 让看似繁复无章的文字遇见它就秒变尴尬😅

工作职责：主打行【如果一个文本文档是一张表格，每一行代表一个记录，每一列代表域，awk就是处理记录专用】

工作流程： 先逐行读取并记录，将行信息整合入$0 => 指定分隔符（默认空格），分割成各个列$1,$2,$3...；

再执行awk 'pattern1 {command1}; pattern2 {command2}...'

⚠️：情况一：如果缺失pattern【也就是利用下面几种运算进行的匹配】，直接执行{command}；
情况二：缺失{command}【在{}中输入的命令，如print】，执行pattern

现在不太懂也没关系，通过实例练习你就明白，下面的实例我会描述对应哪种情况

模拟cat打印行

awk '{print $0}' Thanos.txt | head -n3 
@ 【因为没有pattern，所以直接执行打印命令】～属于情况一

选择列打印

@ 利用.gff文件练习 ～属于情况一
awk '{print $1,$3,$4,$5}' *.gff | head -n5
【当然，这里默认是空格分割，如果想tab分割呢？】
awk '{[print $1 "\t" $3 "\t" $4 "\t" $5 "\t"]}' 
@ 一会练习完下面的cut命令，看看有什么发现？

awk_1.png

列重排

@ 这个功能相当于下面cut的升级版 ～ 属于情况一
@ 可以改变列的顺序，并且可以自定义分隔符
awk '{print $4","$5":"$1":"$3}' *.gff | head -n5
@ 随心所欲，只为你改变！

awk_2.png

转换格式

比如将gff/gtf格式转为bed格式
【注意两点：一、转换后的bed格式分隔符为tab；
二、bed与gff的坐标格式不同：0-坐标系统有BAM、BED、BCFv2、PSL；1-坐标系统有SAM、VCF、GFF、wiggle。需要做出调整】
cat test.gtf | awk '{print $1 "\t" $4-1 "\t" $5}' > test.bed

如果单有这些简单的功能，还算不上Almost最强，加上下面这些你再试试？

刘小泽喊你进阶啦！
- 逻辑运算（<, > , <=, >=, ==, !=）
- 数学运算（+，-， *， /， %）
- 关系运算（与&&，或||，非！）
- 正则运算 （实则就是将你想匹配的东西放在/ /里，然后在它前面加～匹配， !~ 不匹配）
- - 举个小例子1：
```
@ 我想匹配在3号染色体上，长度大于1.5k的注释,看前五行    
@ ～情况二，我们只需要pattern就好
awk '$5-$4 > 15000 && $1 ~/Chr3/' *.gff | head -n5
```

awk_3.png

  > 然后自己再试试匹配1号或2号染色体上，长度小于1.5k的注释，默认输出前10行

- 小例子2:

  ```平均值
  @ 我想看一下3号染色体上编码区（CDS）的平均值 ～情况二
  awk 'BEGIN {len=0;line=0}; $1 ~/Chr3/ && $3 ~/CDS/ { len += ( $5 - $4 );line += 1}; END {print "CDS_mean=" len/line}' *.gff
  =>
  CDS_mean=225.019
  【这里 $1 ~/Chr3/ 中的～是为了匹配正则表达式】
  @ 那么你试试看，1号染色体上外显子exon的平均长度吧
  @ 此处只是用来练习命令，并不表示真实长度。因为重叠的区域（比如起始位点一样，终止位点不同）这样会被记录两次，结果是不准确的；如果要精确统计，需要用编程去重以后来实现
  ```

特殊的变量：
- （上面说的） $Num**表示哪一个字块（这里再次提醒：**$ 0代表当前内容）
- NR代表当前所在行号，想要打印3-4行内容：head -4|tail -2，
  就等同于awk 'NR>=3 && NR<=4 {print $0}'`
- NF表示目前的记录被分割的字段的数目，即Number of Field
  
  例如要输出文件的列数【默认空格分割，这里设为tab分割】：
  awk -F “\t” '{print NF; exit}' test.bed
- FS 指定列分隔符
- OFS 指定列输出分隔符
  
  例如将上游文件的默认分隔符\t更换成，并输出3-5列
  cat test.txt | awk '{FS="\t"; OFS=",";}{print $3,$4,$5}'

常用的还有：

cut：提取

-f 提取指定字段（filed）cut -f 1,3,4,5 GFF3_genes.gff | head
也可以cut -f 1-5
-d: 指定分隔符（默认\t）
-c 截取一定范围的字符，例如cut -c1-3就是截取三个字符

例如，将bed文件第一列的染色体编号去掉
awk '{print $1}' *.bed | cut -c4 反之，要再添加上的话： `awk '{print$ 1}' *.bed | cut -c4 | awk '{print "chr"$1}'`

uniq: 去重

与sort连用，-c 在每列旁边显示该行重复出现的次数
-d ：只输出重复的行

column:格式化输出

一般cut后的结果参差不齐，可以用它对齐输出结果，默认\t
指定分隔符：以逗号为例， -t -s ','
【⚠️：他的使用只是为了好看，不要把它column处理的结果交给下游继续处理，这样会让文本解析速度下降。😋看来Linux还是更习惯参差不齐的文本】

sort: 排序

-k 被排序的部位，包含start 和stop两部分, 二者,隔开；每次使用都要加上这两部分【默认排序第一列】
-n 按数字大小排序
sort -k1,1 -k2,2n my.bed 对第一列首字母按字符排序，对第二列按数值
-r 降序（默认升序）【与-k连用表示对某一列逆序】
-t 指定字块分隔符（默认tab）
-c 检查是否按照某种方式排过顺序【echo $? 返回0表示执行成功】
-V 排序时不用ASCII码方式排列，就显示为正常的排序方式
【比如三个染色体编号chr1, chr2, chr11。不加-V排序结果就是chr1-- chr11 -- chr2; 加了-V就是：chr1 -- chr2 -- chr11】
【一般数据处理过程中不用这个操作，原始的格式处理更快】

join：连接

要求两个文件之间必须有共同点，所以使用join前必须先将文件排序

格式：join -1 <file1_field> -2 <file2_field> <file1> <file2>

共同点通过-1 、 -2传递进来，比如说两个文件的某一列有共同点

例如一个bed文件前三列如下：my.bed
chr1 34 36
chr2 38 40
chr1 25 39
chr3 12 18
一个染色体长度文件为：length.txt
chr1 54362
chr2 35613
chr3 46612

join这两个文件之前先sort
sort -k1,1 my.bed > sorted.bed
sort -c -k1,1 length.txt
以双方第一列为共同点拼接：
join -1 1 -2 1 sorted.bed length.txt > with_length.txt
【假如两个文件没有共同点，比如length.txt少了chr2数据，那么join后的文件也缺少chr2数据】

-a 选项指定哪个文件可以不遵循，在没有共同点时可以单列出来

举个🌰【很大的那种！】：

首先输出一下未排序文件，你可以自己下一个小数据(43M)试试, 人类的比较大(1.2G)我就不下了

wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff【这是拟南芥的基因注释文件，当练手很好用】

先练习下之前学的一些处理：
@0 拿到手数据，可以先ll -h *.gff 看一下文件大小
                less -SN *.gff 看一下文本内容

@1 去除多余的#注释行 与 空行， 打印行号
grep -v "#" *.gff | grep -v "^$" | wc -l
=> 590264行（一般这种使用量大的数据是没有多余行的，但是人类的有。这里只是为了演示用法）

@2 截取gff文件的1-5列，将第二列除去, 输出默认的前10行到test.txt中
cut -f 1,3,4,5 *.gff > test.txt 
cat test.txt | head
=>
Chr1  chromosome  1   30427671
Chr1  gene    3631    5899
Chr1  mRNA    3631    5899
Chr1  protein 3760    5630
Chr1  exon    3631    3913
Chr1  five_prime_UTR  3631    3759
Chr1  CDS 3760    3913
Chr1  exon    3996    4276
Chr1  CDS 3996    4276
Chr1  exon    4486    4605

接下来对截取的test.txt进行处理

@3 想根据第二列的feature进行排序【注意sort -k的使用！】
sort -k 2,2 test.txt | head -n5
=>
Chr1  CDS 3760    3913
Chr1  CDS 3996    4276
Chr1  chromosome  1   30427671
Chr1  exon    3631    3913
Chr1  exon    3996    4276

@4 【进阶】先根据第一列Chr数字大小降序排序，再根据第二列排序
你是不是试过了这个》sort -k 1,1nr -k 2,2 *.txt | head -n5
你会发现输出的结果第一列还是Chr1开头，并没有降序，为什么呢？
=》原因就是，-k 1，1还是根据第一个字段的全部排序，还是根据Chr1的‘C’进行匹配，其实我们只想用第一个字段的数字（也就是第一个字段的第四个字符）进行匹配
=》如何实现？其实sort -k参数 内置了这个功能。使用-k 1.4，1.4 就是根据这种特定方式匹配啦
sort -k 1.4,1.4nr -k 2,2 *.txt | head -n5
=> 
Chr5  CDS 10001590    10001736
Chr5  CDS 10004720    10004824
Chr5  CDS 10004720    10004824
Chr5  CDS 10005070    10005255
Chr5  CDS 10005070    10005255
有没有很好用？！

如果只是想统计一下整体的feature（第三列）情况，可以这样：
cut -f 3 *.gff |sort|uniq -c >feature.txt
=>
 197160 CDS
      7 chromosome
 215909 exon
  34621 five_prime_UTR
  28775 gene
  35386 mRNA
   3911 mRNA_TE_gene
    180 miRNA
    480 ncRNA
  35386 protein
    924 pseudogene
   1274 pseudogenic_exon
    926 pseudogenic_transcript
     15 rRNA
     13 snRNA
     71 snoRNA
    689 tRNA
  30634 three_prime_UTR
   3903 transposable_element_gene
想看哪个feature最多？没问题，一步搞定！
sort -k 1r feature.txt
=>
 215909 exon
 197160 CDS
  35386 protein
  35386 mRNA
  34621 five_prime_UTR
  30634 three_prime_UTR
  28775 gene
   3911 mRNA_TE_gene
   3903 transposable_element_gene
   1274 pseudogenic_exon
    926 pseudogenic_transcript
    924 pseudogene
    689 tRNA
    480 ncRNA
    180 miRNA
     71 snoRNA
     15 rRNA
     13 snRNA
      7 chromosome
【小练习1：】GENCODE下载人类基因组GRCh38注释gff3 ,正好练习下数据库使用

然后统计人类的基因组 feature

想偷懒？给你个机会～ wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz

应该得到如下结果：=>

1237914 exon

735618 CDS

203835 transcript

148007 five_prime_UTR

144591 three_prime_UTR

85439 start_codon

77451 stop_codon

58381 gene

119 stop_codon_redefined_as_selenocysteine

【小练习2:】 相信你能够完成上面的内容，输出了结果。

小前言：我们知道，Linux基于Unix开发，肯定要继承unix的精华，那就是利用小程序的整合去完成大任务。Linux中的管道命令|就是这样一种体现。管道中的数据可以不被写入磁盘，在更高速的内存中进行处理，就像空中绿道，不被道路上拥挤的车流阻拦，你可以在空中自由骑行，肯定效率高很多，心情也更舒畅～

问题来了：<u>你能否用一行管道命令从读取文件开始到输出排序好的feature文件呢？</u>

提示：读取|截取|排序|统计|统计后排序