满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）
衛(wèi)計委在2017年，2019年，2020年（還沒有答案）提供標準數(shù)據(jù)用于腫瘤生信分析的室間質(zhì)評。這樣預知結(jié)果的數(shù)據(jù)自然是不能放過了，本文嘗試參考GATK Best Practice：Somatic SNVs + Indels ，Cnvkit，Manta的pipeline來完成滿分流程分析，也可以使用標準數(shù)據(jù)反向判斷GATK Mutect2的實際準確度,算法優(yōu)劣。
注：本文僅用于學習，距離真正的臨床應用還有相當大距離，歡迎大佬批評指正
**1. 分析流程概覽如下：
2. 本文用到的分析系統(tǒng)及分析流程文件
| SliverWorkspace 2.0.277363 | 提供運行控制平臺/個人版 |
| FFPE SNV CNV SV V1.workflow | 分析流程文件，可以一鍵導入分析平臺(點擊查看操作) 當然可以參照圖片中運行腳本，shell里運行，效果也是一樣 |
| 最終結(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 分析報告模板 |
| 分析結(jié)果（pipeline結(jié)果與標準答案） | result.zip pipeline結(jié)果與標準答案 |
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. 詳細流程分析具體如下（不習慣圖形軟件的使用shell變量+腳本運行也是一樣的）
GATK提供的標準流程從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.
這一步省略其實對最終結(jié)果影響不大。
-
Sort
-
MarkDuplicate 標記重復
-
重新校正堿基質(zhì)量值第一步，BaseRecalibrator：計算所有需要重校正的reads和特征值，然后把這些信息輸出為校準表文件
-
重新校正堿基質(zhì)量值第二步，ApplyBQSR：用第一步得到的校準表文件，重新調(diào)整BAM文件中的堿基質(zhì)量值，并使用這個新的質(zhì)量值重新輸出一個新的BAM文件。
-
Normal 數(shù)據(jù)對著Tumor的步驟重復一遍。序列比對
-
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é)果進行保留和排除，后面使用IGV進行人工篩選。最終輸出的結(jié)果為，${sn}.result.SNV.xls（其實是個csv文件，擴展名改為.xls是便于使用excel打開，很多人都這么干）
-
使用Manta Call SV，第一步，生成分析腳本文件runWorkflow.py
-
運行runWorkflow.py獲取SV突變結(jié)果
-
使用自己過濾的腳本文件，對Manta獲取的SV過濾。如根據(jù)SOMATICSCORE分數(shù)過濾，根據(jù)hg19_refGene.txt提供文件，計算突變基因等等。
-
使用CnvKIt，獲取CNV突變
-
使用CnvKit畫圖
-
使用py腳本文件，對CnvKit輸出結(jié)果過濾。同樣根據(jù)hg19_refGene.txt文件匹配基因，以及發(fā)生拷貝數(shù)變異的區(qū)域的外顯子區(qū)域等。
-
獲取整體QC數(shù)據(jù)，insertSize，depth等
- Samtools獲取堿基測序深度
- Samtools flagstat統(tǒng)計比對信息
- 使用py腳本文件對以上獲取的數(shù)據(jù)做統(tǒng)一處理，獲取整體QC狀態(tài)結(jié)果。
最終結(jié)果：
最終結(jié)果及各個命令/任務運行時間如下：
整個分析過程耗時 3h 58min 29s，比較耗時，這個版本為了邏輯清晰一些，多數(shù)任務都是串行運行，沒有很好利用計算資源。
后續(xù)文章會對整個流程進行優(yōu)化（主要是并行化），最終分析時間在1h 13min左右（提升約3倍）。
GATK 輸出結(jié)果中SNV&INDEL的準確度問題，經(jīng)過反復試驗，不論如何設置過濾參數(shù)，最終的結(jié)果始終會有假陰性問題，這是GATK（4.1.3.0）中個別過濾器的問題，目前的補救措施是將部分GATK過濾器過濾掉的結(jié)果仍然包含在最終結(jié)果中，再使用IGV工具人工過濾（官方文檔也是這么推薦的），判斷該結(jié)果是否可靠。
總結(jié)
以上是生活随笔為你收集整理的满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【新版本】Aspose.Cells 10
- 下一篇： QQ互联第三方登录jar包