如何寫一個腳本（附送一個腳本）

為什麼要寫腳本

為了把時間放在更有意義的事情上

生信人每天會做大量的重複的操作，例如一個流程會用不同的參數跑好幾遍找合适的參數，不斷把一個格式的資料轉換成另一種格式。

是以從某種程度而言，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

最後，歡迎到我的星球和我探讨

我的知識星球

如何寫一個腳本（附送一個腳本）

為什麼要寫腳本

前期準備

如何寫腳本

開頭怎麼寫

中間怎麼寫

最後說說腳本怎麼用

繼續閱讀

1、Linux 指令行使用技巧

VBA調用cmd指令行下執行的指令

在目前位置打開指令行視窗的技巧

unit 1 - redhat Enterprise 8.0 Linux 指令行使用技巧

Windows指令行中使用SSH連接配接Linux

Linux下指令行中的複制和粘貼

1.Linux指令行使用技巧

spec檔案詳解

HK-2000資料采集儀資料庫操作說明

終端環境之tmux

查找檔案中的字元串

拒絕使用者登入:/bin/false和/usr/sbin/nologin

Shell程式設計——sort排序、uniq忽略重複、tr替換壓縮删除、cut指定删除字段、正規表達式元字元sort 指令uniq 指令tr 指令cut 指令正規表達式

Linxu常用指令技巧彙總

《Linux指令行與Shell腳本程式設計大全第2版.布盧姆》pdf

ACS基本配置-權限等級管理