生信Linux文本解决三剑客

作者 : 开心源码 本文共10123个字,预计阅读时间需要26分钟 发布时间: 2022-05-14 共176人阅读

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

前言:

首先呢,在你的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:最后一行

/<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项的proteingrep "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.fasgrep -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 ‘{print1}’ *.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.bedchr1 34 36chr2 38 40chr1 25 39chr3 12 18一个染色体长度文件为:length.txtchr1 54362chr2 35613chr3 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   30427671Chr1  gene    3631    5899Chr1  mRNA    3631    5899Chr1  protein 3760    5630Chr1  exon    3631    3913Chr1  five_prime_UTR  3631    3759Chr1  CDS 3760    3913Chr1  exon    3996    4276Chr1  CDS 3996    4276Chr1  exon    4486    4605

接下来对截取的test.txt进行解决

@3 想根据第二列的feature进行排序【注意sort -k的使用!】sort -k 2,2 test.txt | head -n5=>Chr1  CDS 3760    3913Chr1  CDS 3996    4276Chr1  chromosome  1   30427671Chr1  exon    3631    3913Chr1  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    10001736Chr5  CDS 10004720    10004824Chr5  CDS 10004720    10004824Chr5  CDS 10005070    10005255Chr5  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>


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

说明
1. 本站所有资源来源于用户上传和网络,如有侵权请邮件联系站长!
2. 分享目的仅供大家学习和交流,您必须在下载后24小时内删除!
3. 不得使用于非法商业用途,不得违反国家法律。否则后果自负!
4. 本站提供的源码、模板、插件等等其他资源,都不包含技术服务请大家谅解!
5. 如有链接无法下载、失效或广告,请联系管理员处理!
6. 本站资源售价只是摆设,本站源码仅提供给会员学习使用!
7. 如遇到加密压缩包,请使用360解压,如遇到无法解压的请联系管理员
开心源码网 » 生信Linux文本解决三剑客

发表回复