我前面写到了生信分析人员如何入门linux和perl,后面还会写R和python的总结,但是在这中间我想插入一个脚本实战指南。其实在我前两篇日志里面也重点提到了学习编程语言最重要的就是实战了,也点出了几个关键词。在实际生物信息学数据处理中应用perl和linux,可以借鉴EMBOSS软件套件,fastx-toolkit等基础软件,实现并且模仿该软件的功能。尤其是SMS2/exonerate/里面的一些常见功能,还有DNA2.0 Bioinformatics Toolbox的一些工具。如果你这些名词不懂,请赶快谷歌!!! 它们做了什么,输入文件是什么,输出文件是什么,你都可以用脚本实现!
Tag Archives: 脚本
没必要学shell进阶语法
因为大部分生物信息学软件都是linux版本的,所以生物信息学数据分析工作者必备技能就是linux,但是大部分人只是拿他当个中转站,我以前也是,直到接触了大批量的任务,自动化流程,才明白这里面的水太深了,不过无所谓,凭我个人的观点,其实shell的进阶语法真的不必要!
我也不想去搞明白操作符两边是否加空格的区别是什么了。
如果,还有其它功能需要实现,我们可以把脚本写的更负载,纯粹的用shell,需要搜索更多的shell技巧。
[perl]
## perl nohup.pl deep_count.sh 0
## perl nohup.pl deep_count.sh 1
## perl nohup.pl deep_count.sh 2
$i=1;
open FH,$ARGV[0];
while(<FH>){
chomp;
next unless $.%3==$ARGV[1];
$cmd="nohup $_ &";
print "$cmd\n";
system($cmd);
sleep(10800) if $i%5==4;
$i++;
#exit;
}
[/perl]
Shell里面的各种括号的区别
[],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!
[[ ]] double brackets
(())Double parentheses
{{}}double curly brackets
我们必须要记住的是下面
[] 相当于test,作逻辑判断
$( ) 与` ` (反引号) 都是用来做命令替换用
${ } 吧... 它其实就是用来作变量替换用的啦
(())就是用来计算的,相当于expr函数。
http://tldp.org/LDP/abs/html/index.html
我们首先看看一对的括号
首先[]是用来逻辑判断的,必须有空格
if [ -f binom.py ]
then
echo 'binom.py exists'
fi
或者
nub=$((i%4))
#echo $nub
if [ $nub == 0 ];then
echo "we need to sleep 4 hours"
sleep 14000
fi
这个[]操作符等价于test函数
if test $1 -gt 0
then
echo "$1 number is positive"
fi
但是都必须有空格!!!
参考:http://www.freeos.com/guides/lsst/ch03sec02.html
关于shell的test操作符还有很多http://tldp.org/LDP/abs/html/fto.html
( ) 将command group 置于 sub-shell 去执行,也称 nested sub-shell。
{ } 则是在同一个 shell 内完成,也称为non-named command group。
补充一个: {} 还可以做变量扩展 {5..9} 或者 {abcd}e, 自己运行一下就知道效果啦
这两个差异很小,而且一般用不着,就不讲了。
那么这一对的括号加上了$符号后又变成了上面鬼东西呢?
当然,只有:$( ) 与${ }才是合法的。
在 bash shell 中,$( ) 与` ` (反引号) 都是用来做命令替换用(command substitution)的。
在操作上,用$( ) 或` ` 都无所谓,用$( )的优点是:
1, ` ` 很容易与' ' ( 单引号)搞混乱,尤其对初学者来说
2, 在多层次的复合替换中,` ` 须要额外的跳脱( \` )处理,而$( ) 则比较直观
再让我们看${ } 吧... 它其实就是用来作变量替换用的啦。
一般情况下,$var 与${var} 并没有啥不一样。
但是用${ } 会比较精确的界定变量名称的范围,比方说:
[code][/code]
$ A=B
$ echo $AB
还可以用来截取变量,这个就很多花样啦
# 是去掉左边(在鉴盘上# 在$ 之左边)
% 是去掉右边(在鉴盘上% 在$ 之右边)
单一符号是最小匹配﹔两个符号是最大匹配
然后我们看看两对的括号:
nub=$((i%4)) 等价于$nub=`expr $i % 1` ;
((i++)) 等价于$i=`expr $i + 1` ;
所以(())就是用来计算的,而且里面的变量不需要$来标记啦
(在 $(( )) 中的变量名称,可于其前面加$ 符号来替换,也可以不用)
在(())前面加上$只是为了把计算结果给保存而已。
而两个中括号和两个大括号都是不合法的!
脚本作业-解读NCBI的ftp里面关于人的一些基因信息
为了感谢大家对我博客的关注,我在这里发布一个作业,适合菜鸟做的。里面有十几个类似的问题,大家可以下载数据自行处理,如果是问这些问题,我优先回答!
NCBI的ftp里面关于人的一些基因信息
我在NCBI的ftp服务器里面下载了这些数据,时间是2015年,大多是hg19系列的,文件名如下:
CDS.fa 这个是ensembl中人的CDS碱基序列文件,hg38
entrez2go.gene 这个是有go注释的基因情况,有一万八的基因都有go注释
entrez2name.gene 这个是NCBI的entrez ID号对应着基因名的文件
entrez2pubmed.gene 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号
entrez2refseq2ensembl.gene 这个是NCBI的entrez ID号对应着基因名的refseq的ID号和ensembl数据库的ID号
human_gene_info这个是基因的详细信息,包括基因的起始终止点坐标等等
Protein.fa 这个是ensembl中人的蛋白的氨基酸序列文件,有十万多个蛋白hg38
ref2ensembl.txt 这个是基因名的refseq的ID号和ensembl数据库的ID号
自行去NCBI的ftp服务器里面下载这些数据。
然后好好熟悉这些数据信息,回答一下几个问题:
人总的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。
CD分子的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。它们有没有氨基酸偏好性??
MHC系列基因信息?CCL系列基因信息如何?CXCL系列信息如何?或者你感兴趣的基因家族信息?
现在研究最热门的基因是什么?发表文章最多的前十个基因是什么?
基因长度情况如何?最长的基因多长?最短的基因多少bp,可靠吗?
蛋白质长度情况如何?
每条染色体的基因分别情况?基因在染色体那个地方分别最多?
请用图形展示你的结论!!!
如果你能回答以上问题,证明你的脚本水平不错了。
如果找不到我,看旁边的公告,加入生信菜鸟群,我就在里面!!!
Linux基础之shell脚本的批处理
脚本类似于下面的样子,大家可以读懂之后就仿写
for i in *sra
do
echo $i
/home/jmzeng/bio-soft/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump --split-3 $i
Done
这个脚本是把当前目录下所有的NCBI下载的sra文件都加压开来成测序fastq格式文件
有这些数据,分布在不同的目录,如果是写命令一个个文件处理,很麻烦,如果有几百个那就更麻烦了,所以需要用shell脚本
这样只需要bash这个脚本即可一次性处理所有的数据
还有很多类似的脚本,非常简单的
for i in *fq
do
echo $i
bowtie2 -p 13 -x ../../RNA.fa -U $i -S $i.sam
done
for i in */accepted_hits.bam
do
echo $i
out=`echo $i |cut -d'/' -f 1`_clout
samtools mpileup -guSDf /home/immune/refer_genome/hg19/hg19.fa $i | bcftools view -cvNg - >snp-vcf/$out.vcf
done
while read id
do
echo $id
out=`echo $id |cut -d'/' -f 2`
reads=`echo $id |cut -d'/' -f 3|sed 's/\r//g'`
tophat2 -p 13 -o $out /home/immune/refer_genome/hg19/hg19 $reads
done <$1
等等