搜索
查看: 3809|回复: 0

[linux] linux shell tricks for bioinformatics系列文章之二

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-2-1 12:47:37 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:天地本无心
原文链接:http://mp.weixin.qq.com/s/yaDVh-fYzXxXixlLiGUxUA

上一次写了一个《Linux tricks for bioinformatics系列文章之一》,有读者从后台给我留言指出了一些欠缺的地方,非常感谢这些认真的读者。如果读者想熟悉这些命令,建议将命令手动操作一遍,从管道符的左边到右边依次输入,观察输出,细细体会每一个操作的作用。只有灵活运用这些招式,才能打出一套漂亮的组合拳,将你的数据打得服服帖帖,看到那只向你瞪眼的企鹅了吗?等你来挑战呢。

1. 将fastq格式文件转化为fasta序列文件。在上个礼拜我提供了一种解法,但是后来有小伙伴在后台私信我,在使用的过程当中遇到了一些问题。因为fasta文件的质量行,也有可能是以@开头。当碰到这样的情况的时候,下面的方法就会带来问题。在深表歉意的同时,贡献了另外一种解法。
原先的解法:

sed '/^@/!d;s//>/;N' your.fastq> your.fasta

新解法:

awk '{if(NR%4==1||NR%4==2){print $0}}'  test.fq | sed 's/^@/>/g' |less

2. 有时候,比如说我们需要做一个clustalw,我们需要将多行的fasta文件(multiple line fasta)格式文件,转换成单行的fasta文件(single line fasta)序列文件,怎么办呢。 巧用awk + printf一行搞定。

awk '/^>/ { print n $0;}  !/^>/ {printf "%s", $0, n="\n"}  END {print ""}'  test.fa


3. 在做GO注释的时候,我们一定会用到go2geneid.txt这个文件。以拟南芥为例,它的go2geneid.txt文件如下。
go2geneID.txt文件中,GO号和gene之间是tab分开,gene之间是以逗号分开的。
tricks2-1.jpg
现在我们很好奇那个GO号拥有的基因数目最多,写脚本吗?NO,shell一行搞定。

awk '{{split($0,arrary,",");split($0,table,"\t");}{print table[1],length(arrary)}}'  go2geneid.txt | sort -nr -k 2 | head

结果如下:
tricks2-2.png
所以我们知道拥有基因数目最多的GO的的名称为GO:0005575,这厮居然有26769个基因,查了一下含义,是表示cell component,这也就不足为奇了。

4. 如何统计文件的行数。大家心想,还用说吗? `wc -l`不就好了。但是当这个问题交给awk时,awk的解法是。

awk 'END {print NR}' file.txt

这种情况下,虽然awk可以胜任,但是不建议使用awk了 ,因为已经有更简单的解法。简洁胜于复杂。

5. 当我们想要提取文件的某列时。

cut -f 7 file.txt > out.txt

当然交给万能的awk,则会这么写。

awk -F '\t' '{print $7}'  file.txt > out.txt  

这种时候就要读者自己去权衡了,因为自己习惯并欣赏awk,所以个人反而更喜欢后一种写法。

6.  将GFF里面第三列不是chromosome的行给取出来。

awk '{if ($3 != "chromosome") print $0}'  TAIR10.gff

7.  删除文件里所有的空行

sed 's/^$/d' file

8. 有时候我们不想看到header,直接用文件本身的内容去进行下一步操作,我们需要将header去掉。

sed '1d'  file.txt

这个题目当然也可以用awk去实现,用awk的解法是...
都看了两篇系列文章了,自己去试试呗~~

9.  有时候我们可能会需要使用awk为每一行行头添加字符串。


awk '{print "gene_ID " $0}' file.txt

10. 统计fastq文件里面的reads数目。

grep  -c "^@"  test.fastq

本来准备把这种解法也算一个,但是这种方法依然有一个问题,那就是质量行是有可能是以@开头的,这样做就会导致有点点误差,虽然也不会差太大。

因此应该是用下面的方法。

cat test.fastq | echo $((`wc -l`/4))
同理,也可以很轻松能够地写出统计fasta文件函数的shell命令,就是改成除以2呗。


欢迎到微信公众号订阅我们
生信媛
bio_sxy
sxy.png.jpg

回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-12-1 13:02 , Processed in 0.028535 second(s), 38 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.