全基因組 - 人類基因組變異分析（PacBio) （4）-- DeepVariant

一、SNPs和INDELs變異檢測(cè)

單核苷酸多態(tài)性（Single Nucleotide Polymorphism，SNP）指的是基因組中單個(gè)核苷酸腺嘌呤(A)、胸腺嘧啶(T)、胞嘧啶(C)或鳥(niǎo)嘌呤(G)在物種成員之間或個(gè)體配對(duì)染色體之間的差異, 是最常見(jiàn)也最簡(jiǎn)單的一類造成基因組多樣性的DNA序列變異。

插入缺失（insertion-deletion，InDel），這里一般指小于50bp的變異，即在DNA序列中添加或刪除少量堿基，主要指在基因組某個(gè)位置上發(fā)生較短長(zhǎng)度的線性片段插入(Insert)或者缺失(Deletion)的現(xiàn)象，合稱為InDel。

我們對(duì)下機(jī)數(shù)據(jù)進(jìn)行比對(duì)分析 (pbmm2軟件)，提取全基因組中所有的潛在多態(tài)性SNP位點(diǎn)和小片段插入/缺失InDel位點(diǎn)（DeepVariant軟件），后期再根據(jù)質(zhì)量值、深度、重復(fù)性等因素做進(jìn)一步的過(guò)濾篩選，最終得到高可信度的SNP和InDel數(shù)據(jù)集并注釋1。

SNP和INDEL變異檢測(cè)有助于我們更深入地了解基因組，生物性狀的表現(xiàn)，物種的起源與進(jìn)化，認(rèn)識(shí)基因變異和疾病的之間的聯(lián)系。從測(cè)序數(shù)據(jù)中進(jìn)行準(zhǔn)確的變異檢測(cè)也是生物學(xué)、醫(yī)學(xué)研究和精準(zhǔn)醫(yī)學(xué)的基礎(chǔ)1。

我們對(duì)下機(jī)數(shù)據(jù)進(jìn)行比對(duì)分析，提取全基因組中所有的潛在多態(tài)性SNP位點(diǎn)和小片段插入/缺失InDel位點(diǎn)，再根據(jù)質(zhì)量值、深度、重復(fù)性等因素做進(jìn)一步的過(guò)濾篩選，最終得到高可信度的SNP數(shù)據(jù)集并注釋。

二、變異檢測(cè)工具

目前SNP和INDEL變異檢測(cè)的軟件有很多，比如老牌行業(yè)“金標(biāo)” 由Broad Institute 開(kāi)發(fā) GATK，即Genome Analysis Toolkit 基因組分析工具。在變異軟件綜合評(píng)測(cè)中2，3，DeepVariant軟件在三代測(cè)序數(shù)據(jù)中表現(xiàn)是非常優(yōu)秀的 （圖1，圖2，圖3）。

圖1. 參考文獻(xiàn)2文章摘要
圖2. 參考文獻(xiàn)2文章中評(píng)測(cè)的軟件
圖3. DeepVariant 原文評(píng)測(cè)軟件

PacBio生信分析培訓(xùn)推薦DeepVariant作為SNP和INDEL變異檢測(cè)的軟件，并且對(duì)于小型變異檢測(cè)PacBio官方推薦的也是DeepVariant（圖4）, 所以接下來(lái)我們?cè)敿?xì)介紹下DeepVariant的使用方法。

圖4. PacBio官方推薦

三、DeepVariant簡(jiǎn)介

官方網(wǎng)址：https://github.com/google/deepvariant

DeepVariant是由谷歌Google基于深度卷積神經(jīng)網(wǎng)絡(luò)開(kāi)發(fā)的一款從DNA測(cè)序數(shù)據(jù)中快速較精確識(shí)別堿基變異位點(diǎn)的軟件4。這個(gè)工具在準(zhǔn)確率上和精確度上，比傳統(tǒng)的比對(duì)拼接方法都高出一大截。DeepVariant，把工作量巨大的拼接問(wèn)題（高通量測(cè)序碎片化的結(jié)果拼接成完整的基因序列），轉(zhuǎn)變成了一個(gè)典型的圖像分類問(wèn)題。而圖像分類正是谷歌擅長(zhǎng)的技術(shù)。在2016 PrecisionFDA的Truth Challenge比賽中，DeepVariant獲得了最高SNP性能獎(jiǎng)，PacBio +DeepVariant（Highest SNP Performance）（圖5）。在那之后，Google Brain團(tuán)隊(duì)又將錯(cuò)誤率降低了50%。

圖5. Variant-Calling-Performance-HG003

DeepVarient軟件運(yùn)行運(yùn)行流程如下圖6所示：

圖6. DeepVarient軟件運(yùn)行流程

左邊：篩選候選的變異位點(diǎn)集合；中間：SNN訓(xùn)練樣本；右邊：用訓(xùn)練好的模型判斷Genotype

四、DeepVariant安裝、使用及實(shí)操

1.軟件安裝

DeepVarient官方提供了3種安裝方式：

  • Docker，官方推薦使用的安裝方式。
  • 源代碼構(gòu)建。
  • 二進(jìn)制文件。

使用官方推薦安裝方式pull docker鏡像

#指定DeepVariant版本，這次安裝為最新版本v1.6.0, 2023年10月24日發(fā)布。
$ BIN_VERSION="1.6.0"

#拉取docker鏡像，大小為5.74GB
$ docker pull google/deepvariant:"${BIN_VERSION}"

2.數(shù)據(jù)準(zhǔn)備

  • 樣本參考基因組文件

例如上一節(jié)pbmm2用到的GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz.
參考基因組需要samtools進(jìn)行索引

#如果沒(méi)有安裝samtools，可以使用conda進(jìn)行安裝。
$ conda install -c bioconda samtools

#GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz 已經(jīng)重命名為GRCh38.fa
#samtools進(jìn)行索引
$ cd Human_ref
$ samtools faidx GRCh38.fa

#查看索引后的文件
$ ls
# GRCh38.fa  GRCh38.fa.fai
  • 比對(duì)排序后的BAM文件

例如, 經(jīng)過(guò)pbmm2比對(duì)，排序，索引后的bam文件：

HG002_1.align.bam
HG002_1.align.bam.bai
HG003.align.bam
HG003.align.bam.bai
HG004.align.bam
HG004.align.bam.bai

3.運(yùn)行DeepVariant

#使用方法
$ BIN_VERSION="1.6.0"
$ docker run \
  -v "YOUR_INPUT_DIR":"/input" \
  -v "YOUR_OUTPUT_DIR:/output" \
  google/deepvariant:"${BIN_VERSION}" \    #根據(jù)DeepVariant版本號(hào)來(lái)設(shè)置
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \      #根據(jù)應(yīng)用替換為其中的一種[WGS,WES,PACBIO,ONT_R104,HYBRID_PACBIO_ILLUMINA]  WGS和WES為二代Illumina數(shù)據(jù)
  --ref=/input/YOUR_REF \           #samtools索引的參考基因組
  --reads=/input/YOUR_BAM \     #比對(duì)過(guò)的bam文件
  --output_vcf=/output/YOUR_OUTPUT_VCF \   #輸出的vcf
  --output_gvcf=/output/YOUR_OUTPUT_GVCF \  #輸出的gvcf, 用于合并樣本進(jìn)行變異分析
  --num_shards=$(nproc) \   #CPU核數(shù)
  --logging_dir=/output/logs \  #每一步的日志文檔
  --haploid_contigs="chrX,chrY"   #如果用GRCh38則設(shè)置為"chrX,chrY"，GRCh37則設(shè)置為"X,Y"。

#  --regions "chr20:10,000,000-10,010,000"  可以指定變異檢測(cè)范圍

HG002_1，HG003，HG004實(shí)際樣本運(yùn)行：


#HG002_1
$ nohup docker run \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/input" \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/output" \
  google/deepvariant:"1.6.0" \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=PACBIO \
  --ref=/input/Human_ref/GRCh38.fa \
  --reads=/input/HG002_1.align.bam \
  --output_vcf=/output/HG002_1.vcf.gz \
  --output_gvcf=/output/HG002_1.gvcf.gz \
  --num_shards=20 \
  --logging_dir=/output/HG002_1_DeepV_logs \
  --haploid_contigs="chrX,chrY" &

#HG003
$ nohup docker run \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/input" \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/output" \
  google/deepvariant:"1.6.0" \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=PACBIO \
  --ref=/input/Human_ref/GRCh38.fa \
  --reads=/input/HG003.align.bam \
  --output_vcf=/output/HG003.vcf.gz \
  --output_gvcf=/output/HG003.gvcf.gz \
  --num_shards=20 \
  --logging_dir=/output/HG003_DeepV_logs \
  --haploid_contigs="chrX,chrY" &

#HG004
$ nohup docker run \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/input" \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/output" \
  google/deepvariant:"1.6.0" \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=PACBIO \
  --ref=/input/Human_ref/GRCh38.fa \
  --reads=/input/HG004.align.bam \
  --output_vcf=/output/HG004.vcf.gz \
  --output_gvcf=/output/HG004.gvcf.gz \
  --num_shards=20 \
  --logging_dir=/output/HG004_DeepV_logs \
  --haploid_contigs="chrX,chrY" &


# nohup 掛起不間斷
# &放入后臺(tái)運(yùn)行

4.合并多個(gè)樣本的gvcf文件

多樣本聯(lián)合變異（Joint Calling）需要用到 GLnexus工具。

GLnexus是由DNAnexus開(kāi)發(fā)，用于可擴(kuò)展的gVCF合并和聯(lián)合變異（joint calling）要求群體測(cè)序項(xiàng)目，GL即genotype likelihood之意5。

GATK作為變異檢測(cè)金標(biāo)準(zhǔn)軟件，缺點(diǎn)在于速度很慢。盡管在分析上游的比對(duì)、單樣本HaplotypeCaller檢測(cè)等環(huán)節(jié)已經(jīng)有很多替代品，如Sentieon、Parabricks等，但下游Joint Calling這一步優(yōu)化仍然有限。而這正是整個(gè)GATK流程的最限速的步驟，在GATK中只能通過(guò)分區(qū)的方法來(lái)加速，效果非常有限5。

GLnexus的開(kāi)發(fā)解決了這個(gè)痛點(diǎn)問(wèn)題，在速度上不說(shuō)幾十上百倍的提升，至少也有十多倍。對(duì)于大規(guī)模群體/隊(duì)列而言（主要針對(duì)人類基因組開(kāi)發(fā)），是個(gè)非常好的工具5。

Deepvariant和Clara Parabricks都推薦它來(lái)做聯(lián)合變異5。

#下載GLnexus docker鏡像
$ docker pull quay.io/mlin/glnexus:v1.2.7

#運(yùn)行代碼
$ sudo docker run \
   -v "${DIR}":"/data" \
  quay.io/mlin/glnexus:v1.2.7 \
  /usr/local/bin/glnexus_cli \
  --config DeepVariant \
  --bed "/data/${CAPTURE_BED}" \
  /data/HG004.g.vcf.gz /data/HG003.g.vcf.gz /data/HG002.g.vcf.gz \
  | sudo docker run -i google/deepvariant:${VERSION} bcftools view - \
  | sudo docker run -i google/deepvariant:${VERSION} bgzip -c \
  > ${DIR}/deepvariant.cohort.vcf.gz

幫助文檔：

#查看幫助文檔
$ docker run \
 quay.io/mlin/glnexus:v1.2.7 \
  /usr/local/bin/glnexus_cli  --help

##幫助文檔
Merge and joint-call input gVCF files, emitting multi-sample BCF on standard output.

Options:
  --dir DIR, -d DIR              scratch directory path (mustn't already exist; default: ./GLnexus.DB)
  --config X, -c X               configuration preset name or .yml filename (default: gatk)
  --bed FILE, -b FILE            three-column BED file with ranges to analyze (if neither --range nor --bed: use full length of all contigs)
  --list, -l                     expect given files to contain lists of gVCF filenames, one per line

  --more-PL, -P                  include PL from reference bands and other cases omitted by default
  --squeeze, -S                  reduce pVCF size by suppressing detail in cells derived from reference bands
  --trim-uncalled-alleles, -a    remove alleles with no output GT calls in postprocessing

  --mem-gbytes X, -m X           memory budget, in gbytes (default: most of system memory)
  --threads X, -t X              thread budget (default: all hardware threads)

  --help, -h                     print this help message

Configuration presets:
            Name          CRC32C    Description
            gatk       268790301    Joint-call GATK-style gVCFs
 gatk_unfiltered       484625853    Merge GATK-style gVCFs with no QC filters or genotype revision
          xAtlas      3445896276    Joint-call xAtlas gVCFs
xAtlas_unfiltered       920229760   Merge xAtlas gVCFs with no QC filters or genotype revision
          weCall      2936345659    Joint-call weCall gVCFs
weCall_unfiltered      2838640462   Merge weCall gVCFs with no filtering or genotype revision
     DeepVariant      2857227159    Joint call DeepVariant whole genome sequencing gVCFs
  DeepVariantWGS      2857227159    Joint call DeepVariant whole genome sequencing gVCFs
  DeepVariantWES      4105299981    Joint call DeepVariant whole exome sequencing gVCFs
DeepVariant_unfiltered       116393675  Merge DeepVariant gVCFs with no QC filters or genotype revision
        Strelka2      3838963651    [EXPERIMENTAL] Merge Strelka2 gVCFs with no QC filters or genotype revision

HG002_1，HG003，HG004實(shí)際樣本gvcf合并：

根據(jù)幫助文檔 --config 可以設(shè)置為DeepVariant ，DeepVariant和DeepVariantWGS應(yīng)該是一樣的preset。如果使用GATK或其它軟件進(jìn)行變異分析則設(shè)置相應(yīng)的模型

#joint-calling，最后得到deepvariant.cohort.vcf.gz
$ sudo docker run \
  -v "/mnt/data/home/mli/Desktop/pb_WGS":"/data" \
  quay.io/mlin/glnexus:v1.2.7 \
  /usr/local/bin/glnexus_cli \
  --config DeepVariant \
  /data/HG002_1.gvcf.gz /data/HG003.gvcf.gz /data/HG004.gvcf.gz \
  | docker run -i google/deepvariant:1.6.0 bcftools view - \
  | docker run -i google/deepvariant:1.6.0 bgzip -c \
  > deepvariant.cohort.vcf.gz

#未設(shè)置
  --bed "/data/${CAPTURE_BED}"

五、結(jié)果說(shuō)明

DeepVariant軟件輸出結(jié)果為vcf格式文件（圖7），相信做生物信息的小伙伴都很熟悉了，這里不再贅述1。

HG002_1單個(gè)樣本vcf:

圖7.HG002 vcf 結(jié)果

HG002_1, HG003, HG004 三個(gè)樣本joint calling的deepvariant.cohort.vcf.gz:

圖8.Joint Calling vcf 結(jié)果

特別說(shuō)明：
DeepVariant運(yùn)行結(jié)束后，SNP和INDEL還在一個(gè)vcf文件里，如果為了后續(xù)單獨(dú)分析，我們可以用GATK分離他們，命令如下1：

$ gatk  SelectVariants -R /input/YOUR_REF  \   
    -V /input/YOUR_VCF  \
    -O /output/YOUR_RAW_SNP  \
    -select-type SNP|INDEL

-select-type參數(shù)分別給定SNP和INDEL，將會(huì)分別得到對(duì)應(yīng)變異類型的結(jié)果，輸出仍然是vcf格式的文件。有了這個(gè)結(jié)果后，就可以進(jìn)行后續(xù)的分析了。

參考文獻(xiàn)：

  1. 三代變異檢測(cè)操作說(shuō)明-DeepVariant
  2. Pei, S. et.al. (2021). Benchmarking variant callers in next-generation and third-generation sequencing analysis. Briefings in Bioinformatics.
  3. Poplin, R.et.al. (2018). A universal SNP and small-indel variant caller using deep neural networks. Nature Biotechnology
  4. 量子位-谷歌推出開(kāi)源工具DeepVariant，用深度學(xué)習(xí)識(shí)別基因變異
  5. 生物信息與育種：GLNexus進(jìn)行joint calling時(shí)的"half-calls"（如./0, ./1）問(wèn)題？
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容