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

原创

INSVAST

发布于 2024-07-15 16:19:31

510

发布于 2024-07-15 16:19:31

文章被收录于专栏：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所示，说明了这样一个流程的数据流。

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

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

以图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阶段需要文件，以及可用的文件和文件。该阶段生成一个文件。i-th part_deduped$shard_i.bam part_deduped$shard_i+1.bam part_deduped$shard_i-1.bam part_recal_data$shard_i.table
在QualCal之后，所有文件需要合并为一个单一的校准表文件，用于变异调用。驱动程序和QualCal的选项可以用于执行边界感知合并。part_recal_data$shard_i.table --passthru --merge
分片的Haplotyper阶段需要文件，以及可用的文件和文件；此外，完整合并的校准表需要作为输入。该阶段生成一个文件。i-th part_deduped$shard_i.bam part_deduped$shard_i+1.bam part_deduped$shard_i-1.bam part_output$shard_i.vcf.gz
在Haplotyper之后，所有文件需要合并为一个单一的输出VCF文件。驱动程序和Haplotyper的选项用于执行边界感知合并。part_output$shard_i.vcf.gz --passthru --merge
如果在流程的任何阶段需要合并的输出bam文件，可以使用util二进制文件执行边界感知合并；需要在命令中添加选项，以便util merge不处理读取，而只是按块复制它们。--mergemode=10

重要提示：在合并结果时，务必按照正确的顺序依次输入部分结果，否则结果将不正确。

分发BWA

上述说明中未包含有关使用BWA进行分发对齐的任何信息。为了分发BWA对齐，您可以使用Sentieon®工具中提供的工具为输入的FASTQ文件创建索引文件；然后，您可以使用fqidx命令的结果作为BWA mem的输入，在不同的服务器上处理FASTQ文件的特定部分；您需要确保在BWA命令中包含该选项，因为fqidx的输出包含交错的reads在单个输出中。这个流程的结果与在单次运行上执行的结果是相同的。fqidx -p

 1BWA_K_size=100000000
 2num_splits=10
 3SAMPLE="sample_name" #Sample name SM tag in bam
 4GROUP="read_group_name" #Read Group name RGID tag in bam
 5PLATFORM="ILLUMINA"
 6NT=$(nproc) #number of threads to use in computation, set to all threads available
 7
 8FASTA="/home/regression/references/b37/human_g1k_v37_decoy.fasta"
 9FASTQ_1="WGS-30x_1.fastq.gz"
10FASTQ_2="WGS-30x_2.fastq.gz"
11FASTQ_INDEX="WGS-30x.fastq.gz.index"
12
13#create FASTQ indices
14sentieon fqidx index -K $BWA_K_size -o $FASTQ_INDEX $FASTQ_1 $FASTQ_2
15
16#get the number of runs that the inputs will be split into given the size
17num_K=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f1)
18BWA_K_size=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f2)
19num_K_splits=$(expr $num_K / $num_splits + 1)
20#run multiple BWA on multiple servers
21file_list=""
22for run in $(seq 0 $((num_splits-1))); do
23        region="$((run*num_K_splits))-$((run*num_K_splits+num_K_splits))"
24        #send each of these to a different server
25        sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
26          -K $BWA_K_size -t $NT -p $FASTA \
27          "<sentieon fqidx extract -i $FASTQ_INDEX -r $region $FASTQ_1 $FASTQ_2" | \
28          sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
29        file_list="$file_list sorted_run$run.bam"
30done
31
32#merge the results
33sentieon util merge -o sorted.bam $file_list
34#or perform deduplication on all the intermediate BAM files
35file_list_bam=""
36for file in $file_list; do file_list_bam="$file_list_bam -i $file"; done
37sentieon driver -r $FASTA -t $NT $file_list --algo LocusCollector --fun score_info score.txt.gz
38sentieon driver -r $FASTA -t $NT $file_list --algo Dedup --score_info score.txt.gz deduped.bam

1.当无法创建FASTQ索引文件（版本201808.02+）时进行BWA分发 如果无法创建FASTQ索引文件，则可以使用带有分数选项的util fqidx工具，并同时使用该选项。此方法将输入的FASTQ文件分割成多个读取块的片段，并从每个片段中提取每个元素以在不同的服务器上进行处理，从起始位置到结束位置。-F m/n -K n m-th m 0 n-1 请注意，此功能仅建议在文件存储速度足够快以支持每个fqidx进程同时读取输入FASTQ文件的IT环境中使用。这在云环境中很常见，或者在具有高带宽NFS系统的本地集群环境中使用。

 1BWA_K_size=100000000
 2num_splits=10
 3SAMPLE="sample_name" #Sample name SM tag in bam
 4GROUP="read_group_name" #Read Group name RGID tag in bam
 5PLATFORM="ILLUMINA"
 6NT=$(nproc) #number of threads to use in computation, set to all threads available
 7
 8FASTA="/home/regression/references/b37/human_g1k_v37_decoy.fasta"
 9FASTQ_1="WGS-30x_1.fastq.gz"
10FASTQ_2="WGS-30x_2.fastq.gz"
11
12#run multiple BWA on multiple servers
13file_list=""
14for run in $(seq 0 $((num_splits-1))); do
15        #send each of these to a different server
16        sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
17          -K $BWA_K_size -t $NT -p $FASTA \
18          "<sentieon fqidx extract -F $((run))/$num_splits -K $BWA_K_size $FASTQ_1 $FASTQ_2" | \
19          sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
20        file_list="$file_list sorted_run$run.bam"
21done
22
23#merge the results
24sentieon util merge -o sorted.bam $file_list
25#or perform deduplication on all the intermediate BAM files
26file_list_bam=""
27for file in $file_list; do file_list_bam="$file_list_bam -i $file"; done
28sentieon driver -r $FASTA -t $NT $file_list --algo LocusCollector --fun score_info score.txt.gz
29sentieon 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协议进行访问的远程对象存储位置。

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

重要提示： 在处理每个片段时，输入GVCF的列表需要保持一致的顺序，因为最终输出中的样本顺序将取决于输入顺序，并且合并需要所有部分文件具有相同的样本顺序。使用--shard选项时，GVCFtyper生成的输出VCF文件不是有效的VCF文件，因此在下面描述的合并之前不应使用它们。在处理完所有片段后，您需要运行GVCFtyper的合并命令来合并结果，确保中间VCF输入的顺序与参考基因组一致。输入文件需要在共享位置（例如NFS存储或可通过s3 http/https协议访问的远程对象存储位置）中可用。

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

总体运行时间和资源要求

为了减少整体运行时间，建议使用足够的分片，允许将运行时细粒度分布为多个服务器。分片可以按分片大小或预期创建复杂性。但是，单个分片的最终合并步骤结果无法分发，需要在单个服务器中运行; 此事实设置了用于分布，因为合并可以主导整个运行时。 GVCFtyper 算法是一种非常有效的算法，可能会 I/O 受 GVCF 存储位置性能的限制。因此，建议将 GVCF 输入存储在文件中提供 600 MBps 传输速率的系统。对于极大量样本（>10000）联合队列调用，内存可能会成为一个问题，某些操作系统限制也可能发挥作用。在这种情况下，建议执行以下操作：

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

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

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

1export VCFCACHE_BLOCKSIZE=4096

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

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

大型输出VCF文件的挑战

当运行非常大的队列时，输出VCF文件将包含大量列，其中包含每个样本的基因型信息。这么多列可能使输出文件变得难以处理，例如在完整文件上运行VQSR可能不切实际，因此您可能需要考虑替代VCF存储输出的方法。根据您用于存储输出的方法，将输出按样本组或基因组坐标分割可能会改善大型输出文件的问题；请随时与Sentieon支持团队联系，告知您如何存储联合调用的输出的具体要求。 （1）按基因组区域分割输出 为了使输出VCF文件变小，您可以在特定的基因组子区域（例如单个染色体）上执行片段合并。您可以通过仅合并子集的中间VCF文件来实现此目的。例如，您可以通过仅合并涵盖每个染色体的片段来创建每个染色体的VCF文件：如果片段1-4涵盖chr1，而片段5同时涵盖chr1和chr2，则以下代码将创建一个仅包含chr1变异体的VCF文件：

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

使用上述方法，您可以创建有效的VCF文件，其大小仅为完整VCF文件的一小部分。为了在上述输出上运行VQSR，您需要提供一个包含完整基因组范围内所有变异记录的VCF文件，因为这些记录需要用于对VQSR模型进行适当的校准。然而，VQSR仅需要VCF数据的前8列，因此您无需将每个特定基因组子区域的所有VCF文件连接起来，可以通过提取和连接每个文件的前8列来创建一个包含必要信息的较小的VCF文件。使用下面的代码将创建所需的文件：VarCal VarCal

 1vcf_list=(GVCFtyper_chr1.vcf.gz GVCFtyper_chr2.vcf.gz) # include more VCFs if applicable
 2
 3#extract the first 8 columns from each region to a new compressed VCF
 4mkfifo tmp.fifo
 5sentieon util vcfconvert tmp.fifo GVCFtyper_annotations.vcf.gz &
 6convert_pid=$!
 7exec 3>tmp.fifo #a file descriptor to hold the fifo open
 8bcftools view -h ${vcf_list[0]} | grep "^##" > tmp.fifo
 9bcftools view -h ${vcf_list[0]} | grep "^#CHROM" | cut -f 1-8 > tmp.fifo
10for vcf in ${vcf_list[@]}; do
11  bcftools view -H $vcf | cut -f 1-8 > tmp.fifo
12done
13exec 3>&-
14wait $convert_pid
15rm tmp.fifo

接下来，可以在包含VCF的前八列的合并VCF上运行算法，生成用于SNP和INDEL的分段和重新校准文件。然后，可以直接将VQSR应用于按基因组区域分割的VCF，使用算法进行操作：VarCal ApplyVarCal

 1vcf_list=(GVCFtyper_chr1.vcf.gz GVCFtyper_chr2.vcf.gz) # include more VCFs if applicable
 2
 3#apply variant quality score recalibration
 4for vcf in ${vcf_list[@]}; do
 5  out_vcf=${vcf%%.vcf.gz}_snp-indel-recal.vcf.gz
 6  sentieon driver -r $FASTA -t $nt --algo ApplyVarCal -v $vcf \
 7    --vqsr_model var_type=SNP,tranches_file=${snp_tranches},sensitivity=99.0,recal=${snp_recal} \
 8    --vqsr_model var_type=INDEL,tranches_file=${indel_tranches},sensitivity=99.0,recal=${indel_recal} \
 9    $out_vcf
10done

（2）按样本组分割输出 另外，Sentieon® GVCFtyper合并工具提供了算法选项作为解决这一挑战的潜在方案。在合并步骤中使用算法选项可以生成有效的部分VCF文件，每个文件包含一部分样本，从而将完整的VCF文件分割成较小、更易处理的输出文件。使用算法选项的用法为：--split_by_sample--split_by_sample--split_by_sample

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

算法选项中使用的输入文件是一个以制表符分隔的文件，每行以特定样本组的输出文件开头，后跟以制表符分隔的相应样本组样本列表。您可以使用多行具有相同输出文件的方式，将多个行中的所有样本分组。下面的示例显示了两个样本组，其中样本1到5将输出到group1输出文件，而样本6到8将输出到group2输出文件：split.conf --split_by_sample

1GVCFtyper_file_group1.vcf.gz Sample1 Sample2 Sample3
2GVCFtyper_file_group1.vcf.gz Sample4 Sample5
3GVCFtyper_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并选择您感兴趣的样本。您可以使用以下代码来实现：

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

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

 1#!/bin/bash
 2MAIN=$1
 3GRPS_CONF=$2
 4SAMPLES=$3
 5REGIONS=$4
 6BCF=/home/release/other_tools/bcftools-1.3/bcftools
 7if [ $# -eq 4 ] || [ $# -eq 3 ]; then
 8  if [ $REGIONS == "" ]; then
 9    BED_ARG=""
10  else
11    BED_ARG="-R $REGIONS"
12  fi
13  #parse input files from group config file
14  GRPS=$(grep -v "^#" $GRPS_CONF| cut -f1 | sort | uniq)
15  $BCF view -h $MAIN | grep -v '^#CHROM' | grep -v '^##bcftools' > out.vcf
16  hdr=$($BCF view -h $MAIN | grep '^#CHROM')
17  hdr="$hdr\tFORMAT"
18  arg="<($BCF view $BED_ARG -H -I $MAIN | cut -f -8)"
19  col=9
20  for g in $GRPS; do
21    s=$($BCF view --force-samples -s $SAMPLES -h $g 2>/dev/null | grep '^#CHROM' | \
22       cut -f 10-)
23    [ -z "$s" ] && continue
24    hdr="$hdr\t$s"
25    c="<($BCF view --force-samples -s $SAMPLES $BED_ARG -H -I $g 2>/dev/null | \
26      cut -f $col-)"
27    arg="$arg $c"
28    col=10
29  done
30echo -e "$hdr" >> out.vcf
31eval paste "$arg" >> out.vcf
32else
33  echo "usage $0 main_vcf_file group_config_file csv_sample_list [bed_file]"
34fi

云环境的特殊考虑事项

在云环境中，通常没有能够容纳大量GVCF输入的NFS存储。对于联合基因分型，Sentieon® GVCFtyper允许托管在对象存储位置（如AWS S3）或通过HTTP访问的GVCF输入文件。然而，对于大型队列（100+），使用对象存储输入的GVCFtyper命令的内存需求可能会过高，因此不建议使用此方法访问输入文件。在云环境中推荐的方法是根据节点处理的分片，将部分GVCF输入文件下载到计算节点上。可以使用bcftools进行GVCF输入的部分下载，但是需要在bcftools命令中添加--no-version选项，以确保不同分片的标头不会因差异而导致GVCFtyper拒绝合并它们。

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

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

下载shard1的部分GVCF输入文件。
启动shard1的GVCFtyper，并与此同时并行下载shard2的部分GVCF输入文件。
在上一步完成之后，启动shard2的GVCFtyper，并与此同时使用bcftools下载shard3。

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

Shell示例

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

  1# Sample file for distributing DNAseq pipeline onto 4 1GBase shards
  2# Each stage command can be distributed to a different server for faster processing,
  3# but the user needs to make sure that the necessary files are present in each machine
  4
  5FASTA="/home/b37/human_g1k_v37_decoy.fasta"
  6FASTQ_1="WGS-30x_1.fastq.gz"
  7FASTQ_2="WGS-30x_2.fastq.gz"
  8FASTQ_INDEX="WGS-30x.fastq.gz.index"
  9KNOWN1="/home/b37/1000G_phase1.indels.b37.vcf.gz"
 10KNOWN2="/home/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz"
 11DBSNP="/home/b37/dbsnp_138.b37.vcf.gz"
 12
 13#######################################
 14# BWA mapping, distributed on 4 servers
 15#######################################
 16BWA_K_size=100000000
 17num_srvr=4
 18
 19#get the number of runs that the inputs will be split into given the size
 20num_K=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f1)
 21BWA_K_size=$(sentieon fqidx query -i $FASTQ_INDEX | cut -d' ' -f2)
 22num_K_srvr=$(expr $num_K / $num_srvr + 1)
 23#run multiple BWA on multiple servers
 24file_list=""
 25for run in $(seq 0 $((num_srvr-1))); do
 26        region="$((run*num_K_srvr))-$((run*num_K_srvr+num_K_srvr))"
 27        #send each of these to a different server
 28        sentieon bwa mem -R "@RG\tID:$GROUP\tSM:$SAMPLE\tPL:$PLATFORM" \
 29          -K $BWA_K_size -t $NT -p $FASTA \
 30          "< sentieon fqidx extract -i $FASTQ_INDEX -r $region $FASTQ_1 $FASTQ_2" | \
 31          sentieon util sort -r $FASTA -t $NT --sam2bam -o sorted_run$run.bam -i -
 32        file_list="$file_list sorted_run$run.bam"
 33done
 34
 35#merge the results
 36sentieon util merge -o sorted.bam $file_list
 37
 38#######################################
 39# define 4 shards
 40#######################################
 41SHARD_0="--shard 1:1-249250621 --shard 2:1-243199373 --shard 3:1-198022430 \
 42  --shard 4:1-191154276 --shard 5:1-118373300"
 43SHARD_1="--shard 5:118373301-180915260 --shard 6:1-171115067 --shard 7:1-159138663 \
 44  --shard 8:1-146364022 --shard 9:1-141213431 --shard 10:1-135534747 \
 45  --shard 11:1-135006516 --shard 12:1-49085594"
 46SHARD_2="--shard 12:49085595-133851895 --shard 13:1-115169878 --shard 14:1-107349540 \
 47   --shard 15:1-102531392 --shard 16:1-90354753 --shard 17:1-81195210 \
 48   --shard 18:1-78077248 --shard 19:1-59128983 --shard 20:1-63025520 \
 49   --shard 21:1-48129895 --shard 22:1-51304566 --shard X:1-118966714"
 50SHARD_3="--shard X:118966715-155270560 --shard Y:1-59373566 --shard MT:1-16569 \
 51   --shard GL000207.1:1-4262 --shard GL000226.1:1-15008 --shard GL000229.1:1-19913 \
 52   --shard GL000231.1:1-27386 --shard GL000210.1:1-27682 --shard GL000239.1:1-33824 \
 53   --shard GL000235.1:1-34474 --shard GL000201.1:1-36148 --shard GL000247.1:1-36422 \
 54   --shard GL000245.1:1-36651 --shard GL000197.1:1-37175 --shard GL000203.1:1-37498 \
 55   --shard GL000246.1:1-38154 --shard GL000249.1:1-38502 --shard GL000196.1:1-38914 \
 56   --shard GL000248.1:1-39786 --shard GL000244.1:1-39929 --shard GL000238.1:1-39939 \
 57   --shard GL000202.1:1-40103 --shard GL000234.1:1-40531 --shard GL000232.1:1-40652 \
 58   --shard GL000206.1:1-41001 --shard GL000240.1:1-41933 --shard GL000236.1:1-41934 \
 59   --shard GL000241.1:1-42152 --shard GL000243.1:1-43341 --shard GL000242.1:1-43523 \
 60   --shard GL000230.1:1-43691 --shard GL000237.1:1-45867 --shard GL000233.1:1-45941 \
 61   --shard GL000204.1:1-81310 --shard GL000198.1:1-90085 --shard GL000208.1:1-92689 \
 62   --shard GL000191.1:1-106433 --shard GL000227.1:1-128374 \
 63   --shard GL000228.1:1-129120 --shard GL000214.1:1-137718 \
 64   --shard GL000221.1:1-155397 --shard GL000209.1:1-159169 \
 65   --shard GL000218.1:1-161147 --shard GL000220.1:1-161802 \
 66   --shard GL000213.1:1-164239 --shard GL000211.1:1-166566 \
 67   --shard GL000199.1:1-169874 --shard GL000217.1:1-172149 \
 68   --shard GL000216.1:1-172294 --shard GL000215.1:1-172545 \
 69   --shard GL000205.1:1-174588 --shard GL000219.1:1-179198 \
 70   --shard GL000224.1:1-179693 --shard GL000223.1:1-180455 \
 71   --shard GL000195.1:1-182896 --shard GL000212.1:1-186858 \
 72   --shard GL000222.1:1-186861 --shard GL000200.1:1-187035 \
 73   --shard GL000193.1:1-189789 --shard GL000194.1:1-191469 \
 74   --shard GL000225.1:1-211173 --shard GL000192.1:1-547496 \
 75   --shard NC_007605:1-171823 --shard hs37d5:1-35477943"
 76
 77#######################################
 78# Locus collector
 79#######################################
 80#Locus Collector for shard 000
 81$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_0 \
 82   --algo LocusCollector --fun score_info .part_deduped.bam000.score.gz \
 83   2> collect000.log
 84#Locus Collector for shard 001
 85$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_1 \
 86   --algo LocusCollector --fun score_info .part_deduped.bam001.score.gz \
 87   2> collect001.log
 88#Locus Collector for shard 002
 89$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_2 \
 90   --algo LocusCollector --fun score_info .part_deduped.bam002.score.gz \
 91   2> collect002.log
 92#Locus Collector for shard 003
 93$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_3 \
 94   --algo LocusCollector --fun score_info .part_deduped.bam003.score.gz \
 95   2> collect003.log
 96#Locus Collector for shard with unmapped reads (NO_COOR)
 97$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam --shard NO_COOR \
 98   --algo LocusCollector --fun score_info .part_deduped.bam004.score.gz \
 99   2> collect004.log
100
101#######################################
102# Dedup using all score.gz files
103#######################################
104#Dedup for shard 000
105$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_0 \
106   --algo Dedup --score_info .part_deduped.bam000.score.gz \
107   --score_info .part_deduped.bam001.score.gz \
108   --score_info .part_deduped.bam002.score.gz \
109   --score_info .part_deduped.bam003.score.gz \
110   --score_info .part_deduped.bam004.score.gz \
111   --rmdup .part_deduped000.bam 2> dedup000.log
112#Dedup for shard 001
113$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_1 \
114   --algo Dedup --score_info .part_deduped.bam000.score.gz \
115   --score_info .part_deduped.bam001.score.gz \
116   --score_info .part_deduped.bam002.score.gz \
117   --score_info .part_deduped.bam003.score.gz \
118   --score_info .part_deduped.bam004.score.gz \
119   --rmdup .part_deduped001.bam 2> dedup001.log
120#Dedup for shard 002
121$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_2 \
122   --algo Dedup --score_info .part_deduped.bam000.score.gz \
123   --score_info .part_deduped.bam001.score.gz \
124   --score_info .part_deduped.bam002.score.gz \
125   --score_info .part_deduped.bam003.score.gz \
126   --score_info .part_deduped.bam004.score.gz \
127   --rmdup .part_deduped002.bam 2> dedup002.log
128#Dedup for shard 003
129$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam $SHARD_3 \
130   --algo Dedup --score_info .part_deduped.bam000.score.gz \
131   --score_info .part_deduped.bam001.score.gz \
132   --score_info .part_deduped.bam002.score.gz \
133   --score_info .part_deduped.bam003.score.gz \
134   --score_info .part_deduped.bam004.score.gz \
135   --rmdup .part_deduped003.bam 2> dedup003.log
136#Dedup for shard with unmapped reads (NO_COOR)
137$SENTIEON_FOLDER/bin/sentieon driver -t 16 -i sorted.bam --shard NO_COOR \
138   --algo Dedup --score_info .part_deduped.bam000.score.gz \
139   --score_info .part_deduped.bam001.score.gz \
140   --score_info .part_deduped.bam002.score.gz \
141   --score_info .part_deduped.bam003.score.gz \
142   --score_info .part_deduped.bam004.score.gz \
143   --rmdup .part_deduped004.bam 2> dedup004.log
144#Merge bam files from all shards into final output
145$SENTIEON_FOLDER/bin/sentieon util merge -i .part_deduped000.bam \
146   -i .part_deduped001.bam -i .part_deduped002.bam \
147   -i .part_deduped003.bam -i .part_deduped004.bam \
148   -o deduped.bam --mergemode=10 2> dedup_merge.log
149
150#######################################
151# BQSR
152#######################################
153#QualCal for shard 000
154$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
155   -i .part_deduped001.bam $SHARD_0 --algo QualCal -k $KNOWN1 -k $KNOWN2 \
156   .part_recal_data000.table 2> bqsr000.log
157#QualCal for shard 001
158$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
159   -i .part_deduped001.bam -i .part_deduped002.bam $SHARD_1 --algo QualCal \
160   -k $KNOWN1 -k $KNOWN2 .part_recal_data001.table 2> bqsr001.log
161#QualCal for shard 002
162$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped001.bam \
163   -i .part_deduped002.bam -i .part_deduped003.bam $SHARD_2 --algo QualCal \
164   -k $KNOWN1 -k $KNOWN2 .part_recal_data002.table 2> bqsr002.log
165#QualCal for shard 003
166$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped002.bam \
167   -i .part_deduped003.bam $SHARD_3 --algo QualCal -k $KNOWN1 -k $KNOWN2 \
168   .part_recal_data003.table 2> bqsr003.log
169#QualCal for shard with unmapped reads (NO_COOR)
170$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped004.bam \
171   --shard NO_COOR --algo QualCal -k $KNOWN1 -k $KNOWN2 \
172   .part_recal_data004.table 2> bqsr004.log
173#Merge table files into complete calibration table
174$SENTIEON_FOLDER/bin/sentieon driver --passthru --algo QualCal \
175   --merge recal_data.table .part_recal_data000.table .part_recal_data001.table \
176   .part_recal_data002.table .part_recal_data003.table .part_recal_data004.table \
177   2> bqsr_merge.log
178
179#######################################
180# Variant Calling
181#######################################
182#Haplotyper for shard 000
183$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
184   -i .part_deduped001.bam -q recal_data.table $SHARD_0 --algo Haplotyper \
185   -d $DBSNP .part_output000.vcf.gz 2> hc000.log
186#Haplotyper for shard 001
187$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped000.bam \
188   -i .part_deduped001.bam -i .part_deduped002.bam  -q recal_data.table \
189   $SHARD_1 --algo Haplotyper -d $DBSNP .part_output001.vcf.gz 2> hc001.log
190#Haplotyper for shard 002
191$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped001.bam \
192   -i .part_deduped002.bam -i .part_deduped003.bam -q recal_data.table \
193   $SHARD_2 --algo Haplotyper -d $DBSNP .part_output002.vcf.gz 2> hc002.log
194#Haplotyper for shard 003
195$SENTIEON_FOLDER/bin/sentieon driver -t 16 -r $FASTA -i .part_deduped002.bam \
196   -i .part_deduped003.bam -q recal_data.table $SHARD_3 --algo Haplotyper \
197   -d $DBSNP .part_output003.vcf.gz 2> hc003.log
198#There is no need to do variant calling on the NO_COOR shard, as there will
199# be no variants on unmapped reads
200#Merge output files into output VCF
201$SENTIEON_FOLDER/bin/sentieon driver --passthru --algo Haplotyper \
202   --merge output.vcf.gz .part_output000.vcf.gz .part_output001.vcf.gz \
203   .part_output002.vcf.gz .part_output003.vcf.gz 2> hc_merge.log
204以下是一个用于创建分片的Shell示例。generate_shards.sh
205determine_shards_from_bam(){
206    local bam step tag chr len pos end
207    bam="$1"
208    step="$2"
209    pos=1
210    samtools view -H $bam |\
211    while read tag chr len; do
212        [ $tag == '@SQ' ] || continue
213        chr=$(expr "$chr" : 'SN:\(.*\)')
214        len=$(expr "$len" : 'LN:\(.*\)')
215        while [ $pos -le $len ]; do
216            end=$(($pos + $step - 1))
217            if [ $pos -lt 0 ]; then
218                start=1
219            else
220                start=$pos
221            fi
222            if [ $end -gt $len ]; then
223                echo -n "$chr:$start-$len "
224                pos=$(($pos-$len))
225                break
226            else
227                echo "$chr:$start-$end"
228                pos=$(($end + 1))
229            fi
230        done
231    done
232    echo "NO_COOR"
233}
234
235determine_shards_from_dict(){
236    local bam step tag chr len pos end
237    dict="$1"
238    step="$2"
239    pos=1
240    cat $dict |\
241    while read tag chr len UR; do
242        [ $tag == '@SQ' ] || continue
243        chr=$(expr "$chr" : 'SN:\(.*\)')
244        len=$(expr "$len" : 'LN:\(.*\)')
245        while [ $pos -le $len ]; do
246            end=$(($pos + $step - 1))
247            if [ $pos -lt 0 ]; then
248                start=1
249            else
250                start=$pos
251            fi
252            if [ $end -gt $len ];then
253                echo -n "$chr:$start-$len "
254                pos=$(($pos-$len))
255                break
256            else
257                echo "$chr:$start-$end"
258                pos=$(($end + 1))
259            fi
260        done
261    done
262    echo "NO_COOR"
263}
264
265
266determine_shards_from_fai(){
267    local bam step tag chr len pos end
268    fai="$1"
269    step="$2"
270    pos=1
271    cat $fai |\
272    while read chr len other; do
273        while [ $pos -le $len ]; do
274            end=$(($pos + $step - 1))
275            if [ $pos -lt 0 ]; then
276                start=1
277            else
278                start=$pos
279            fi
280            if [ $end -gt $len ]; then
281                echo -n "$chr:$start-$len "
282                pos=$(($pos-$len))
283                break
284            else
285                echo "$chr:$start-$end"
286                pos=$(($end + 1))
287            fi
288        done
289    done
290    echo "NO_COOR"
291}
292
293if [ $# -eq 2 ]; then
294    filename=$(basename "$1")
295    extension="${filename##*.}"
296    if [ "$extension" = "fai" ]; then
297        determine_shards_from_fai $1 $2
298    elif [ "$extension" = "bam" ]; then
299        determine_shards_from_bam $1 $2
300    elif [ "$extension" = "dict" ]; then
301        determine_shards_from_dict $1 $2
302    fi
303else
304    echo "usage $0 file shard_size"
305fi

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

生物基因