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

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

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

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比对脚本

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

#!/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比对进程都已经开始啦：

拍照看看工作中的CPU：

输出的文件如下：

$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了。