满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）
衛計委在2017年，2019年，2020年（還沒有答案）提供標準數據用于腫瘤生信分析的室間質評。這樣預知結果的數據自然是不能放過了，本文嘗試參考GATK Best Practice：Somatic SNVs + Indels ，Cnvkit，Manta的pipeline來完成滿分流程分析，也可以使用標準數據反向判斷GATK Mutect2的實際準確度,算法優劣。
注：本文僅用于學習，距離真正的臨床應用還有相當大距離，歡迎大佬批評指正
**1. 分析流程概覽如下：
2. 本文用到的分析系統及分析流程文件
| SliverWorkspace 2.0.277363 | 提供運行控制平臺/個人版 |
| FFPE SNV CNV SV V1.workflow | 分析流程文件，可以一鍵導入分析平臺(點擊查看操作) 當然可以參照圖片中運行腳本，shell里運行，效果也是一樣 |
| 最終結果過濾腳本（python2.7 )及編譯版本 | Illumina_pt2.bed 等用到的bed，intelval等文件 SnvAnnotationFilter.py SNV過濾腳本 CnvAnnotationFilter.py CNV過濾腳本 SvAnnotationFilter.py SV 過濾腳本 QcProcessor.py 獲取整體QC數據的腳本 report_template.docx 分析報告模板 |
| 分析結果（pipeline結果與標準答案） | result.zip pipeline結果與標準答案 |
3. 本文用到的原始文件，用fastqc查看質量狀態是clean data，Q值均高于30，這里就不需要去接頭和QC了。
百度網盤：鏈接: 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. 本文用到的環境變量（目錄/程序/文件/數值/字符）reference文件和數據庫體積過于龐大請自行下載安裝（如: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 | 數值 |
| cutoff.cnv_min | -0.5 | 數值 |
| cutoff.cnvdep | 1000 | 數值 |
| cutoff.depth | 500 | 數值 |
| cutoff.event | 2 | 數值 |
| cutoff.nvaf | 0.005 | 數值 |
| cutoff.sv | 200 | 數值 |
| cutoff.TLOD | 16 | 數值 |
| cutoff.vaf | 0.01 | 數值 |
| envis.scatter | 8 | 數值 |
| envis.threads | 32 | 數值 |
| 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. 詳細流程分析具體如下（不習慣圖形軟件的使用shell變量+腳本運行也是一樣的）
GATK提供的標準流程從Normal/Tumor的fastq文件開始 map>reorder>sort>mark duplicate>recalibrator>apply BQSR
- 流程輸入文件
分析流程輸入文件，這里使用變量${sn}表示樣本編號，對室間質評文件名做了調整。
-
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.
這一步省略其實對最終結果影響不大。
-
Sort
-
MarkDuplicate 標記重復
-
重新校正堿基質量值第一步，BaseRecalibrator：計算所有需要重校正的reads和特征值，然后把這些信息輸出為校準表文件
-
重新校正堿基質量值第二步，ApplyBQSR：用第一步得到的校準表文件，重新調整BAM文件中的堿基質量值，并使用這個新的質量值重新輸出一個新的BAM文件。
-
Normal 數據對著Tumor的步驟重復一遍。序列比對
-
Normal Reorder
-
Normal Sort
-
Normal MarkDuplicate
- Normal BaseRecalibrator
-
Normal ApplyBQSR
-
Tumor GetPileupSummaries
-
Normal GetPileupSummaries
-
CalculateContamination
-
Mutect2
-
FilterMutectCalls 使用GATK提供的過濾器，過濾得到的突變
-
將過濾后的文件轉換為Annovar注釋所需要的格式
-
使用Annovar對結果注釋
-
使用自己寫的腳本對注釋后的結果過濾，比如按照室間質評要求，過濾掉突變頻率低于1%，測序深度低于500的突變。對GATK某些過濾器過濾掉的結果進行保留和排除，后面使用IGV進行人工篩選。最終輸出的結果為，${sn}.result.SNV.xls（其實是個csv文件，擴展名改為.xls是便于使用excel打開，很多人都這么干）
-
使用Manta Call SV，第一步，生成分析腳本文件runWorkflow.py
-
運行runWorkflow.py獲取SV突變結果
-
使用自己過濾的腳本文件，對Manta獲取的SV過濾。如根據SOMATICSCORE分數過濾，根據hg19_refGene.txt提供文件，計算突變基因等等。
-
使用CnvKIt，獲取CNV突變
-
使用CnvKit畫圖
-
使用py腳本文件，對CnvKit輸出結果過濾。同樣根據hg19_refGene.txt文件匹配基因，以及發生拷貝數變異的區域的外顯子區域等。
-
獲取整體QC數據，insertSize，depth等
- Samtools獲取堿基測序深度
- Samtools flagstat統計比對信息
- 使用py腳本文件對以上獲取的數據做統一處理，獲取整體QC狀態結果。
最終結果：
最終結果及各個命令/任務運行時間如下：
整個分析過程耗時 3h 58min 29s，比較耗時，這個版本為了邏輯清晰一些，多數任務都是串行運行，沒有很好利用計算資源。
后續文章會對整個流程進行優化（主要是并行化），最終分析時間在1h 13min左右（提升約3倍）。
GATK 輸出結果中SNV&INDEL的準確度問題，經過反復試驗，不論如何設置過濾參數，最終的結果始終會有假陰性問題，這是GATK（4.1.3.0）中個別過濾器的問題，目前的補救措施是將部分GATK過濾器過濾掉的結果仍然包含在最終結果中，再使用IGV工具人工過濾（官方文檔也是這么推薦的），判斷該結果是否可靠。
總結
以上是生活随笔為你收集整理的满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【新版本】Aspose.Cells 10
- 下一篇： QQ互联第三方登录jar包