在使用Linux做NGS数据处理的过程中,批量处理是提高工作效率的好方法。
做好配置文件
简单地讲,我们需要先制作一个配置文件,这里命名为config,然后config文件如以下所示:
$ cat config
SRR1039510_1.fastq.gz SRR1039510_2.fastq.gz
SRR1039511_1.fastq.gz SRR1039511_2.fastq.gz
SRR1039512_1.fastq.gz SRR1039512_2.fastq.gz
命令脚本文件
然后把我们的代码写入到一个shell脚本中,这里命名为qc.sh(注:代码中的$1是后面要传入的config)
$ cat qc.sh
#!/bin/bash
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
trim_galore -q 25 --phred33 \
--length 36 --stringency 3 --paired \
-o ./ $fq1 $fq2
done
提交至后台
最后再提交至后台
nohup bash qc.sh config &
最后的最后,要学会通过top查看命令是否成功提交了,如果提交成功,服务器会一个一个地处理数据,这样我们就可以忙别的事情了,等到数据处理得差不多再看处理结果。
进阶版
假如你觉得上面一个一个地处理太!慢!了!那你可以看看下面的进阶处理方法
做好配置文件
$ cat config
SRR1039510_1.fastq.gz SRR1039510_2.fastq.gz
SRR1039511_1.fastq.gz SRR1039511_2.fastq.gz
SRR1039512_1.fastq.gz SRR1039512_2.fastq.gz
......这里省略若干行......
命令脚本文件
$ cat qc.sh
#!/bin/bash
number1=$2
number2=$3
cat $1 | while read id
do
if((i%$number1==$number2))
then
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
trim_galore -q 25 --phred33 \
--length 36 --stringency 3 --paired \
-o ./ $fq1 $fq2
fi ## end for number1
i=$((i+1))
done
提交至后台
for i in {0..2}
do
(nohup bash qc.sh config 3 $i 1>log.$i.txt 2>&1 & )
done
最后的最后,要学会通过top查看命令是否成功提交了,如果提交成功,服务器会批量处理数据,向这里的例子,每次就同时处理3个数据了,当然前提是服务器的资源足够。
错误案例
配置:
$ cat config
SRR1039510 SRR1039510_1.fastq.gz SRR1039510_2.fastq.gz
SRR1039511 SRR1039511_1.fastq.gz SRR1039511_2.fastq.gz
SRR1039512 SRR1039512_1.fastq.gz SRR1039512_2.fastq.gz
......这里省略若干行......
脚本(有问题):
$ cat qc.sh
#!/bin/bash
number1=$2
number2=$3
cat $1 | while read id
do
if [ ! -f ok.trim.$sample.status ]
then
touch ok.trim.$sample.status
echo "start trim for $sample" `date`
arr=(${id})
sample=${arr[0]}
fq1=${arr[1]}
fq2=${arr[2]}
trim_galore -q 25 --phred33 \
--length 36 --stringency 3 --paired \
-o ./ $fq1 $fq2
echo "end trim for $sample" `date`
fi
done
提交命令:
for i in {0..2}
do
(nohup bash qc.sh config 3 $i 1>log.$i.txt 2>&1 & )
done
top之后仍然可以看到任务在运行,3个(不要以为这样就ok了)
我们检查一下生成的文件,发现生成的文件却只有一个样本的
$ ll -th
-rw-rw-r-- 1 hcguo hcguo 1.2G Jun 18 17:45 SRR1039510_1_trimmed.fq.gz
-rw-rw-r-- 1 hcguo hcguo 1.6K Jun 18 17:45 log.0.txt
-rw-rw-r-- 1 hcguo hcguo 1.6K Jun 18 17:45 log.2.txt
-rw-rw-r-- 1 hcguo hcguo 1.6K Jun 18 17:45 log.1.txt
-rw-rw-r-- 1 hcguo hcguo 544 Jun 18 17:39 SRR1039510_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 hcguo hcguo 0 Jun 18 17:39 ok.trim..status
仔细看top命令时,发现3个命令其实处理的是同一个样本,也就是说我们提交的命令有问题,这个时候我们应该进一步查看一下log日志或者nohup.out,看看问题所在(这里不做演示了)