跑BWA比对测试一下酷睿I9的CPU

Posted on by

最近给学员买了一个CPU，非常的霸气：

image-20191023185532496

查看测试数据，配对的肿瘤外显子数据，fq文件如下：

```3.5G 20:10 NPC10F-N_1.fastq.gz
3.6G 20:09 NPC10F-N_2.fastq.gz
3.2G 20:09 NPC10F-T_1.fastq.gz
3.3G 20:08 NPC10F-T_2.fastq.gz
2.7G 20:08 NPC15F-N_1.fastq.gz
2.8G 20:07 NPC15F-N_2.fastq.gz
2.8G 20:07 NPC15F-T_1.fastq.gz
2.9G 20:06 NPC15F-T_2.fastq.gz
2.8G 20:06 NPC29F-N_1.fastq.gz
2.9G 20:05 NPC29F-N_2.fastq.gz
2.4G 20:05 NPC29F-T_1.fastq.gz
2.5G 20:04 NPC29F-T_2.fastq.gz
2.0G 20:04 NPC34F-N_1.fastq.gz
2.0G 20:03 NPC34F-N_2.fastq.gz
2.2G 20:03 NPC34F-T_1.fastq.gz
2.3G 20:03 NPC34F-T_2.fastq.gz
2.1G 20:02 NPC37F-N_1.fastq.gz
2.1G 20:02 NPC37F-N_2.fastq.gz
2.2G 20:02 NPC37F-T_1.fastq.gz
2.2G 20:01 NPC37F-T_2.fastq.gz


都是标准的肿瘤外显子测序数据，T按照道理应该是比N测序数据量大，但是这些样本并不是。

使用**简单的shell脚本**构建配置文件如下：

NPC10F-N /data//NPC10F-N_1.fastq.gz /data//NPC10F-N_2.fastq.gz
NPC10F-T /data//NPC10F-T_1.fastq.gz /data//NPC10F-T_2.fastq.gz
NPC15F-N /data//NPC15F-N_1.fastq.gz /data//NPC15F-N_2.fastq.gz
NPC15F-T /data//NPC15F-T_1.fastq.gz /data//NPC15F-T_2.fastq.gz
NPC29F-N /data//NPC29F-N_1.fastq.gz /data//NPC29F-N_2.fastq.gz
NPC29F-T /data//NPC29F-T_1.fastq.gz /data//NPC29F-T_2.fastq.gz
NPC34F-N /data//NPC34F-N_1.fastq.gz /data//NPC34F-N_2.fastq.gz
NPC34F-T /data//NPC34F-T_1.fastq.gz /data//NPC34F-T_2.fastq.gz
NPC37F-N /data//NPC37F-N_1.fastq.gz /data//NPC37F-N_2.fastq.gz
NPC37F-T /data//NPC37F-T_1.fastq.gz /data//NPC37F-T_2.fastq.gz


这个取决于每个人的脚本习惯，我的脚本是；

ls /data/project/DNA/tumor-NT-NPC-WES/fastqData/_1.fastq.gz > 1
ls /data/project/DNA/tumor-NT-NPC-WES/fastqData/2.fastq.gz > 2
cut -d”/“ -f 7 1 |cut -d”“ -f 1 > 0
paste 0 1 2 > config.clean


### 一个简单的bwa比对脚本

虽然简单，但是技术含量有点高！

```sh
#!/bin/bash
# usage: for i in {0..9};do (nohup bash step1-bwa.sh ./ config.fq 10 $i & );done
analysis_dir=$1
config_file=$2
number1=$3
number2=$4
INDEX=/data/bigbiosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38
cat $config_file |while read id
do
 arr=($id)
 fq1=${arr[1]}
 fq2=${arr[2]}
 sample=${arr[0]}

if((i%$number1==$number2))
 then

#####################################################
################ Step 1 : Alignment #################
#####################################################
if [ ! -f ok.step1-Alignment.$sample.status ]; then
 start=$(date +%s.%N)
 echo bwa $sample `date` >> $sample.log
 bwa mem -t 5 -M -R "@RG\tID:$sample\tSM:$sample\tLB:WES\tPL:Illumina" $INDEX $fq1 $fq2 1>>$sample.sam 2>>/dev/null
 if [ $? -eq 0 ]; then
 echo "bwa succeed" $sample >> $sample.log
 touch ok.step1-Alignment.$sample.status
 else
 echo "failed" $sample >> $sample.log
 fi
 echo bwa $sample `date` >> $sample.log 
 dur=$(echo "$(date +%s.%N) - $start" | bc)
 printf "Execution time for BWA : %.6f seconds" $dur >> $sample.log
 echo >> $sample.log

fi

fi
 i=$((i+1))

done

批量提交的用法是：for i in {0..9};do (nohup bash step1-bwa.sh ./ config.fq 10 $i & );done 很简单：

可以看到全部的10个样本的bwa比对进程都已经开始啦：

image-20191023184635485

输出的文件如下：

$ls -lh *sam
-rw-rw-r-- 1 vip1 vip1 33G 10月 23 19:36 NPC10F-N.sam
-rw-rw-r-- 1 vip1 vip1 30G 10月 23 19:33 NPC10F-T.sam
-rw-rw-r-- 1 vip1 vip1 26G 10月 23 19:30 NPC15F-N.sam
-rw-rw-r-- 1 vip1 vip1 26G 10月 23 19:25 NPC15F-T.sam
-rw-rw-r-- 1 vip1 vip1 27G 10月 23 19:20 NPC29F-N.sam
-rw-rw-r-- 1 vip1 vip1 25G 10月 23 19:15 NPC29F-T.sam
-rw-rw-r-- 1 vip1 vip1 2.3G 10月 23 18:43 NPC34F-N.sam
-rw-rw-r-- 1 vip1 vip1 23G 10月 23 19:15 NPC34F-T.sam
-rw-rw-r-- 1 vip1 vip1 21G 10月 23 19:11 NPC37F-N.sam
-rw-rw-r-- 1 vip1 vip1 22G 10月 23 19:12 NPC37F-T.sam

可以看到9个样本早就成功了，但是其中一个样本明显输出的sam文件小很多，因为我们前面脚本写的好，所以日志都在，简单检查一下：

$grep Execution *log
NPC10F-N.log:Execution time for BWA : 3506.270858 seconds
NPC10F-T.log:Execution time for BWA : 3247.848839 seconds
NPC15F-N.log:Execution time for BWA : 3101.098587 seconds
NPC15F-T.log:Execution time for BWA : 2815.466299 seconds
NPC29F-N.log:Execution time for BWA : 2473.056571 seconds
NPC29F-T.log:Execution time for BWA : 2183.428048 seconds
NPC34F-N.log:Execution time for BWA : 251.685388 seconds
NPC34F-T.log:Execution time for BWA : 2184.577157 seconds
NPC37F-N.log:Execution time for BWA : 1976.718491 seconds
NPC37F-T.log:Execution time for BWA : 2000.084872 seconds

可以看到，对样本NPC34F-N来说，运行了才251秒，也就是说仅仅是载入了bwa的参考基因组而已，简单推测，应该是该样本的fastq文件有问题，可能是拷贝过程损坏。

最简单的就是看看这些fastq文件走fastqc质控是否成功，代码如下：

 cut -f 3 ../config.fq |xargs fastqc -t 10 -O ./
 cut -f 2 ../config.fq |xargs fastqc -t 10 -O ./

有趣的是，样本都成功进行了fastqc质控啊！后来检查配置，才发现，是因为我们的CPU是36线程，但是我提交了10个任务，每个任务调用5个线程，导致资源分配不足，部分任务被kill了。

​

Comments are closed.