linux shell tricks for bioinformatics系列文章之三： 内附部分sRNAs分析pipeline
1. 如何取出一个字符穿的前50个字符
cat your.fa | grep -v ">" | cut -c 1-50
cat your.fa | grep -v ">" | sed -r 's/(.{50}).*/\1/g
当然还可以用awk，但此处就不写了。
2. 昨天在生信菜鸟群里，有一个人问，她想把一个Paired end序列的5'端前50个base和3'端的后50个base取出来做mapping, 想在bowtie2里面找参数，据我所知是没有这样的参数的。不过用shell也就一行代码的事。
下面是取出5'端前50个base的命令：
awk '{if(NR%4==2||NR%4==0) {print substr($0,1,50)} else {print $0} }' your.fq
需要注意的是由于第二行的序列和第四行的质量值是一一对应的，所以如果序列只取前50个，表示质量值的那一行也只需要取前50个。
3. 同理取出3'端的后50个base， 然后将结果依然保存为fastq的方法是：
awk '{if(NR%4==2||NR%4==0) {print substr($0,length($0)-5 + 1)} }' your.fq
4. 假如说，我们现在有一个sRNAs的数据，拿到之后呢第一步就是去接头。取几行序列出来，用clustwl之类的工具，就能很轻易地找到sRNAs的adapter序列。
clustal sRNAs.fa
5. 去掉接头序列后，就应该去掉含有测序的ambiguous碱基，通常用N表示。也就是说，如果序列中含有N碱基，这条reads就应该被抛弃。
sed 'N;s/\n/\t/' sRNAs.fa | sed -e '/\t.*N/d' | tr "\t" "\n"
稍微地解释一下这段代码，fasta是以两行为单位，一行header，一行sequence, 所以我们先用sed将两行变一行，换行符变为\t分隔符; 然后将序列中含有N的reads给删掉,注意此时header中有N字符是可以容许的，然后末了再用tr将\t分隔符又又换换行符输出。
6. 去掉“N"这种ambiguous碱基之后，我们有必要做长度筛选，因为植物里面的sRNAs一般是18nt-30nt，而去掉adator之后，可能有些reads特别短了，还有些很长，可能是来自于structural RNAs的降解产物，所以应该去除这些过长或者过短的，只保留18-30的reads进行分析。
sed 'N;s/\n/\t/' test.fa | awk 'BEGIN{OFS="\t";} {if (length($2)>=18 && length($2)<=30) {print $1,"\n",$2}}' | tr -d "\t"
同理，先将两行变一行，然后将序列长度大于等于18以及小于等于30的reads留下来。然后先输出header, 加上换行符之后，再输出序列，然后用tr将序列前面的分隔符去掉。
7. 对于去掉接头，去掉“N”碱基，去掉过短或者过长的sRNAs序列之后，我们下一步就该看看sRNAs有那些特征了，比如一个比较重要的特征就是5'端的碱基频率，因为不同的5’端碱基，决定了sRNAs可能与那种AGO蛋白相结合，然后发挥生理作用，因此这个特征也是非常有意义的。
cat test.fa | grep -v ">" |cut -c 1-1 | sort| uniq -c| awk '{print $2,$1}' | tr " " "\t"
比如我们看看08年戚一军发的那片cell文章里与ago1结合的sRNA的特征。这个图是用excel画的简图。
8. 对于去完接头，以及作了一些前处理的sRNAs序列文件，虽然还没有进行去掉结构RNA等步骤。但这时候的你，已经迫不及待地看看sRNAs的长度分布了，怎么办呢，shell一行搞定。
cat sRNAs.fa | awk 'NR%2==0' | awk '{print length($0)}'|sort |uniq -c|sed 's/^[][ ]*//g'|awk '{print $2,$1}'
拿到上述代码的输入文件之后，用excel打开，一步操作就可以出来如下的长度分布图。虽然图片不好看，但是却能清楚地表达数据本身的信息，即sRNAs的长度分布图的峰值是21nt左右。
突然发现excel也能够提供类似于heatmap的样式，比较好玩，因此放在这儿。
9. 递归地删除空目录
find . -depth -type d -empty -exec rmdir -v {} \;
由此引申，可以递归地删除分析过程中产生的中间文件，以节省存储空间。
10. 监控文件的变化，特别是当大文件比对，然后你想看看output到底有没有变化时。
watch ls -al your.sam