Sentieon | 应用教程: Sentieon分布模式

# 一、介绍

**本文档描述了如何利用Sentieon®基因组学工具的分片能力将DNAseq®流程分布到多台服务器上；将其他流程（如TNseq®）进行分布遵循相同原则，因为所有Sentieon®基因组学工具都具有相同的内置分布式处理能力。** 这种分布的目标是为了减少流程的总运行时间，以更快地生成结果；然而，这种分布也会带来一些额外的开销，使计算成本增加。

利用分布能力，流程的每个阶段被分成小任务；每个任务处理基因组的一部分，并可以在不同的服务器上并行运行。每个任务生成一个部分结果，需要按顺序合并为最终的单一输出；这种合并需要仔细进行，以确保考虑到边界并生成与没有分片运行的流程相同的结果。
 
 分布的执行框架不在本文档的范围内，用户需要在保持正确的数据依赖关系的同时，分发数据/文件并启动正确的进程。

---

# 二、分片和分片化

我们将基因组分成许多连续且不重叠的部分，每个部分称为一个分片（shard）。每个分片可以定义为单个染色体contig的名称，或者是遵循contig:shard_start-shard_end约定的染色体的一部分。特殊的分片NO_COOR用于所有没有坐标的未映射读取。

Sentieon®二进制文件支持将分片分布到多个服务器，并且可以通过添加一个或多个带参数的分片选项在单个命令中处理多个分片。在单个命令行中使用多个选项时，这些分片需要按照参考染色体列表是连续的；例如，一个命令可以包含一个覆盖chr2结尾的分片和一个覆盖chr3开头的分片，但不能同时包含一个覆盖chr2结尾和chr22开头的分片。`--shard SHARD --shard`
 
 您可以参考附录中的脚本，根据与参考相关联的dict或输入bam文件中的bam文件表头，为基因组创建分片的示例文件。我们建议使用100M作为分片大小。`generate_shards.sh`

---

# 三、分布式处理的数据流和数据/文件依赖关系

按照推荐的工作流程，DNAseq®的流程将一对fastq文件通过以下阶段进行处理：BWA对齐生成sorted.bam，去重生成dedup.bam，BQSR生成recal.table，Haplotyper生成output.vcf.gz文件。如图1所示，说明了这样一个流程的数据流。

![](https://img2024.cnblogs.com/blog/3218664/202407/3218664-20240715174653658-1270344885.png#align=center)

<center>图1 DNAseq®流程的典型数据流程</center>

为了将上述流程分布到多个服务器上，每个阶段被划分为处理分片数据的命令，这些命令需要来自特定分片以及相邻分片的输入；但是，有两个例外情况：Dedup阶段需要来自所有分片上的LocusCollector命令的所有score文件，而Haplotyper阶段需要一个完整的合并校正表。

![](https://img2024.cnblogs.com/blog/3218664/202407/3218664-20240715174706227-1019694470.png#align=center)

<center>图2 以4个分片进行分布式处理的DNAseq®流程的数据流程，不处理未映射的读取</center>

以图2为例，说明了一个以4个分片进行分布的流水线的数据流程；这并不是一个典型的使用案例，因为使用推荐的分片大小为100M碱基会导致需要超过30个分片。在图2的示例中，各个阶段需要以下输入并生成以下输出：

- 分片的LocusCollector阶段（去重1）需要sorted.bam作为输入。该阶段生成一个`i-th part_deduped.bam$shard_i.score.gz`文件。
 
 - 分片的Dedup阶段（去重2）需要sorted.bam输入，以及所有LocusCollector阶段的所有结果。该阶段生成一个`i-th part_deduped$shard_i.bam`文件。
 
 - 分片的QualCal阶段需要`part_deduped$shard_i.bam` 文件，以及 `part_deduped$shard_i+1.bam` 文件和 `part_deduped$shard_i-1.bam` 文件（如果可用）。该阶段产生一个 `part_recal_data$shard_i.table` 文件。
 
 - 在QualCal之后，所有 `part_recal_data$shard_i.table` 文件需要合并为一个单一的校准表文件，用于变异检测。driver 中的 `--passthru` 选项和 QualCal 的` --merge` 选项可用于执行边界感知的合并。
 
-  分片的Haplotyper阶段需要 `part_deduped$shard_i.bam` 文件，以及 `part_deduped$shard_i+1.bam` 文件和 `part_deduped$shard_i-1.bam` 文件（如果可用）；此外，完整的合并后的重校准表需要作为输入。该阶段产生一个 `part_output$shard_i.vcf.gz` 文件。
 
 - 在Haplotyper之后，所有 `part_output$shard_i.vcf.gz `文件需要合并为一个单一的输出 VCF 文件。driver 中的 `--passthru` 选项和 Haplotyper 的 `--merge` 选项用于执行边界感知的合并。
 
 - 如果在流程的任何阶段您需要合并后的输出 BAM 文件，可以使用 util 二进制文件来执行边界感知的合并；您需要在命令中添加 `--mergemode=10` 选项，以便 util merge 不处理读段，而仅按块进行复制。
 
**重要提示：** 在合并结果时，务必按照正确的顺序依次输入部分结果，否则结果将不正确。

---

# 四、分发BWA

上述指令不包含任何关于使用 BWA 分布比对的信息。为了分布 BWA 比对，您可以使用 Sentieon® 工具中包含的 fqidx 工具为输入的 FASTQ 文件创建索引文件；然后您可以提取 FASTQ 文件的特定部分在不同的服务器上进行处理，使用 fqidx 命令的结果作为 BWA mem 的输入；您需要确保 BWA 命令中包含 `-p` 选项，因为 fqidx 的输出在单个输出中包含交错的读段。此流程的结果与在单次运行中执行的结果完全一致。

```
 BWA_K_size=100000000
num_splits=10
SAMPLE="sample_name" #Sample name SM tag in bam
GROUP="read_group_name" #Read Group name RGID tag in bam
PLATFORM="ILLUMINA"
NT=$(nproc) #number of threads to use in computation, set to all threads available

FASTA="/home/regression/references/b37/human_g1k_v37_decoy.fasta"
FASTQ_1="WGS-30x_1.fastq.gz"
FASTQ_2="WGS-30x_2.fastq.gz"
FASTQ_INDEX="WGS-30x.fastq.gz.index"

#create FASTQ indices
sentieon fqidx index -K $BWA_K_size -o $FASTQ_INDEX $FASTQ_1 $FASTQ_2

#get the number of runs that the inputs will be split into given the size
num_K=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f1)
BWA_K_size=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f2)
num_K_splits=$(expr $num_K / $num_splits + 1)
#run multiple BWA on multiple servers
file_list=""
for run in $(seq 0 $((num_splits-1))); do
       region="$((run*num_K_splits))-$((run*num_K_splits+num_K_splits))"
        #send each of these to a different server
        sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
          -K $BWA_K_size -t $NT -p $FASTA \
          "<sentieon fqidx extract -i $FASTQ_INDEX -r $region $FASTQ_1 $FASTQ_2" | \
          sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
        file_list="$file_list sorted_run$run.bam"
done

#merge the results
sentieon util merge -o sorted.bam $file_list
#or perform deduplication on all the intermediate BAM files
file_list_bam=""
for file in $file_list; do file_list_bam="$file_list_bam -i $file"; done
sentieon driver -r $FASTA -t $NT $file_list --algo LocusCollector --fun score_info score.txt.gz
sentieon driver -r $FASTA -t $NT $file_list --algo Dedup --score_info score.txt.gz deduped.bam
```

**当无法创建 FASTQ 索引文件时分布 BWA (版本 201808.02+)**

在无法创建 FASTQ 索引文件的情况下，您可以使用 util fqidx 配合比例选项 `-F m/n`，同时使用 `-K` 选项。这将输入的 FASTQs 划分为 n 个 读段块的拆分，然后从每个拆分中提取每第 m 个元素，以便在不同的服务器上处理，其中 m 从 0 到 n-1。

请注意，此功能仅建议在文件存储速度足够快以支持每个fqidx进程同时读取输入FASTQ文件的IT环境中使用。这在云环境中很常见，或者在具有高带宽NFS系统的本地集群环境中使用。

```
 BWA_K_size=100000000
num_splits=10
SAMPLE="sample_name" #Sample name SM tag in bam
GROUP="read_group_name" #Read Group name RGID tag in bam
PLATFORM="ILLUMINA"
NT=$(nproc) #number of threads to use in computation, set to all threads available

FASTA="/home/regression/references/b37/human_g1k_v37_decoy.fasta"
FASTQ_1="WGS-30x_1.fastq.gz"
FASTQ_2="WGS-30x_2.fastq.gz"

#run multiple BWA on multiple servers
file_list=""
for run in $(seq 0 $((num_splits-1))); do
       #send each of these to a different server
       sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
          -K $BWA_K_size -t $NT -p $FASTA \
         "<sentieon fqidx extract -F $((run))/$num_splits -K $BWA_K_size $FASTQ_1 $FASTQ_2" | \
         sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
        file_list="$file_list sorted_run$run.bam"
done

#merge the results
sentieon util merge -o sorted.bam $file_list
#or perform deduplication on all the intermediate BAM files
file_list_bam=""
for file in $file_list; do file_list_bam="$file_list_bam -i $file"; done
sentieon driver -r $FASTA -t $NT $file_list --algo LocusCollector --fun score_info score.txt.gz
sentieon driver -r $FASTA -t $NT $file_list --algo Dedup --score_info score.txt.gz deduped.bam
```

---

# 五、大规模队列联合调用需要GVCFtyper做分布式处理

针对包含超过1000个样本的大规模联合调用，我们推荐在GVCFtyper中设置`--genotype_mode multinomial`。虽然默认的基因分型模式在较小样本集中理论上更准确，但多项式模式在大样本集中同样准确，并且在大量样本中的扩展性更好。

在运行大量样本（>3000）的联合队列调用时，建议使用类似上述描述的分发技术进行分发；然而，需要考虑以下几点挑战：

- 联合分型的整体运行时间。
   
 - 运行分布式 GVCFtyper 命令的机器资源需求。

- 存储和访问大量GVCF输入文件的物流管理。
 -  输出VCF文件的大小。
 
一般推荐使用Sentieon®的GVCFtyper的分片功能，将不同的基因组片段分别处理在不同的机器上。以下命令通过完整列表访问所有 GVCF 输入，并假设这些输入存放在通用位置，例如 NFS 存储或可通过 S3、HTTP/HTTPS 协议访问的远程对象存储：

```
#individual shard calling in machine 1
sentieon driver -r $FASTA -t $NT --shard $shard1 --shard $shard2 \
  --algo GVCFtyper $gvcf_argument GVCFtyper-shard_1.vcf.gz
#individual shard calling in machine 2
sentieon driver -r $FASTA -t $NT --shard $shard3 \
  --algo GVCFtyper $gvcf_argument GVCFtyper-shard_2.vcf.gz
…
```

**重要提示：**

在处理每个片段时，输入GVCF的列表需要保持一致的顺序，因为最终输出中的样本顺序将取决于输入顺序，而合并要求所有部分文件具有相同的样本顺序。
 
 使用`--shard`选项时，GVCFtyper生成的输出VCF文件不是有效的VCF文件，因此在按如下所述方式完成合并之前，不应直接使用它们。

在处理完所有片段后，您需要运行GVCFtyper的合并命令来合并结果，确保中间VCF输入的顺序与参考基因组一致。输入文件需要在共享位置（例如NFS存储或可通过S3、http/https协议访问的远程对象存储位置）中可用。

```
#merge step
sentieon driver --passthru -t $nt --algo GVCFtyper --merge joint_final.vcf.gz \
  GVCFtyper-shard_1.vcf.gz GVCFtyper-shard_2.vcf.gz …
```

---

# 六、总体运行时间和资源要求

为了减少整体运行时间，建议使用足够的分片，允许将运行时细粒度分布为多个 服务器。分片可以按分片大小或预期创建复杂性。但是，单个分片的最终合并步骤结果无法分发，需要在单个服务器中运行; 这一事实设定了分布式分析所能使用的服务器数量下限，因为合并步骤可能会主导整体运行时间。

GVCFtyper 算法是一种非常有效的算法，其性能很可能受限于 GVCF 输入文件存储位置的 I/O 表现。因此，建议将 GVCF 输入存储在文件中提供 600 MBps 传输速率的系统。
 
 对于极大量样本 （>10000） 联合队列分型，内存可能会成为一个问题，某些操作系统限制也会产生影响。在这种情况下，建议执行以下操作：

-  将操作系统打开文件限制设置为足够大的数量，以允许软件以打开足够的文件句柄。通过`ulimit -n NUM`命令完成。
 
 - 将操作系统堆栈限制设置为足够大的数量，以允许软件为操作分配足够的内存，因为内存与样本数，并且可能太大而无法容纳在默认情况下分配的堆栈中。 通过`ulimit -s NUM`命令完成。
 
 - 使用包含输入 GVCF 列表的文件，以防止命令长度超过操作系统的参数列表长度限制。您可以使用以下命令执行此操作：

```
sentieon driver ... --algo GVCFtyper output.vcf.gz - < input_files.txt
```

- 按照 jemalloc 应用说明中所述使用 jemalloc 内存分配，并通过以下方式更新 vcf 缓存大小，将以下代码添加到脚本中：

```
export VCFCACHE_BLOCKSIZE=4096
```

- 将分片大小设置为较小的数字，即50M。此外，使用 GVCFtyper 命令中的驱动程序选项`--traverse_param`，以确保所有线程都得到充分利用。该命令将变为：

```
sentieon driver -r $FASTA -t $NT --shard $shard1 --shard $shard2 \
 --traverse_param 10000/200 \
 --algo GVCFtyper GVCFtyper-shard_1.vcf.gz - < input_files.txt
```

## 1. 大型输出VCF文件的挑战

当运行超大规模队列时，输出 VCF 文件将包含海量的样本基因型信息列。这会导致输出文件变得难以处理，例如在完整文件上运行 VQSR 会变得不切实际，因此您可能需要寻求 VCF 存储以外的替代方案。

根据您存储输出的方法，按样本组或基因组坐标拆分输出可以缓解超大输出文件的问题；有关如何存储联合分型输出的具体需求，请随时联系 Sentieon 技术支持。
 
**（1）按基因组区域分割输出**

为了减小输出 VCF 文件的大小，您可以在特定的基因组子区域（例如单个染色体）上执行片段合并。您可以通过仅合并子集的中间VCF文件来实现此目的。例如，您可以通过仅合并涵盖每个染色体的片段来创建每个染色体的VCF文件：如果片段1-4涵盖chr1，而片段5同时涵盖chr1和chr2，则以下代码将创建一个仅包含chr1变异体的VCF文件：

```

#merge the necessary intermediate sharded VCFs
sentieon driver --passthru -t $nt --algo GVCFtyper --merge \
  GVCFtyper_chr1_tmp.vcf.gz \
  GVCFtyper-shard_1.vcf.gz GVCFtyper-shard_2.vcf.gz GVCFtyper-shard_3.vcf.gz \
  GVCFtyper-shard_4.vcf.gz GVCFtyper-shard_5.vcf.gz
#remove variants not in the proper contig and index the VCF
bcftools view GVCFtyper_chr1_tmp.vcf.gz  -r chr1 \
  -o - | sentieon util vcfconvert - GVCFtyper_chr1.vcf.gz
```

使用上述方法可以创建合法的 VCF 文件，其大小仅为完整 VCF 的一小部分。
为了在上述输出上运行 VQSR，您需要为 VarCal 算法提供一个包含全基因组所有变异记录的 VCF 文件，因为这是正确校准 VQSR 模型所必需的。然而，VarCal 仅需要 VCF 数据的前 8 列，因此您无需拼接所有基因组子区域的完整 VCF，只需提取并拼接每个文件的关键信息（前 8 列）来创建一个较小的 VCF 文件即可。使用以下代码可以创建所需文件：

```
vcf_list=(GVCFtyper_chr1.vcf.gz GVCFtyper_chr2.vcf.gz) # include more VCFs if applicable

#extract the first 8 columns from each region to a new compressed VCF
mkfifo tmp.fifo
sentieon util vcfconvert tmp.fifo GVCFtyper_annotations.vcf.gz &
 convert_pid=$!
 exec 3>tmp.fifo #a file descriptor to hold the fifo open
 bcftools view -h ${vcf_list[0]} | grep "^##" > tmp.fifo
 bcftools view -h ${vcf_list[0]} | grep "^#CHROM" | cut -f 1-8 > tmp.fifo
for vcf in ${vcf_list[@]}; do
  bcftools view -H $vcf | cut -f 1-8 > tmp.fifo
done
exec 3>&-
wait $convert_pid
rm tmp.fifo
```

随后，可以在包含前 8 列的合并 VCF 上运行 VarCal 算法，为 SNP 和 INDEL 生成 tranches 和 recal 文件。之后可以使用 ApplyVarCal 算法直接对按基因组区域拆分的 VCF 应用 VQSR：

```
vcf_list=(GVCFtyper_chr1.vcf.gz GVCFtyper_chr2.vcf.gz) # include more VCFs if applicable

#apply variant quality score recalibration
for vcf in ${vcf_list[@]}; do
  out_vcf=${vcf%%.vcf.gz}_snp-indel-recal.vcf.gz
  sentieon driver -r $FASTA -t $nt --algo ApplyVarCal -v $vcf \
    --vqsr_model var_type=SNP,tranches_file=${snp_tranches},sensitivity=99.0,recal=${snp_recal} \
    --vqsr_model var_type=INDEL,tranches_file=${indel_tranches},sensitivity=99.0,recal=${indel_recal} \
    $out_vcf
done
```

**（2）按样本组分割输出**

或者，Sentieon® GVCFtyper 合并步骤包含` --split_by_sample` 算法选项作为应对这一挑战的潜在方案。`--split_by_sample` 选项在合并步骤中使用，可生成合法的局部 VCF，每个文件包含样本的子集，从而将完整的 VCF 文件分隔为更小、更易于管理的输出文件。其用法如下：

```
sentieon driver --passthru -t $nt --algo GVCFtyper --merge \
  --split_by_sample split.conf GVCFtyper_main.vcf.gz \
  GVCFtyper-shard_1.vcf.gz GVCFtyper-shard_2.vcf.gz …
```

`--split_by_sample`选项使用的 `split.conf` 输入文件是一个以制表符分隔的文件，每行以特定样本组的输出文件开头，后跟该组对应样本的制表符分隔列表。您可以使用多行对应同一个输出文件，这会将多行中的所有样本归为一组。以下示例展示了两个样本组，其中样本 1 到 5 将输出到 group1 文件，而样本 6 到 8 将输出到 group2 文件：

```
GVCFtyper_file_group1.vcf.gz Sample1 Sample2 Sample3
GVCFtyper_file_group1.vcf.gz Sample4 Sample5
GVCFtyper_file_group2.vcf.gz Sample6 Sample7 Sample8
```

上述GVCFtyper合并命令将生成以下输出：

- `GVCFtyper_main.vcf.gz`：一个包含所有VCF信息，包括INFO列在内的部分VCF文件。不包含样本信息。此文件对于运行VQSR非常有用，因为它包含了变异校准所需的所有必要信息。
 
 - `GVCFtyper_file_group1.vcf.gz`：一个包含至INFO列、FORMAT列以及group1中所有样本的所有列的部分VCF文件。
 
 - `GVCFtyper_file_group2.vcf.gz`：一个包含至INFO列、FORMAT列以及group2中所有样本的所有列的部分VCF文件。

然后，您可以使用bcftools合并部分VCF并选择您感兴趣的样本。您可以使用以下代码来实现：

```
bash extract.sh GVCFtyper_main.vcf.gz \
  split.conf Samle1,Sample4,Sample7
```

在该脚本中，第三个参数是一个逗号分隔的感兴趣样本列表，而`extract.sh`脚本如下所示：

```
#!/bin/bash
MAIN=$1
GRPS_CONF=$2
SAMPLES=$3
REGIONS=$4
 BCF=/home/release/other_tools/bcftools-1.3/bcftools
if [ $# -eq 4 ] || [ $# -eq 3 ]; then
  if [ $REGIONS == "" ]; then
    BED_ARG=""
  else
    BED_ARG="-R $REGIONS"
  fi
  #parse input files from group config file
  GRPS=$(grep -v "^#" $GRPS_CONF| cut -f1 | sort | uniq)
  $BCF view -h $MAIN | grep -v '^#CHROM' | grep -v '^##bcftools' > out.vcf
  hdr=$($BCF view -h $MAIN | grep '^#CHROM')
  hdr="$hdr\tFORMAT"
  arg="<($BCF view $BED_ARG -H -I $MAIN | cut -f -8)"
  col=9
  for g in $GRPS; do
    s=$($BCF view --force-samples -s $SAMPLES -h $g 2>/dev/null | grep '^#CHROM' | \
       cut -f 10-)
    [ -z "$s" ] && continue
    hdr="$hdr\t$s"
    c="<($BCF view --force-samples -s $SAMPLES $BED_ARG -H -I $g 2>/dev/null | \
      cut -f $col-)"
    arg="$arg $c"
    col=10
  done
echo -e "$hdr" >> out.vcf
eval paste "$arg" >> out.vcf
else
  echo "usage $0 main_vcf_file group_config_file csv_sample_list [bed_file]"
fi
```
此处提供的方案允许在 `GVCFtyper_main.vcf.gz` 上轻松运行 VQSR 以创建重校准后的文件，随后您可以在 extract.sh 脚本中使用该文件以获得 VQSR 后的结果。

## 2. 云环境的特殊考虑事项

在云环境中，通常没有可以容纳大量 GVCF 输入的 NFS 存储。对于联合分型，Sentieon® GVCFtyper 允许使用托管在对象存储位置（如 AWS S3）或通过 HTTP 访问位置的 GVCF 输入文件。然而，对于大规模队列（100+），使用对象存储输入的 GVCFtyper 命令对内存的需求可能过高，因此不建议采用这种访问输入的方法。

云环境中的推荐方法是根据节点处理的分片，将部分 GVCF 输入文件下载到计算节点。GVCF 输入的部分下载可以使用 bcftools 完成，但务必在 bcftools 命令中添加 `--no-version` 选项，以确保不同分片的表头不会因差异过大而导致 GVCFtyper 合并步骤拒绝它们：

```

#use bcftools to download shard1 GVCF partial inputs and create the index
bcftools view --no-version -r $shard1_csv_intervals -o - $URL_sample.g.vcf.gz | \
  sentieon util vcfconvert - ${sample}_s1.g.vcf.gz
#or download multiple files in parallel using xargs using a list
cat list_inputs.txt | xargs -P $NUM_PARALLEL_DOWNLOAD -I @ | \
  sh -c "bcftools view -r $shard1_csv_intervals -o - $S3BUCKET_INPUTS/@ | \
  sentieon util vcfconvert - @"
```

为了获得额外的效率，可以使用“瀑布式”方法，如图3和图4所示。在这种方法中，计算节点将按顺序处理多个分片，同时并行下载下一个分片的部分 GVCF 输入与当前分片的 GVCFtyper 处理，从而通过 CPU 密集型和 I/O 密集型进程的结合更有效地共享资源。流程如下：

-  下载shard1的部分GVCF输入文件。
 
-  启动shard1的GVCFtyper，并与此同时并行下载shard2的部分GVCF输入文件。
 
 - 在上一步完成之后，启动shard2的GVCFtyper，并与此同时使用bcftools下载shard3。
 
![](https://img2024.cnblogs.com/blog/3218664/202407/3218664-20240715174422138-40386455.png#align=center)

<center>图3 在云环境中将GVCFtyper分发到多个服务器</center>

![](https://img2024.cnblogs.com/blog/3218664/202407/3218664-20240715174433536-2107229297.png#align=center)

<center>图4 在云环境中将GVCFtyper分发到多个服务器，下载下一个分片的输入并与上一个分片的处理并行进行的瀑布式方法的详细说明</center>

---

# 七、Shell示例

以下是用于分发命令的Shell示例，使用1G碱基的分片大小，仅选用于演示目的。实际运行中推荐的分片大小是100M碱基。

```
# Sample file for distributing DNAseq pipeline onto 4 1GBase shards
# Each stage command can be distributed to a different server for faster processing,
# but the user needs to make sure that the necessary files are present in each machine

FASTA="/home/b37/human_g1k_v37_decoy.fasta"
FASTQ_1="WGS-30x_1.fastq.gz"
FASTQ_2="WGS-30x_2.fastq.gz"
FASTQ_INDEX="WGS-30x.fastq.gz.index"
KNOWN1="/home/b37/1000G_phase1.indels.b37.vcf.gz"
KNOWN2="/home/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz"
DBSNP="/home/b37/dbsnp_138.b37.vcf.gz"

#######################################
# BWA mapping, distributed on 4 servers
#######################################
BWA_K_size=100000000
num_srvr=4

#get the number of runs that the inputs will be split into given the size
num_K=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f1)
BWA_K_size=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f2)
num_K_srvr=$(expr $num_K / $num_srvr + 1)
#run multiple BWA on multiple servers
 file_list=""
for run in $(seq 0 $((num_srvr-1))); do
        region="$((run*num_K_srvr))-$((run*num_K_srvr+num_K_srvr))"
         #send each of these to a different server
        sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
          -K $BWA_K_size -t $NT -p $FASTA \
          "< sentieon fqidx extract -i $FASTQ_INDEX -r $region $FASTQ_1 $FASTQ_2" | \
          sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
       file_list="$file_list sorted_run$run.bam"
done

#merge the results
 sentieon util merge -o sorted.bam $file_list

#######################################
# define 4 shards
#######################################
SHARD_0="--shard 1:1-249250621 --shard 2:1-243199373 --shard 3:1-198022430 \
  --shard 4:1-191154276 --shard 5:1-118373300"
SHARD_1="--shard 5:118373301-180915260 --shard 6:1-171115067 --shard 7:1-159138663 \
  --shard 8:1-146364022 --shard 9:1-141213431 --shard 10:1-135534747 \
  --shard 11:1-135006516 --shard 12:1-49085594"
SHARD_2="--shard 12:49085595-133851895 --shard 13:1-115169878 --shard 14:1-107349540 \
   --shard 15:1-102531392 --shard 16:1-90354753 --shard 17:1-81195210 \
   --shard 18:1-78077248 --shard 19:1-59128983 --shard 20:1-63025520 \
   --shard 21:1-48129895 --shard 22:1-51304566 --shard X:1-118966714"
SHARD_3="--shard X:118966715-155270560 --shard Y:1-59373566 --shard MT:1-16569 \
   --shard GL000207.1:1-4262 --shard GL000226.1:1-15008 --shard GL000229.1:1-19913 \
   --shard GL000231.1:1-27386 --shard GL000210.1:1-27682 --shard GL000239.1:1-33824 \
   --shard GL000235.1:1-34474 --shard GL000201.1:1-36148 --shard GL000247.1:1-36422 \
   --shard GL000245.1:1-36651 --shard GL000197.1:1-37175 --shard GL000203.1:1-37498 \
   --shard GL000246.1:1-38154 --shard GL000249.1:1-38502 --shard GL000196.1:1-38914 \
   --shard GL000248.1:1-39786 --shard GL000244.1:1-39929 --shard GL000238.1:1-39939 \
   --shard GL000202.1:1-40103 --shard GL000234.1:1-40531 --shard GL000232.1:1-40652 \
   --shard GL000206.1:1-41001 --shard GL000240.1:1-41933 --shard GL000236.1:1-41934 \
   --shard GL000241.1:1-42152 --shard GL000243.1:1-43341 --shard GL000242.1:1-43523 \
   --shard GL000230.1:1-43691 --shard GL000237.1:1-45867 --shard GL000233.1:1-45941 \
   --shard GL000204.1:1-81310 --shard GL000198.1:1-90085 --shard GL000208.1:1-92689 \
   --shard GL000191.1:1-106433 --shard GL000227.1:1-128374 \
   --shard GL000228.1:1-129120 --shard GL000214.1:1-137718 \
   --shard GL000221.1:1-155397 --shard GL000209.1:1-159169 \
   --shard GL000218.1:1-161147 --shard GL000220.1:1-161802 \
   --shard GL000213.1:1-164239 --shard GL000211.1:1-166566 \
   --shard GL000199.1:1-169874 --shard GL000217.1:1-172149 \
   --shard GL000216.1:1-172294 --shard GL000215.1:1-172545 \
   --shard GL000205.1:1-174588 --shard GL000219.1:1-179198 \
   --shard GL000224.1:1-179693 --shard GL000223.1:1-180455 \
   --shard GL000195.1:1-182896 --shard GL000212.1:1-186858 \
   --shard GL000222.1:1-186861 --shard GL000200.1:1-187035 \
   --shard GL000193.1:1-189789 --shard GL000194.1:1-191469 \
   --shard GL000225.1:1-211173 --shard GL000192.1:1-547496 \
   --shard NC_007605:1-171823 --shard hs37d5:1-35477943"

#######################################
# Locus collector
#######################################
#Locus Collector for shard 000
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_0 \
   --algo LocusCollector --fun score_info .part_deduped.bam000.score.gz \
   2> collect000.log
#Locus Collector for shard 001
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_1 \
   --algo LocusCollector --fun score_info .part_deduped.bam001.score.gz \
   2> collect001.log
#Locus Collector for shard 002
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_2 \
   --algo LocusCollector --fun score_info .part_deduped.bam002.score.gz \
   2> collect002.log
#Locus Collector for shard 003
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_3 \
   --algo LocusCollector --fun score_info .part_deduped.bam003.score.gz \
   2> collect003.log
#Locus Collector for shard with unmapped reads (NO_COOR)
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam --shard NO_COOR \
   --algo LocusCollector --fun score_info .part_deduped.bam004.score.gz \
    2> collect004.log

#######################################
# Dedup using all score.gz files
#######################################
#Dedup for shard 000
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_0 \
   --algo Dedup --score_info .part_deduped.bam000.score.gz \
   --score_info .part_deduped.bam001.score.gz \
   --score_info .part_deduped.bam002.score.gz \
   --score_info .part_deduped.bam003.score.gz \
   --score_info .part_deduped.bam004.score.gz \
   --rmdup .part_deduped000.bam 2> dedup000.log
#Dedup for shard 001
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_1 \
   --algo Dedup --score_info .part_deduped.bam000.score.gz \
   --score_info .part_deduped.bam001.score.gz \
   --score_info .part_deduped.bam002.score.gz \
   --score_info .part_deduped.bam003.score.gz \
   --score_info .part_deduped.bam004.score.gz \
   --rmdup .part_deduped001.bam 2> dedup001.log
#Dedup for shard 002
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_2 \
   --algo Dedup --score_info .part_deduped.bam000.score.gz \
   --score_info .part_deduped.bam001.score.gz \
   --score_info .part_deduped.bam002.score.gz \
   --score_info .part_deduped.bam003.score.gz \
   --score_info .part_deduped.bam004.score.gz \
   --rmdup .part_deduped002.bam 2> dedup002.log
#Dedup for shard 003
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_3 \
   --algo Dedup --score_info .part_deduped.bam000.score.gz \
   --score_info .part_deduped.bam001.score.gz \
   --score_info .part_deduped.bam002.score.gz \
   --score_info .part_deduped.bam003.score.gz \
   --score_info .part_deduped.bam004.score.gz \
   --rmdup .part_deduped003.bam 2> dedup003.log
#Dedup for shard with unmapped reads (NO_COOR)
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam --shard NO_COOR \
   --algo Dedup --score_info .part_deduped.bam000.score.gz \
   --score_info .part_deduped.bam001.score.gz \
   --score_info .part_deduped.bam002.score.gz \
   --score_info .part_deduped.bam003.score.gz \
   --score_info .part_deduped.bam004.score.gz \
   --rmdup .part_deduped004.bam 2> dedup004.log
#Merge bam files from all shards into final output
$SENTIEON_FOLDER/bin/sentieon util merge -i .part_deduped000.bam \
   -i .part_deduped001.bam -i .part_deduped002.bam \
   -i .part_deduped003.bam -i .part_deduped004.bam \
   -o deduped.bam --mergemode=10 2> dedup_merge.log

#######################################
# BQSR
#######################################
#QualCal for shard 000
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
   -i .part_deduped001.bam $SHARD_0 --algo QualCal -k $KNOWN1 -k $KNOWN2 \
   .part_recal_data000.table 2> bqsr000.log
#QualCal for shard 001
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
   -i .part_deduped001.bam -i .part_deduped002.bam $SHARD_1 --algo QualCal \
   -k $KNOWN1 -k $KNOWN2 .part_recal_data001.table 2> bqsr001.log
#QualCal for shard 002
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped001.bam \
   -i .part_deduped002.bam -i .part_deduped003.bam $SHARD_2 --algo QualCal \
   -k $KNOWN1 -k $KNOWN2 .part_recal_data002.table 2> bqsr002.log
#QualCal for shard 003
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped002.bam \
   -i .part_deduped003.bam $SHARD_3 --algo QualCal -k $KNOWN1 -k $KNOWN2 \
   .part_recal_data003.table 2> bqsr003.log
#QualCal for shard with unmapped reads (NO_COOR)
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped004.bam \
   --shard NO_COOR --algo QualCal -k $KNOWN1 -k $KNOWN2 \
   .part_recal_data004.table 2> bqsr004.log
#Merge table files into complete calibration table
$SENTIEON_FOLDER/bin/sentieon driver --passthru --algo QualCal \
   --merge recal_data.table .part_recal_data000.table .part_recal_data001.table \
   .part_recal_data002.table .part_recal_data003.table .part_recal_data004.table \
   2> bqsr_merge.log

#######################################
# Variant Calling
#######################################
#Haplotyper for shard 000
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
   -i .part_deduped001.bam -q recal_data.table $SHARD_0 --algo Haplotyper \
   -d $DBSNP .part_output000.vcf.gz 2> hc000.log
#Haplotyper for shard 001
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
   -i .part_deduped001.bam -i .part_deduped002.bam  -q recal_data.table \
   $SHARD_1 --algo Haplotyper -d $DBSNP .part_output001.vcf.gz 2> hc001.log
#Haplotyper for shard 002
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped001.bam \
   -i .part_deduped002.bam -i .part_deduped003.bam -q recal_data.table \
   $SHARD_2 --algo Haplotyper -d $DBSNP .part_output002.vcf.gz 2> hc002.log
#Haplotyper for shard 003
$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped002.bam \
   -i .part_deduped003.bam -q recal_data.table $SHARD_3 --algo Haplotyper \
  -d $DBSNP .part_output003.vcf.gz 2> hc003.log
#There is no need to do variant calling on the NO_COOR shard, as there will
# be no variants on unmapped reads#Merge output files into output VCF
$SENTIEON_FOLDER/bin/sentieon driver --passthru --algo Haplotyper \
   --merge output.vcf.gz .part_output000.vcf.gz .part_output001.vcf.gz \
   .part_output002.vcf.gz .part_output003.vcf.gz 2> hc_merge.log
```
以下是一个用于创建分片的Shell示例。

```
determine_shards_from_bam(){
    local bam step tag chr len pos end
    bam="$1"
    step="$2"
    pos=1
    samtools view -H $bam |\
    while read tag chr len; do
        [ $tag == '@SQ' ] || continue
        chr=$(expr "$chr" : 'SN:$.*$')
        len=$(expr "$len" : 'LN:$.*$')
        while [ $pos -le $len ]; do
            end=$(($pos + $step - 1))
            if [ $pos -lt 0 ]; then
                start=1
            else
                start=$pos
            fi
            if [ $end -gt $len ]; then
               echo -n "$chr:$start-$len "
                pos=$(($pos-$len))
                break
            else
                echo "$chr:$start-$end"
                pos=$(($end + 1))
            fi
        done
    done
    echo "NO_COOR"
}

determine_shards_from_dict(){
    local bam step tag chr len pos end
    dict="$1"
    step="$2"
    pos=1
    cat $dict |\
    while read tag chr len UR; do
        [ $tag == '@SQ' ] || continue
        chr=$(expr "$chr" : 'SN:$.*$')
        len=$(expr "$len" : 'LN:$.*$')
        while [ $pos -le $len ]; do
           end=$(($pos + $step - 1))
           if [ $pos -lt 0 ]; then
                start=1
            else
                start=$pos
            fi
            if [ $end -gt $len ];then
                echo -n "$chr:$start-$len "
                pos=$(($pos-$len))
                break
            else
                echo "$chr:$start-$end"
                pos=$(($end + 1))
            fi
        done
    done
    echo "NO_COOR"
}

determine_shards_from_fai(){
    local bam step tag chr len pos end
    fai="$1"
    step="$2"
    pos=1
    cat $fai |\
   while read chr len other; do
        while [ $pos -le $len ]; do
            end=$(($pos + $step - 1))
           if [ $pos -lt 0 ]; then
                start=1
            else
                start=$pos
            fi
            if [ $end -gt $len ]; then
                echo -n "$chr:$start-$len "
               pos=$(($pos-$len))
                break
            else
                echo "$chr:$start-$end"
                pos=$(($end + 1))
            fi
        done
    done
   echo "NO_COOR"
}

if [ $# -eq 2 ]; then
    filename=$(basename "$1")
    extension="${filename##*.}"
    if [ "$extension" = "fai" ]; then
        determine_shards_from_fai $1 $2
    elif [ "$extension" = "bam" ]; then
        determine_shards_from_bam $1 $2
    elif [ "$extension" = "dict" ]; then
        determine_shards_from_dict $1 $2
    fi
else
    echo "usage $0 file shard_size"
fi
```

[**想了解更多Sentieon软件应用教程，可以点击此处进行跳转**](https://doc.insvast.com/doc/10/)