01

生信分析人员数据处理脚本实战

我前面写到了生信分析人员如何入门linux和perl,后面还会写R和python的总结,但是在这中间我想插入一个脚本实战指南。其实在我前两篇日志里面也重点提到了学习编程语言最重要的就是实战了,也点出了几个关键词。在实际生物信息学数据处理中应用perl和linux,可以借鉴EMBOSS软件套件,fastx-toolkit等基础软件,实现并且模仿该软件的功能。尤其是SMS2/exonerate/里面的一些常见功能,还有DNA2.0 Bioinformatics Toolbox的一些工具。如果你这些名词不懂,请赶快谷歌!!! 它们做了什么,输入文件是什么,输出文件是什么,你都可以用脚本实现!

Continue reading

22

没必要学shell进阶语法

因为大部分生物信息学软件都是linux版本的,所以生物信息学数据分析工作者必备技能就是linux,但是大部分人只是拿他当个中转站,我以前也是,直到接触了大批量的任务,自动化流程,才明白这里面的水太深了,不过无所谓,凭我个人的观点,其实shell的进阶语法真的不必要!

当然,只是我一家之言!
我实在是不想去背诵大括号,小括号,中括号以及双重括号到底区别是什么!
http://www.bio-info-trainee.com/?p=1018  [],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

我也不想去搞明白操作符两边是否加空格的区别是什么了。

if((i%5==0)) 来判断变量是否被一个数整除
i=$((i+1))来表示变量自增。
这些东西真的很诡异!
如果你有qsub,condor等任务提交系统,那么你只需要熟悉他们就可以了,但大部分散兵游勇的生物信息学家并没有集群,所以压根不会接触任务提交系统,就需要些自动化脚本了!
受限制与机器的cpu以及内存数,需要判断提交了多少任务,等待多久再执行,所以会把一个简单的自动化脚本写的很复杂!
比如下面这个脚本:cat >download_hg38_from_UCSC.sh

for i in $(seq 1 22) X Y M;
do echo $i;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg38.fa;
done
rm -fr chr*.fa
可以下载hg38基因组的fasta文件,但是是分染色体一个个下载的!
再比如下面这个,批量做GSEA分析的脚本:
while read id
do
echo $id
gene=`echo $id |awk '{print $1}'`
probe=`echo $id |awk '{print $2}'`
echo $i
do_GSEA $probe $gene; ##这里是我自己定义的一个function,就不贴出来了
if((i%5==0))
then
sleep 10  ##重点就在这里,每次提交的任务有限制,所以需要休息,不然机器的cpu负载太高!
fi
i=$((i+1))
done <$1

如果,还有其它功能需要实现,我们可以把脚本写的更负载,纯粹的用shell,需要搜索更多的shell技巧。

但是事实上并没有这个必要,我们现在有了更方便的脚本语言,比如我所擅长的perl
我写一个nohup提交任务的脚本!
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2

[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,写了很久,总是报错,但是用perl,一分钟我就写完了,所以,最好是用自己熟悉的一种语法最好!
25

Shell里面的各种括号的区别

[],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

[[ ]] double brackets

(())Double parentheses

{{}}double curly brackets

我们必须要记住的是下面

[] 相当于test,作逻辑判断

$( ) 与` ` (反引号) 都是用来做命令替换用

${ } 吧... 它其实就是用来作变量替换用的啦

(())就是用来计算的,相当于expr函数。

参考:http://sayle.net/book/

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` ;

所以(())就是用来计算的,而且里面的变量不需要$来标记啦

(在 $(( )) 中的变量名称,可于其前面加$ 符号来替换,也可以不用)

在(())前面加上$只是为了把计算结果给保存而已。

而两个中括号和两个大括号都是不合法的!

 

01

脚本作业-解读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,可靠吗?

蛋白质长度情况如何?

每条染色体的基因分别情况?基因在染色体那个地方分别最多?

请用图形展示你的结论!!!

 

如果你能回答以上问题,证明你的脚本水平不错了。

如果找不到我,看旁边的公告,加入生信菜鸟群,我就在里面!!!

19

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脚本

Linux基础之shell脚本的批处理254

这样只需要bash这个脚本即可一次性处理所有的数据

Linux基础之shell脚本的批处理282

还有很多类似的脚本,非常简单的

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

 

等等