搜索
查看: 2165|回复: 0

[linux] linux shell tricks for bioinformatics之三: 内附部分sRNAs分析pipeline

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-2-21 19:33:42 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:天地本无心
原文链接:
https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247483765&idx=1&sn=1cd4544c1abd3f6259934e03b1d069dc&chksm=e9e028d4de97a1c254fb6a3ffc675656d9fc80a7f0c1b00cf5fc6e0fb4876f4cb4c36e806294&scene=38&key=08fb5c62c8a3c99d6edbefd0441d52eddad86ff9bc0dfe1114313e957d84f557f4c12b764d7cb8de18e04d9a9ff79f25bf632274658b39c79dcc8cbdae90bdde5317e379da4d4afd9c7463c6c82e163b&ascene=0&uin=MjA1MzEwNzcwMg%3D%3D&devicetype=android-23&version=26050434&nettype=WIFI&abtest_cookie=AQABAAgAAQBGhh4AAAA%3D&pass_ticket=b2OZ3bwVKQX%2BSf8H5XU5deONOJ5V0hP9ypUMxL24AW%2B8Bkxsw7VOI3D%2FRrD81X1e&wx_header=1


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)-50 + 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画的简图。

1.png

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左右。

2.png
突然发现excel也能够提供类似于heatmap的样式,比较好玩,因此放在这儿。

3.png

9. 递归地删除空目录

find . -depth -type d -empty -exec rmdir -v {} \;

由此引申,可以递归地删除分析过程中产生的中间文件,以节省存储空间。

10. 监控文件的变化,特别是当大文件比对,然后你想看看output到底有没有变化时。

watch  ls  -al  your.sam


欢迎到微信公众号订阅我们
生信媛
bio_sxy




上一篇:基因组分析:共线性作图
下一篇:我谨慎地使用了rm -fr,为什么还是导致了灾难?
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-21 10:03 , Processed in 0.051425 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.