满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）
衛(wèi)計(jì)委在2017年，2019年，2020年（還沒有答案）提供標(biāo)準(zhǔn)數(shù)據(jù)用于腫瘤生信分析的室間質(zhì)評。這樣預(yù)知結(jié)果的數(shù)據(jù)自然是不能放過了，本文嘗試參考GATK Best Practice：Somatic SNVs + Indels ，Cnvkit，Manta的pipeline來完成滿分流程分析，也可以使用標(biāo)準(zhǔn)數(shù)據(jù)反向判斷GATK Mutect2的實(shí)際準(zhǔn)確度,算法優(yōu)劣。
注：本文僅用于學(xué)習(xí)，距離真正的臨床應(yīng)用還有相當(dāng)大距離，歡迎大佬批評指正
**1. 分析流程概覽如下：
2. 本文用到的分析系統(tǒng)及分析流程文件
| SliverWorkspace 2.0.277363 | 提供運(yùn)行控制平臺(tái)/個(gè)人版 |
| FFPE SNV CNV SV V1.workflow | 分析流程文件，可以一鍵導(dǎo)入分析平臺(tái)(點(diǎn)擊查看操作) 當(dāng)然可以參照圖片中運(yùn)行腳本，shell里運(yùn)行，效果也是一樣 |
| 最終結(jié)果過濾腳本（python2.7 )及編譯版本 | Illumina_pt2.bed 等用到的bed，intelval等文件 SnvAnnotationFilter.py SNV過濾腳本 CnvAnnotationFilter.py CNV過濾腳本 SvAnnotationFilter.py SV 過濾腳本 QcProcessor.py 獲取整體QC數(shù)據(jù)的腳本 report_template.docx 分析報(bào)告模板 |
| 分析結(jié)果（pipeline結(jié)果與標(biāo)準(zhǔn)答案） | result.zip pipeline結(jié)果與標(biāo)準(zhǔn)答案 |
3. 本文用到的原始文件，用fastqc查看質(zhì)量狀態(tài)是clean data，Q值均高于30，這里就不需要去接頭和QC了。
百度網(wǎng)盤：鏈接: https://pan.baidu.com/s/1t9R14aQNP6Xk4tq1wFWJpg 提取碼: u5yp
| B1701_R1.fq.gz | 4.85G | 07d3cdccee41dbb3adf5d2e04ab28e5b |
| B1701_R2.fq.gz | 4.77G | c2aa4a8ab784c77423e821b9f7fb00a7 |
| B1701NC_R1.fq.gz | 3.04G | 4fc21ad05f9ca8dc93d2749b8369891b |
| B1701NC_R2.fq.gz | 3.11G | bc64784f2591a27ceede1727136888b9 |
4. 本文用到的環(huán)境變量（目錄/程序/文件/數(shù)值/字符）reference文件和數(shù)據(jù)庫體積過于龐大請自行下載安裝（如:ftp.broadinstitute.org/Annovar等等）
| sn | 1701 | 字符 |
| data | /opt/data | 目錄 |
| result | /opt/result | 目錄 |
| tools.java | /opt/jdk1.8.0_162/bin/java | 程序 |
| tools.bgzip | /opt/tabix-0.2.6/bgzip | 程序 |
| tools.bwa | /opt/bwa/bwa | 程序 |
| tools.cnvkit | /usr/local/bin/cnvkit.py | 程序 |
| tools.fastqc | /opt/FastQC/fastqc | 程序 |
| tools.gatk | /opt/ref/gatk-4.1.3.0/gatk | 程序 |
| tools.manta | /opt/manta/dist/bin/configManta.py | 程序 |
| tools.samtools | /opt/samtools/samtools | 程序 |
| tools.tabix | /opt/tabix-0.2.6/tabix | 程序 |
| tools.annovar | /opt/ref/annovar/table_annovar.pl | 程序 |
| tools.convertor | /opt/ref/annovar/convert2annovar.pl | 程序 |
| cutoff.cnv_max | 0.5 | 數(shù)值 |
| cutoff.cnv_min | -0.5 | 數(shù)值 |
| cutoff.cnvdep | 1000 | 數(shù)值 |
| cutoff.depth | 500 | 數(shù)值 |
| cutoff.event | 2 | 數(shù)值 |
| cutoff.nvaf | 0.005 | 數(shù)值 |
| cutoff.sv | 200 | 數(shù)值 |
| cutoff.TLOD | 16 | 數(shù)值 |
| cutoff.vaf | 0.01 | 數(shù)值 |
| envis.scatter | 8 | 數(shù)值 |
| envis.threads | 32 | 數(shù)值 |
| envis.memory | 32G | 字符 |
| refs.1000G | /opt/ref/1000G_phase1.indels.hg19.vcf | 文件 |
| refs.af_only | /opt/ref/af-only-gnomad.raw.sites.hg19.vcf.gz | 文件 |
| refs.bed | /opt/ref/projects/Illumina_pt2.bed | 文件 |
| refs.interval | /opt/ref/projects/Illumina_pt2.interval_list | 文件 |
| refs.dbsnp | /opt/ref/dbsnp_138.hg19.vcf | 文件 |
| refs.dict | /opt/ref/ucsc.hg19.dict | 文件 |
| refs.gene | /opt/ref/hg19_refGene.txt | 文件 |
| refs.hum | /opt/ref/ucsc.hg19.fa | 文件 |
| refs.mills | /opt/ref/Mills_and_1000G_gold_standard.indels.hg19.vcf | 文件 |
| refs.small.exac | /opt/ref/small_exac_common_3_b37.vcf | 文件 |
5. 詳細(xì)流程分析具體如下（不習(xí)慣圖形軟件的使用shell變量+腳本運(yùn)行也是一樣的）
GATK提供的標(biāo)準(zhǔn)流程從Normal/Tumor的fastq文件開始 map>reorder>sort>mark duplicate>recalibrator>apply BQSR
- 流程輸入文件
分析流程輸入文件，這里使用變量${sn}表示樣本編號，對室間質(zhì)評文件名做了調(diào)整。
-
Tumor 比對，管道操作給samtools，直接輸出bam格式文件。
-
Reorder Bam
GATK文檔中這樣的
Not to be confused with SortSam which sorts a SAM or BAM file with a valid sequence dictionary, ReorderSam reorders reads in a SAM/BAM file to match the contig ordering in a provided reference file, as determined by exact name matching of contigs. Reads mapped to contigs absent in the new reference are dropped. Runs substantially faster if the input is an indexed BAM file.
這一步省略其實(shí)對最終結(jié)果影響不大。
-
Sort
-
MarkDuplicate 標(biāo)記重復(fù)
-
重新校正堿基質(zhì)量值第一步，BaseRecalibrator：計(jì)算所有需要重校正的reads和特征值，然后把這些信息輸出為校準(zhǔn)表文件
-
重新校正堿基質(zhì)量值第二步，ApplyBQSR：用第一步得到的校準(zhǔn)表文件，重新調(diào)整BAM文件中的堿基質(zhì)量值，并使用這個(gè)新的質(zhì)量值重新輸出一個(gè)新的BAM文件。
-
Normal 數(shù)據(jù)對著Tumor的步驟重復(fù)一遍。序列比對
-
Normal Reorder
-
Normal Sort
-
Normal MarkDuplicate
- Normal BaseRecalibrator
-
Normal ApplyBQSR
-
Tumor GetPileupSummaries
-
Normal GetPileupSummaries
-
CalculateContamination
-
Mutect2
-
FilterMutectCalls 使用GATK提供的過濾器，過濾得到的突變
-
將過濾后的文件轉(zhuǎn)換為Annovar注釋所需要的格式
-
使用Annovar對結(jié)果注釋
-
使用自己寫的腳本對注釋后的結(jié)果過濾，比如按照室間質(zhì)評要求，過濾掉突變頻率低于1%，測序深度低于500的突變。對GATK某些過濾器過濾掉的結(jié)果進(jìn)行保留和排除，后面使用IGV進(jìn)行人工篩選。最終輸出的結(jié)果為，${sn}.result.SNV.xls（其實(shí)是個(gè)csv文件，擴(kuò)展名改為.xls是便于使用excel打開，很多人都這么干）
-
使用Manta Call SV，第一步，生成分析腳本文件runWorkflow.py
-
運(yùn)行runWorkflow.py獲取SV突變結(jié)果
-
使用自己過濾的腳本文件，對Manta獲取的SV過濾。如根據(jù)SOMATICSCORE分?jǐn)?shù)過濾，根據(jù)hg19_refGene.txt提供文件，計(jì)算突變基因等等。
-
使用CnvKIt，獲取CNV突變
-
使用CnvKit畫圖
-
使用py腳本文件，對CnvKit輸出結(jié)果過濾。同樣根據(jù)hg19_refGene.txt文件匹配基因，以及發(fā)生拷貝數(shù)變異的區(qū)域的外顯子區(qū)域等。
-
獲取整體QC數(shù)據(jù)，insertSize，depth等
- Samtools獲取堿基測序深度
- Samtools flagstat統(tǒng)計(jì)比對信息
- 使用py腳本文件對以上獲取的數(shù)據(jù)做統(tǒng)一處理，獲取整體QC狀態(tài)結(jié)果。
最終結(jié)果：
最終結(jié)果及各個(gè)命令/任務(wù)運(yùn)行時(shí)間如下：
整個(gè)分析過程耗時(shí) 3h 58min 29s，比較耗時(shí)，這個(gè)版本為了邏輯清晰一些，多數(shù)任務(wù)都是串行運(yùn)行，沒有很好利用計(jì)算資源。
后續(xù)文章會(huì)對整個(gè)流程進(jìn)行優(yōu)化（主要是并行化），最終分析時(shí)間在1h 13min左右（提升約3倍）。
GATK 輸出結(jié)果中SNV&INDEL的準(zhǔn)確度問題，經(jīng)過反復(fù)試驗(yàn)，不論如何設(shè)置過濾參數(shù)，最終的結(jié)果始終會(huì)有假陰性問題，這是GATK（4.1.3.0）中個(gè)別過濾器的問題，目前的補(bǔ)救措施是將部分GATK過濾器過濾掉的結(jié)果仍然包含在最終結(jié)果中，再使用IGV工具人工過濾（官方文檔也是這么推薦的），判斷該結(jié)果是否可靠。
總結(jié)
以上是生活随笔為你收集整理的满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【新版本】Aspose.Cells 10
- 下一篇： QQ互联第三方登录jar包