linux shell tricks for bioinformatics系列文章之二
上一次写了一个《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之间是以逗号分开的。
现在我们很好奇那个GO号拥有的基因数目最多，写脚本吗？NO，shell一行搞定。
awk '{{split($0,arrary,",");split($0,table,"\t");}{print table[1],length(arrary)}}' go2geneid.txt | sort -nr -k 2 | head
结果如下：
所以我们知道拥有基因数目最多的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呗。
为了方便使用电脑的读者，本文也同步发表在知乎专栏。有兴趣的可以移步知乎
知乎链接： https://zhuanlan.zhihu.com/pybio