查看原文
其他

linux shell tricks for bioinformatics系列文章之三: 内附部分sRNAs分析pipeline

2017-02-09 天地本无心 生信媛


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



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存