002-生信Linux文本处理三剑客

慢慢看,憋着急!很有用!

前言:

首先呢,在你的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编辑


文本查看:

查看命令

  1. head, tail 查看文件头尾

    -n 查看指定行 (默认值为10行)

  2. cat 将文件全部打出来

  3. less 逐页显示文本

    -S 规则输出

    -N显示行号

  4. 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

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

编辑器命令:

q:退出

g:第一行 G:最后一行

/ 向下👇搜索 ? 向上👆搜索

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 被排序的部位,包含startstop两部分, 二者,隔开;每次使用都要加上这两部分【默认排序第一列】
  • -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中的管道命令|就是这样一种体现。管道中的数据可以不被写入磁盘,在更高速的内存中进行处理,就像空中绿道,不被道路上拥挤的车流阻拦,你可以在空中自由骑行,肯定效率高很多,心情也更舒畅~

**问题来了:**你能否用一行管道命令从读取文件开始到输出排序好的feature文件呢?


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

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related