首先我们先回顾一下,上次推文主要介绍了bash编程的基础知识,还有一些基本的代码规范。如果你还没读过上次的推文,请不要犹豫先点击下面的链接。
回顾完之后,这次的推文主要学习如何在bash中写更好的loops,还有一些更加高级的shell bash编程知识。事不宜迟,开始今天的学习。
首先下载一些测试数据:
fastq-dump --split-files -X 10000 SRR1553607fastq-dump --split-files -X 10000 SRR1972917
接着假s为我们需要将reads的长度缩短到20 bp。我们使用 cutadapt-l20
。现在让我们先编写一个糟糕的脚本:
for name in *.fastq; do echo "cutadapt -l 20 $name -o $name.trimmed.fq"done
以便小伙伴们更好地了解发生了什么。我们这里先打印命令而不是执行命令,我们的脚本将产生以下内容:
cutadapt -l 20 SRR1553607_1.fastq -o SRR1553607_1.fastq.trimmed.fqcutadapt -l 20 SRR1553607_2.fastq -o SRR1553607_2.fastq.trimmed.fqcutadapt -l 20 SRR1972917_1.fastq -o SRR1972917_1.fastq.trimmed.fqcutadapt -l 20 SRR1972917_2.fastq -o SRR1972917_2.fastq.trimmed.fq
为什么说这个循环(loop)是一个糟糕的例子呢?第一,我们通过文件名(*fq)进行模式匹配,这样一些不是我们想处理,但是又有相同文件名的文件也会被处理。第二,此代码不断在文件名中添加扩展名,每个生成的文件现在都以我们不期待的结尾 .fastq.trimmed.fq
。
要解决扩展名的问题,我们需要调用更复杂的bash构造,bash的替换运算符 %.*
:
for name in *.fastq; do echo "cutadapt -l 20 $name -o ${name%.*}.trimmed.fq"done
可是现在的代码显得更复杂,可读性更低了。
首先,我们需要养成一个习惯,永远不要在 *
匹配的文件“模式”(例如 *.fastq
或 *.bam
等)上运行命令。因为文件的处理顺序可能与期望的不符。另外运行时可能会增加一些你不想运行的文件;这个糟糕的习惯最终会导致一些棘手的问题。
一个好的习惯是,我们需要整理出我们要处理文件的“根”,换而言之就是数据之间用于独特标识的那一部分。以上面的测试数据为例子,它们的“根“就是:
SRR1553607SRR1972917
将上面的根存进去 ids.txt
中,然后我们使用更好的写命令或者循环的工具 parallel
:
cat ids.txt | parallel echo cutadapt -l 20 {}_1.fastq -o {}_1.trimmed.fqcat ids.txt | parallel echo cutadapt -l 20 {}_2.fastq -o {}_2.trimmed.fq
##结果显示cutadapt -l 20 SRR1553607_1.fastq -o SRR1553607_1.trimmed.fqcutadapt -l 20 SRR1972917_1.fastq -o SRR1972917_1.trimmed.fqcutadapt -l 20 SRR1553607_2.fastq -o SRR1553607_2.trimmed.fqcutadapt -l 20 SRR1972917_2.fastq -o SRR1972917_2.trimmed.fq
这个代码看起来就是更加的清楚明了,该代码根据我们给予的“根”,使用 {}
进行匹配,指明了对应的输入和生成文件。
当我们用编程语言编写一个 for
loop时,我们正在构建一个迭代的命令式:我们要求计算机首先完成一个工作,然后循环到最后。但通过GNU Parallel编写命令时,我们遵循所谓的描述性功能编程。就是,我们尝试用模式描述我们想要的内容,然后让计算机填写该模式并输入完整命令。
GNU Parallel 是一个非常好用文件并行的工具。
假设有一个名为的文件 ids.txt
,其中包含:
ABC
假设我们要输出:
Hello AHello BHello C
通过文件输入:
cat ids.txt | parallel echo Hello {}
在命令行中通过用3个冒号( :::
)来指定输入:
parallel echo Hello {} ::: A B C
最后,当用四个冒号( ::::
)分隔时,您也可以在文件末尾传递文件
parallel echo Hello {} :::: ids.txt
上面的每个命令产生相同的输出。
假如我们获取所有的排列组合:
parallel echo Hello {1} and {2} ::: A B ::: 1 2
Hello A and 1Hello A and 2Hello B and 1Hello B and 2
获取一一对应的组合,使用 --link
:
parallel --link echo Hello {1} and {2} ::: A B ::: 1 2
Hello A and 1Hello B and 2
更多详细的关于GNU parallel的内容,可以查阅我之前的推文:
一个好的脚本是应该自带说明manual的。例如,一个脚本需要运行的参数,参数的使用说明等。
下面给大家一个模板例子:
bash getdata.sh
*** This script needs arguments to work! ***
Usage: getdata.sh PRJNUM COUNT LIMT
Parameters: PRJNUM = SRA Bioproject number COUNT = how many sequencing runs to download LIMIT = how many reads to extract per sequencing run
Example: getdata.sh PRJN2234 1000 5
下面是如何编写该说明manual的代码:
if [ $# -eq 0 ]; then echo echo "*** This script needs arguments to work! ***" echo echo "Usage:" echo " getdata.sh PRJNUM COUNT LIMT" echo echo "Parameters:" echo " PRJNUM = SRA Bioproject number" echo " COUNT = how many sequencing runs to download" echo " LIMIT = how many reads to extract per sequencing run" echo echo "Example:" echo " getdata.sh PRJN2234 1000 5" echo exit 1fi
Bash有一个输入流( stdin
)和两个输出流( stdout
和 stderr
)。
通常命令的输出将进入标准输出( stdout
),错误消息将变为标准错误( stderr
)。
默认情况下,两者stdout和stderr都被定向到终端。例如,我可以输入:
ls * foo > B.txt
因为f不存在它输出:
ls: foo: No such file or directory
更加好的方式是使用 2>
,将标准错误存储起来:
ls * foo > B.txt 2> err.txt
如果你遇到错误,则可以调查错误信息文档以获取消息。
通常,我们必须在bash中操作文件名以删除其中的各个部分。也许我们想要删除目录名称,或者仅保留文件名,或者仅保留不带扩展名的文件名,或者删除扩展名等等。
下面让我看一些例子:
FILE=/A/B/C.txt.gzecho $FILE
如预期打印:
/A/B/C.txt.gz
从名称中删除目录,并仅使用basenameshell命令保留文件名:
FILE=/A/B/C.txt.gzNAME=$(basename ${FILE})echo $NAME
打印:
C.txt.gz
要切断最右边的扩展名:
FILE=/A/B/C.txt.gzCHOP=${FILE%.*}echo $CHOP
它将打印
/A/B/C.txt
现在只获取扩展名:
FILE=/A/B/C.txt.gzCHOP=${FILE##*.}echo $CHO
它打印:
gz
用反引号将其括起来:
VALUE=`ls -1 | wc -l`echo "The number of files is $VALUE"
要将默认值分配给变量,请使用以下结构:
FOO=${VARIABLE:-default}
例如,要将 LIMIT
变量设置为第一个参数, $1
或者 1000
默认值如果未指定该参数:
LIMIT=${1:-1000}
编写一个脚本的最好的办法是先将需要运行的代码打印出来,而不是直接运行所有的代码:
echo fastq $SOMETHING
将每一步的命令打印到屏幕可以让我们更加直观的检查每一行的代码。如果整个流程的代码看起来都没问题,就ji执行命令,然后bash再次将它们通过管道传递给命令。
今天的学习就到这里结束了,希望本推文对大家有所帮助。
本文整理参考于:biostar中writing-better-scripts的内容,有条件的小伙伴可以自行购买和下载。https://www.biostarhandbook.com/books/scripting/writing-better-scripts.html