查看原文
其他

Linux shell trick for bioinformatics 系列文章之一

2016-12-25 天地本无心 生信媛

    由一个纯生物背景的人,转到生物信息学领域之后。刚开始学会写脚本的时候,仿佛开启了一扇通往新世界的大门。好多以前觉得很难完成的事,写个脚本轻松完成。随着时间的推移,攒下了越来越多的小脚本。伴随着技能树的拓展,知道了很多事根本不用写脚本,直接用shell一行搞定。

    并且用shell来处理这些小事的时候,还有一个好处就是通过管道符将这些linux操作串联以来,形成一套pipeline。 诚然,写脚本也可以,但是用linux代码更加符合自己对简单,优雅的审美要求,所以告别那些无用的dirty code,拥抱简洁强大的Shell。


1.   将一个文件按行倒序排列,即将第一行变为倒数第一行,第二行变为倒数第二行,一次类推。我们有一个稍微笨一些的方法:


awk 'BEGIN {x=0}{x=x+1}{print x,$0}'  yourfile.txt | sort -nr -k 1| sed 's/^.* //g' |le 

或者我们采用一个优雅至极的方法;

tac  yourfile.txt

就是linux常用命令cat反过来就是倒序输出。

      

2.  对于Rfam数据库下下来的noncoding RNAs的fasta文件,注释信息就在fasta文件的header里面,如果我们想要提取所有的sequence的注释信息。


grep "^>" rfam.fasta | tr -d ">" |le


3.  用Trinity坐组装的之后得到的fasta格式的转录组文件,我们要做注释的话,有些header可能过长,所以我们要将header给截短些,仅保留最关键的信息。


cut -f 1 -d " " trinity.fasta |le


4.  统计fastq的所有碱基的数目。


awk 'BEGIN{sum=0;}{if(NR%4==2){sum+=length($0);}}END{print sum;}' your.fq


5. SAM文件中根据某一行的数值进行筛选,并且保留header.  以拟南芥这种模式植物为例,和拟南芥的记忆组比对得到的SAM文件有9行header,我们想保留header,再根据某行值进行过滤。


awk '$3<1000 || NR <= 9' your.SAM


6.  将fastq格式文件转化为fasta格式。


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


7.  当我们得到基因的注释之后,得到了two column的基因注释table, 但是里面有很多NA值。但我们不想看到NA,怎么办呢。方法很多。


perl -ne 'print unlesss /NA/'  file.txt

awk  '{if(!/NA/){print $0}}'  file.txt

awk '$0 !~ /NA/ {print}'  file.txt

grep -v "NA"  file.txt

sed -e '/*NA*/d'  file.txt


8.  递归地查找当前目录下所有以"fasta"为后缀的文件,统计其行数。如果想统计fasta文件的所有碱基数目,请参考统计fastq碱基数目那一条。


find . -name "*.fasta " | xargs wc -l


9.  当你有两个gene list的文件,要找出仅在gene.file2中存在的行,不在gene.file1中出现的gene用来做GO分析。


grep -Fxv -f gene.file1 gene.file2


10.  在某个目录下有很多文件,你想看看你最感兴趣的基因名字出现在哪个文件里,文件很多,并且子文件夹还有很多文件,即在一个目录下递归地查找匹配某字符串的文件。怎么办呢。


grep -Hrn "shengxinyuan" .




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

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