为什么要写脚本

为了把时间放在更有意义的事情上

生信人每天会做大量的重复的操作，例如一个流程会用不同的参数跑好几遍找合适的参数，不断把一个格式的数据转换成另一种格式。

所以从某种程度而言，wet（和各种试剂打交道）和dry（和计算机打交道）实验人员是差不多的，都是有一个想法，然后通过实验进行验证。验证的过程中有一部分的探索性的，另一部分是大量的重复性工作。

dry相对于wet实验的一个好处就在于可以通过写脚本让计算机去自动化完成那些重复性的工作。另外，虽然前面使用管道能够很方便地处理数据，但是由于我们往往不会记录这些过程，很容易导致实验不可重复，出了问题找不到原因，因此使用脚本还能对操作进行记录，提高鲁棒性（robust）。

前期准备

由于这个脚本是用于处理mapping-by-sequencing，为了便于实践，推荐你们创建如下文件夹， 并下载SHOREmap提供的demo数据

# 进入家目录

cd ~

# 创建项目文件夹目录

mkdir mbs-2017-4 && cd mbs-2017-4

mkdir -p {data/{seq,reference},analysis,script}

cd data/seq

# 下载测序数据

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads1.fq.gz &

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads2.fq.gz &

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads1.fq.gz &

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads2.fq.gz &

# 解压缩

ls | while read id; do gunzip $id;done &

# 下载拟南芥参考基因组和注释文件

cd ../reference

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_chr_all.fas &

wget -c -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_GFF3_genes.gff &

# 准备样本信息文件

cd ../../script

# 注意中间的空格为"\t"

vi sample.txt

BC.fg.reads1 R1 data/seq/BC.fg.reads1.fq.gz

BC.fg.reads2 R2 data/seq/BC.fg.reads1.fq.gz

BC.bg.reads1 R1 data/seq/BC.bg.reads1.fq.gz

BC.bg.reads2 R2 data/seq/BC.bg.reads1.fq.gz

如何写脚本

我会根据我写的一个脚本讲解，然后你们可以根据自己的需求优化这个脚本，写出属于自己的脚本。 一般第一个脚本我们都需要运行”hello world“，用来向编程之神祷告，让自己能够更容易的学习编程。在Linux中新建一个hello_world.sh，然后输入如下内容。

#! /bin/bash

echo "hello world"

运行方式有两种

# 方法1

bash hello_world.sh

# 方法2

chmod u+x hello_world.sh

./hello_world.sh

下面这个脚本是为了帮助我自动完成重测序的工作。由于目前用的服务器性能堪忧，8G数据的alignmeng+snp calling的工作需要跑一天时间。按照初学时期的做法，就是手动输入命令完成一个阶段的任务，然后继续下一个阶段。这极大降低了效率，万一程序在半夜完成，你无法马上继续下一阶段；而且不可能时刻盯着，你还要其他活要干，所以我写了下面这个初级脚本用来完成任务。 这个脚本存放在script目录下。

#! /bin/bash

set -u

set -e

set -o pipefail

# set work path

PATH=/bin:/usr/bin:/usr/sbin:/usr/local/bin:/usr/local/sbin:~/bin

work_dir=$1

sample_info=$2

reference=${work_dir}/data/reference/TAIR10_chr_all.fas

# alignment

bwa index $reference

sample_names=($(cut -f 1 "$sample_info" | uniq))

mkdir -p ${work_dir}/analysis/align

cd ${work_dir}/data/seq

for sample in ${sample_names[@]}

do

# create an output file from the sample name

result_file="${sample}.sam"

bwa mem $reference ${sample}_R1.fastq ${sample}_R2.fastq > ${work_dir}/analysis/align/$result_file

#echo "$result_file"

done

#snp calling

cd ${work_dir}/analysis/align

samfiles=$(ls *sam)

echo $samfiles

mkdir -p ${work_dir}/analysis/snp

for file in ${samfiles[@]}

do

output=${file%.*}

samtools view -b -o ${output}.bam ${file}

samtools sort -o ${output}.sorted.bam ${output}.bam

samtools index ${output}.sorted.bam

samtools mpileup -u -t DP -f $reference ${output}.sorted.bam | bcftools call -vm -Ov > ../snp/${output}.vcf

done

# convert vcf4.2 to vcf 4.1

cd ${work_dir}/analysis/snp

infiles=$(ls *vcf)

for infile in ${infiles[@]}

do

bcftools view --max-alleles 2 -O v ${infile} | \

sed "s/##fileformat=VCFv4.2/##fileformat=VCFv4.1/" | \

sed "s/(//" | \

sed "s/)//" | \

sed "s/,Version=\"3\">/>/" | \

bcftools view -O z > ${infile%%.*}.dg.vcf.gz

done

开头怎么写

#!/bin/bash

set -e

set -u

set -o pipeline

所谓良好的开端是成功的一半，拖延症主要原因就在于无法开始。无论你知不知道你要写啥，建立一个any_name.sh(any name就是你随便取一个名字，类似于any key）的文件然后把下面4行代码写进去就对了。

第1行是shebang，是告知系统用什么解释器来运行这个程序（假设它可执行）。第2~3行的目的是让脚本更加“敏感“，默认状态下一行shell脚本会一行一行往下运行，即便出错也不会终止，这3行代码就可以及时终止，避免rm -rf $NULL/*的惨剧发生。

中间怎么写

写完了开头之后，我们就根据自己的预期，将日常中不断重复的工作写到脚本中。但是为了脚本具有普遍适用性，所以要用到

变量和命令参数

条件语句

循环语句和文件名替换

变量和命令参数

# set work path

PATH=/bin:/usr/bin:/usr/sbin:/usr/local/bin:/usr/local/sbin:~/bin

work_dir=$1

sample_info=$2

reference=${work_dir}/data/reference/TAIR10_chr_all.fas

所谓变量就是“指鹿为马”，比如我说reference=/genome/A.thalina/TAIR10_chr_all.fa,那么当你echo "$reference"的时候就会显示后面的内容，而不是referece。

$ reference=/genome/A.thalina/TAIR10_chr_all.fa

$ echo $reference

/genome/A.thalina/TAIR10_chr_all.fa

命令参数就是命令后接的部分，例如cat some.file的命令参数就是some.file。我们可以利用这个功能从外界传入自定义的内容。shell脚本以$0,$1,$2,$3..分别代表命令自身，第一个参数，第二个参数，第三个参数等。 新建一个argument_test.sh，输入如下内容。

#! /bin/bash

echo "command is $0"

echo "first argument is $1"

echo "second argument is $2"

echo "third argument is $3"

运行结果：

$ bash argument_test.sh A B C

command is argument_test.sh

first argument is A

second argument is B

third argument is C

这里我传入了我项目的根目录和存放样本信息的文件。

条件语句

在我的脚本其实并没有用到条件语句，毕竟只是第一版，能用就行。正常情况下最好加上条件语句，提高脚本的适用性。

if [条件]

then

commands

else

comgmands

fi

[条件]可以用bash下的命令，例如

# 单个条件

if grep "pattern" some.file >/dev/null

then

echo "found 'pattern' in some.file"

fi

# 多个条件

# 或 || ; 与&& ;非 !

if grep 'pattern1' some.file > /dev/null &&

grep 'pattern2' some.file > dev/null

then

echo "found 'pattern1' and 'pattern2' in some.file"

fi

除了bash下的命令外，还可以用test或[]判断条件是否成立，如下

# -f 判断是否为文件

if test -f some.file

then

....

fi

# 或[] 记住里面的左右一定要有空格

if [ -f some.file ]

then

....

fi

# 这里的或是-o, 与是-a ,非还是!

if [ -f some.file -a -x some.fie ]

then

....

fi

《鸟叔的Linux的私房菜》第三版的380页中提高了许多判断方式，一般常用的如下：

测试的标志

代表意义

-d

是否为文件夹

-f

是否为文件

-e

文件是否存在

-h

是否为连接

-w

能否写入

-x

能否执行

str1 = str2

字符是否相同

str1 != str2

字符是否不相同

-eq

数值是否相等

-ne

数值是否不等

-lt

小于

-gt

大于

-le

小于等于

-ge

大于等于

循环语句和文件名替换

一般而言，创建一个pipline处理文件需要三个步骤：

选择目标文件

使用指定命令处理目标文件

追踪处理后的输出文件

而循环语句就是在选择目标文件后，用于重复相同的指令处理目标文件，还可以记录输出文件。

其中选择目标文件有两种方式：1. 提供记录目标文件信息的文本，2.从含目标文件的文件夹内筛选。 方法1：

# 假设你现在在mbs-2017-4/script下

sample_info=sample.txt

sample_names=($(cut -f 1 "$sample_info" ))

echo "${sample_names[@]}"

BC.fg.reads1 BC.fg.reads2 BC.bg.reads1 BC.bg.reads2

说明 ${sample_names[@]}的表示显示所有变量内容，我们可以选择特定的元素，将@替换成0,1,2,3（以0开始），此外 ${#sample_names[@]}表示共有多个元素， ${!sample_names[@]}表示各个元素的索引。

在将目标文件赋予变量后，就可以使用for循环对目标文件应用相应的命令。在我的脚本中，因为bwa比对软件需要两个输出文件，因此就需要对第一行的内容去重，然后在for循环中添加后缀。也就是

# alignment

sample_names=($(cut -f 1 "$sample_info" | cut -d '.' -f 1,2 | uniq))

mkdir -p ${work_dir}/analysis/align

cd ${work_dir}/data/seq

for sample in ${sample_names[@]}

do

# create an output file from the sample name

result_file="${sample}.sam"

bwa mem $reference ${sample}.read1.fastq ${sample}.read2.fastq > ${work_dir}/analysis/align/$result_file

echo "$result_file"

done

方法2比较粗暴，直接利用通配符列出文件。

sample_names=$(ls *.fq)

for sample in ${sample_names[@]}

do

最后说说脚本怎么用

这个脚本的主要目的是完成重测序的任务，目前还不太完善，需要注意一下几点：

比对的参数需要根据具体情况修改，

要求具有特定的目录构造

缺少日志信息输出，不利于调试

对样本信息的格式有一定的要求，需要对sample_names=($(cut -f 1 "$sample_info" | cut -d '.' -f 1,2 | uniq))修改，不够人性化

后续寻找mutation需要其他程序支持。

用法很简单：提供一个样本信息的文本（三列，样本名，R1/R2，所在文件目录），和文件的根目录在哪里，假设你在scripts下有了sample.txt和mbs.sh，则

chmod u+x mbs.sh

./mbas sample.txt /path/to/project

最后，欢迎到我的星球和我探讨

我的知识星球