Introduction 介绍

GDC DNA-Seq分析管道识别全外显子组测序（WXS）和全基因组测序（WGS）数据中的体细胞变异。 通过比较正常和肿瘤样品比对中的等位基因频率，注释每个突变，并将来自多个病例的突变聚合到一个项目文件中来鉴定体细胞变体.

首先将内容比对到参考基因组上，接着进行co-cleaning以提高比对质量. 然后通过4种不同的方法分别识别体细胞突变并注释，最后所有pipline识别的突变将聚合进一个MAF文件中.

DNA-Seq 分析主要通过以下6个步骤:

• Genome Alignment
• Alignment Co-Cleaning
• Somatic Variant Calling
• Variant Annotation
• Mutation Aggregation
• Aggregated Mutation Masking

Data Processing Steps 数据处理步骤

Pre-Alignment 比对前处理

在比对之前，提交给GDC的BAM文件被读取、拆分并转换为FASTQ格式。 未通过Illumina chastity测试的reads将被删除。 请注意，此过滤步骤与使用基本质量分数修剪reads不同.

Alignment Workflow 比对流程

DNA-Seq 分析从 比对流程 开始. 取用 BWA 两种算法中的一种将reads组与参考基因组进行比对 [1]. 如果平均读取长度大于或等于70bp，则使用BWA-MEM，否则使用BWA-aln. 每个读取组分别与参考基因组对齐，并使用 Picard Tools SortSam and MergeSamFiles合并每个组分所有的reads. 然后标记可能作为PCR重复的reads以避免下游寻找突变的过程出错.

QUALITY CONTROL 质控

在比对前后收集质控指标并进行检查，以识别低质量数据。使用 FASTQC 对未比对的reads统计GC含量、平均读长及质量分数等基本指标. GDC对比对的reads收集质量指标，包括samtools idxstat 和 flagstatQuality. WGS及WXS的比对信息将通过Picard CollectMultipleMetrics 收集. WGS覆盖度通过picard CollectWgsMetrics 收集，而WXS覆盖度通过 picard CollectHsMetrics收集.

可以使用以下参数通过API访问每个文件端点的质量控制度量 expand=analysis.metadata.read_groups,analysis.metadata.read_groups.read_group_qcs 点击 这里 查看示例.

REFERENCE GENOME 参考基因组

使用人参考基因组GRCh38.d1.vd1 \进行所有比对。 诱饵病毒序列包括在参考基因组中以防止读数错误影响已知存在于人样品中的病毒的读数。 包括十种类型的人类病毒基因组：人巨细胞病毒（CMV），爱泼斯坦 - 巴尔病毒（EBV），乙型肝炎（HBV），丙型肝炎（HCV），人类免疫缺陷病毒（HIV），人类疱疹病毒8（HHV-8） ），人T淋巴细胞病毒1（HTLV-1），Merkel细胞多瘤病毒（MCV），Simian空泡病毒40（SV40）和人乳头瘤病毒（HPV）。 GDC使用的参考序列可以这里下载.

DNA-Seq 比对的命令行参数

请注意，由于管道开发和改进正在进行，从GDC门户下载的文件中的版本号可能会有所不同.

STEP 1: 使用BIOBAMBAM将BAMS转换为FASTQS - BIOBAMBAM2 2.0.54

• Shell

bamtofastq \
collate=1 \
exclude=QCFAIL,SECONDARY,SUPPLEMENTARY \
filename= <input.bam> \
gz=1 \
inputformat=bam
level=5 \
outputdir= <output_path> \
outputperreadgroup=1 \
outputperreadgroupsuffixF=_1.fq.gz \
outputperreadgroupsuffixF2=_2.fq.gz \
outputperreadgroupsuffixO=_o1.fq.gz \
outputperreadgroupsuffixO2=_o2.fq.gz \
outputperreadgroupsuffixS=_s.fq.gz \
tryoq=1 \

STEP 2: BWA 比对 - BWA 0.7.15 - SAMTOOLS 1.3.1

如果平均读取长度大于或等于70bp:

• Shell

bwa mem \
-t 8 \
-T 0 \
-R <read_group> \
<reference> \
<fastq_1.fq.gz> \
<fastq_2.fq.gz> |
samtools view \
-Shb
-o <output.bam> -

如果平均读取长度小于70bp:

• Shell

bwa aln -t 8 <reference> <fastq_1.fq.gz> > <sai_1.sai> &&
bwa aln -t 8 <reference> <fastq_2.fq.gz> > <sai_2.sai> &&
bwa sampe -r <read_group> <reference> <sai_1.sai> <sai_2.sai> <fastq_1.fq.gz> <fastq_2.fq.gz> | samtools view -Shb -o <output.bam> -

如果质量分数编码为Illumina 1.3或1.5，请使用带有“-l”标志的BWA aln.

STEP 3: BAM SORT 排序 - PICARD 2.6.0

• Shell

java -jar picard.jar SortSam \
CREATE_INDEX=true \
INPUT=<input.bam> \
OUTPUT=<output.bam> \
SORT_ORDER=coordinate \
VALIDATION_STRINGENCY=STRICT

STEP 4: BAM MERGE 合并 - PICARD 2.6.0

• Shell

java -jar picard.jar MergeSamFiles \
ASSUME_SORTED=false \
CREATE_INDEX=true \                 
[INPUT= <input.bam>]  \
MERGE_SEQUENCE_DICTIONARIES=false \
OUTPUT= <output_path> \
SORT_ORDER=coordinate \
USE_THREADING=true \
VALIDATION_STRINGENCY=STRICT

STEP 5: MARK DUPLICATES 标记重复 - PICARD 2.6.0

• Shell

java -jar picard.jar MarkDuplicates \
CREATE_INDEX=true \
INPUT=<input.bam> \
VALIDATION_STRINGENCY=STRICT

Co-cleaning Workflow

Co-cleaning workflow 进一步提高了比对质量. Co-cleaning 需要以单独的管道进行，因为他需要使用与同一患者相关的多个BAM文件（即肿瘤BAM和正常组织BAM）. 这两个步骤都可以通过 GATK 实现.

INDEL 局部重比对

使用 IndelRealigner 进行插入和缺失的局部重新排列. 该步骤可以找到包含跨BAM文件的错误比对区域，这通常由相对于参考基因组的插入 - 缺失（indel）突变引起. indel突变的错位（通常可被错误地评分为替换）降低了下游调用步骤的准确性。.

碱基质量分数重新校准

使用 BaseRecalibrator 执行基本质量得分重新校准（BQSR）步骤. 此步骤基于可检测的和系统的错误调整基本质量分数，还提高了下游变体调用算法的准确性. 请注意，原始质量分数保存在co-cleaned的BAM文件的OQ字段中。 如果需要将BAM文件转换为FASTQ格式，则应使用这些分数.

DNA-Seq Co-Cleaning 命令行参数

STEP 1: REALIGNTARGETCREATOR 重比对创建

• Shell

java -jar GenomeAnalysisTK.jar \
-T RealignerTargetCreator \
-R <reference>
-known <known_indels.vcf>
[ -I <input.bam> ]
-o <realign_target.intervals>

STEP 2: INDELREALIGNER 插入重比对

• Shell

java -jar GenomeAnalysisTK.jar \
-T IndelRealigner \
-R <reference> \
-known <known_indels.vcf> \
-targetIntervals <realign_target.intervals> \
--noOriginalAlignmentTags \
[ -I <input.bam> ] \
-nWayOut <output.map>

STEP 3: BASERECALIBRATOR 碱基重校准; DBSNP V.144

• Shell

java -jar GenomeAnalysisTK.jar \
-T BaseRecalibrator \
-R <reference> \
-I <input.bam> \
-knownSites <dbsnp.vcf>
-o <bqsr.grp>

STEP 4: PRINTREADS 输出reads

• Shell

java -jar GenomeAnalysisTK.jar \
-T PrintReads \
-R <reference> \
-I <input.bam> \
--BQSR <bqsr.grp> \
-o <output.bam>

Somatic Variant Calling Workflow 体细胞突变寻找流程

比对和 co-cleaned 后的 BAM 文件通过 Somatic Mutation Calling Workflow 处理 tumor-normal pairs. 使用4个单独管道进行突变查找:

• MuSE [2]
• MuTect2 [3]
• VarScan2 [4]
• SomaticSniper [5]

每个管道都会将突变信息记录在VCF格式文件中. 每个可用字段的详细信息，请参阅GDC VCF Format 文档. 此时，在DNA-Seq管道中，所有下游分析都分为四个独立的路径，这些路径对应于它们各自的突变调用管道.

PIPELINE DESCRIPTIONS 管道描述

4个分离的突变寻找管道实现了GDC数据的统一化. 目前没有The一个管道被统一认为是最佳的，所以这是由研究人员负责选择的最适合数据的管道。有关管道的一些细节如下所示.

MuTect2 pipeline 使用“正常组”来识别其他种系突变. 该小组由来自数千名TCGA血液正常基因组个体产生，所述个体以高可信度被评估为无癌症。 该方法允许将更高水平的置信度分配给MuTect2管道调用的体细胞变体.

其他三个管道的基本轮廓可以在这里找到:

• VarScan2 pipeline
• MuSE pipeline
• SomaticSniper pipeline

INDELS 插入和缺失

MuTect2 和 VarScan 检测出来的插入和缺失会被记录在GDC VCF 文件中.

GERMLINE VARIANTS 种系突变

此时，种系突变被故意排除为协调数据. GDC不建议使用以前检测到并存储在Legacy Archive中的种系突变，因为它们不符合GDC高质量数据标准.

Variant Call 命令行参数

MUSE

MuSEv1.0rc_submission_c039ffa; dbSNP v.144

Step 1: MuSE call

• Shell

MuSE call \
-f <reference> \
-r <region> \                                   
<tumor.bam> \
<normal.bam> \
-O <intermediate_muse_call.txt>

Step 2: MuSE sump

• Shell

MuSE sump \
-I <intermediate_muse_call.txt> \                           
-E \                                
-D <dbsnp_known_snp_sites.vcf> \
-O <muse_variants.vcf>  

Note: -E is used for WXS data and -G can be used for WGS data.

MUTECT2

GATK nightly-2016-02-25-gf39d340; dbSNP v.144

• Shell

java -jar GenomeAnalysisTK.jar \
-T MuTect2 \
-R <reference> \
-L <region> \
-I:tumor <tumor.bam> \
-I:normal <normal.bam> \
--normal_panel <pon.vcf> \                        
--cosmic <cosmic.vcf> \
--dbsnp <dbsnp.vcf> \
--contamination_fraction_to_filter 0.02 \                   
-o <mutect_variants.vcf> \
--output_mode EMIT_VARIANTS_ONLY \
--disable_auto_index_creation_and_locking_when_reading_rods

SOMATICSNIPER

Somatic-sniper v1.0.5.0

• Shell

bam-somaticsniper \
-q 1 \
-L \
-G \
-Q 15 \
-s 0.01 \
-T 0.85 \
-N 2 \
-r 0.001 \
-n NORMAL \
-t TUMOR \
-F vcf \
-f ref.fa \
<tumor.bam> \
<normal.bam> \
<somaticsniper_variants.vcf>

VARSCAN

Step 1: Mpileup; Samtools 1.1

• Shell

samtools mpileup \
-f <reference> \
-q 1 \
-B \
<normal.bam> \
<tumor.bam> >
<intermediate_mpileup.pileup>

Step 2: Varscan Somatic; Varscan.v2.3.9

• Shell

java -jar VarScan.jar somatic \
<intermediate_mpileup.pileup> \
<output_path> \
--mpileup      1 \
--min-coverage 8 \
--min-coverage-normal 8 \
--min-coverage-tumor 6 \
--min-var-freq 0.10 \
--min-freq-for-hom 0.75 \
--normal-purity 1.0 \
--tumor-purity 1.00 \
--p-value 0.99 \
--somatic-p-value 0.05 \
--strand-filter 0 \
--output-vcf

Step 3: Varscan ProcessSomatic; Varscan.v2.3.9

• Shell

java -jar VarScan.jar processSomatic \
<intermediate_varscan_somatic.vcf> \
--min-tumor-freq 0.10 \
--max-normal-freq 0.05 \
--p-value 0.07

Variant Call 注释流程

在 Somatic Annotation Workflow 中使用 Variant Effect Predictor (VEP) v84 [6] 和VEP GDC插件对原始VCF文件进行注释.

VEP使用VCF文件中的坐标和等位基因来推断每个变体的生物学背景，包括每个突变的位置，其生物学后果（移码/沉默突变）和受影响的基因。 有关详细信息，请参阅 GDC VCF Format 上的文档。 VCF文件中的突变也与来自外部变异数据库的已知突变匹配， 以下数据库用于VCF注释:

• GENCODE v.22
• sift v.5.2.2
• ESP v.20141103
• polyphen v.2.2.2
• dbSNP v.146
• Ensembl genebuild v.2014-07
• Ensembl regbuild v.13.0
• HGMD public v.20154
• ClinVar v.201601

由于许可限制，COSMIC不用于GDC VEP工作流程中的注释.

除注释外， False Positive Filter 用于标记VarScan和SomaticSniper输出中的低质量变体。 SomaticSniper中SSQ <25的变体也被删除.

Tumor-Only Variant Calling Workflow 肿瘤突变寻找流程

Tumor only突变寻找 是对没有正常人作参照的样本进行处理。该方法利用了大多数肿瘤样品中存在的正常细胞污染的特点。这些calling是使用GATK4中包含的MuTect2版本进行的。Tumor-only variant call 文件可以在 GDC Portal 中以关键字 "Workflow Type: GATK4 MuTect2" 找到.

Tumor-Only Variant Call 命令行参数

• Source

GATK4 v4.0.4.0

## 1\. Generate OXOG metrics:

java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
CollectSequencingArtifactMetrics \
-I Tumor_Sample_Alignment.bam \
-O <job_identifier> \
--FILE_EXTENSION .txt \
-R GRCh38.d1.vd1.fa  ## Only chr1-22 + XYM

## 2\. Generate pileup summaries on tumor sample:

java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
GetPileupSummaries
-I Tumor_Sample_Alignment.bam \
-O <job_identifier>.targeted_sequencing.table \
-V af-only-gnomad-common-biallelic.grch38.main.vcf.gz \ # Germline reference from gnomad
-L intervals.bed \ ## Only chr1-22 + XYM
-R GRCh38.d1.vd1.fa

## 3\. Calculate contamination on tumor sample

java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
CalculateContamination \
-I <job_identifier>.targeted_sequencing.table \ # From step 2
-O <job_identifier>.targeted_sequencing.contamination.table

## 4\. Find tumor sample name from BAM

java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
GetSampleName \
-I Tumor_Sample_Alignment.bam \
-O <job_identifier>.targeted_sequencing.sample_name

## 5\. Run MuTect2 using only tumor sample on chromosome level (25 commands with different intervals)

java -Djava.io.tmpdir=/tmp/job_tmp_3 -d64 -jar -Xmx3G -XX:+UseSerialGC \
/bin/gatk-4.0.4.0/gatk-package-4.0.4.0-local.jar \
Mutect2 \
-R GRCh38.d1.vd1.fa \
-L chr4:1-190214555 \ # Specify chromosome
-I Tumor_Sample_Alignment.bam \
-O 3.mt2.vcf \
-tumor <tumor_sample_name> \ # From step 4
--af-of-alleles-not-in-resource 2.5e-06 \
--germline-resource af-only-gnomad.hg38.vcf.gz \ # Germline reference from gnomad
-pon gatk4_mutect2_4136_pon.vcf.gz # New panel of normal created by 4136 TCGA curated normal samples, using GATK4

## After this step, all chromosome level VCFs are merged into one.

## 6\. Sort VCF with Picard

java -d64 -XX:+UseSerialGC -Xmx16G -jar /usr/local/bin/picard.jar \
SortVcf \
SEQUENCE_DICTIONARY=GRCh38.d1.vd1.dict \
OUTPUT=<job_identifier>.targeted_sequencing.mutect2.tumor_only.sorted.vcf.gz \
I=merged_multi_gatk4_mutect2_tumor_only_calling.vcf \ # From step 5
CREATE_INDEX=true

## 7\. Filter variant calls from MuTect
java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
FilterMutectCalls \
-O <job_identifier>.targeted_sequencing.mutect2.tumor_only.contFiltered.vcf.gz \
-V <job_identifier>.targeted_sequencing.mutect2.tumor_only.sorted.vcf.gz \ # From step 6
--contamination-table <job_identifier>.targeted_sequencing.contamination.table \ # From step 3
-L intervals.bed

## 8\. Filter variants by orientation bias
java -d64 -XX:+UseSerialGC -Xmx3G -jar /gatk/gatk.jar \
FilterByOrientationBias \
-O <job_identifier>.targeted_sequencing.tumor_only.gatk4_mutect2.raw_somatic_mutation.vcf.gz \ # final output
-P <job_identifier>.pre_adapter_detail_metrics.txt \ # From step 1
-V <job_identifier>.targeted_sequencing.mutect2.tumor_only.contFiltered.vcf.gz \ # From step 7
-L intervals.bed \
-R GRCh38.d1.vd1.fa \
-AM G/T \
-AM C/T

Tumor-Only Variant 注释流程

在用MuTect2进行单肿瘤突变calling后，一系列的过滤器将被用以最小化可下载的VCF文件中的种系突变的release。 在所有案例中, GDC根据等位基因频率，映射质量，体细胞/种系概率和拷贝数自定义一组过滤器。在一些案例中，过滤之前会有额外的突变聚类步骤 .

PureCN R-package [7] [8] 用于根据肿瘤纯度，倍性，污染，拷贝数和杂合性缺失通过体细胞/种系状态和克隆性对变体进行分类。 使用该包执行以下内容步骤:

• Interval Capture : 使用FASTA和BED文件坐标生成间隔文件.
• GC-Normalization : 计算GC标准化肿瘤/正常覆盖率数据.
• Normal DB Creation : 用规范化的覆盖文件及 panel-of-normals VCF生成 normal database
• Somatic Variant Calling : 对先前识别的突变进行分类

请注意，如果没有足够的数据来生成目标捕获工具包特定的正常数据库，则不会执行PureCN。 在极少数情况下，PureCN可能找不到数字解决方案。 如果未执行PureCN或找不到解决方案，则会在VCF标头中指出。使用这些管道注释的VCF文件可在 GDC Portal 键入关键字 "Workflow Type: GATK4 MuTect2 Annotation" 找到.

Somatic Aggregation Workflow 体细胞聚合流程

Somatic Aggregation Workflow从多个VCF文件生成一个MAF文件;有关文件结构的详细信息，请参阅 GDC MAF Format 格式指南. 在此步骤中，为每个项目的每个突变调用管道生成一个MAF文件，并包含此项目中的所有可用案例。

Masked Somatic Aggregation Workflow

由于存在种系突变，Somatic Aggregation Workflow生成的MAF文件是受控访问的。通过删除可能包含种系突变信息的列和变体，可以修改开放式MAF文件以进行公开发布。 有关用于删除变体的条件的详细信息，请参阅 GDC MAF Format 格式。

虽然这些标准管道导致开放式MAF文件中的一些真正的阳性体细胞突变过度过滤，但它们会阻止个人可识别的种系突变信息公开。 如果遗漏某些体细胞突变是一个问题，GDC建议研究人员探索受控和开放获取的MAF文件。

File Access and Availability 文件访问及可用性

来自GDC DNA-Seq分析管道的文件可在BAM，VCF和MAF格式的GDC数据门户中获得。 下面列出了所有可用数据类型及其各自文件格式的说明

Pipline source ：GDC