不合适的人类参考基因组使mapping质量为0

之前抄着gatk的Best Practices做了个somatic的流程,用的是gatk这个页面发布的ucsc.hg19.fasta(md5sum a244d8a32473650b25c6e8e1654387d6)。
在测试某赢给的标准品时,发现有一个经他们wes检验的位点,我这边无法call出。其实对于gatk这些疑似bug,我也见怪不怪了,之前查过gatk论坛,发现这情况不在少数,官方维护者似乎也给不出明确的答复,我倒是遇见过一些由于soft cliping原因出不来的位点(但这种位点是否要保留,我也不太清楚,这里贴个我请教过bcftool的链接),总之gatk用个--alleles强制call,肯定是能出的。如果有想彻底弄明白的朋友,可以看看这两篇帖子
Why do a clear expected variant not show up in the Mutect2 vcf file – GATK
Mutect2 missing a low AF known deletion – GATK
回过头发现还是无法call出,到igv查看此位点,发现mapping质量为0,即匹配到参考基因组多个位置上。

chr6_30856775.png

而在运行mutect2时,在控制台上其实可以看到
mutect2默认过滤.png

显然,mutect2运行时,默认开启的过滤器将其过滤了,关闭mapq相关过滤器-DF MappingQualityAvailableReadFilter -DF MappingQualityNotZeroReadFilter -DF MappingQualityReadFilter
这样操作,结果当然是能凑出来。但是,由此也让人十分好奇这种大面积mapq为0的原因。
由于我igv用的不是太熟,不清楚有没有办法基于igv查是对比到了哪些序列,于是是在igv上复制了其中一条reads,在blat上比对。
比如CTGGGCATGCAGGACCGGACCATCCCAGACAGTGACATCTCTGCTTCCAGCTCCTGGTCAGATTCCACTGCCGGCCGCCACAGCAGGTACTTGGCACACCTGGCACACTTGTAGCTGCCCCGAGAGGAGCTCCTGGGACCTCTACTTCCC
blat结果chr6.png

可以看到,这段序列匹配到了(其它一些十几bp的软剪切就不说了)

>chr6_qbl_hap6
>chr6_mcf_hap5
>chr6_mann_hap4
>chr6_dbb_hap3
>chr6_cox_hap2
>chr6_apd_hap1
>chr6

根据小道消息(https://springwood.me/ngs-somatic-mutation/),只保留常染色体与性染色体,结果是全比上去了,甚至深度还和突变频率还都提高不少。

chr6_30856775_限定染色体前后.png

那么,我们可以简单粗暴地只要常染色体和性染色体吗?我想答案是否定的。
这一点,其实varsome给出了他们的答案Which reference genome is being used to align the reads? (varsome.com)
image.png

在学习gatk此网页后(Human genome reference builds - GRCh38 or hg38 - b37 - hg19 – GATK (broadinstitute.org),也就可以了解其用意了,核心思想总之是,去除具有高度多态性的重叠区域,因为这些区域会影响比对质量。
那为什么不把chr_random和chr_un这些都删掉,反正我设计bed也没打算测它们,或许还能提高运行速度?
这时候就得考虑实际实验时候的脱靶了,chr_random和chr_un都是在人体内真实存在的序列,假设一个不小心,脱靶测到了这段,比对时候各归其位还好说。但倘若所用的参考序列里删了这些chr_random和chr_un,本该比到chr_random和chr_un的序列强行比对上保留的常染色体和性染色体,会有假阳的风险。如果觉得不太明白可以看biostars上的这个问题

到这里,问题也就解决了。后面在ucsc晃悠时,我才发现这个ucsc有对这个问题做出解决(http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/analysisSet/)

image.png

而且liheng大佬在17年就说这个事了(而我在22年还在全网狂搜.....),不过大佬推荐的是千人基因组的(https://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz),ucsc还专门贴了链接,知耻而后勇了可能是。
还有一点奇怪的是,gatk存档的ucsc参考基因组md5为 a244d8a32473650b25c6e8e1654387d6
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz 806c02398f5ac5da8ffd6da2d1d5d1a9
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/latest/hg19.fa.gz 7707462fc100c7d987c075bc146b16ae
和ucsc的两个版本都和其对不上,我个人猜测是初版--迭代版(gatk)--最新版的关系。
另外,云分析平台DNAnexus也有关于参考基因组的说明与推荐
目前看来,我用bwa的话直接下载
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/analysisSet/hg19.p13.plusMT.no_alt_analysis_set.bwa_index.tar.gz
可能是最省事的,或者下载human_g1k_v37与hs37d5。
这边也看到一个比较全面的关于人类参考基因组的帖子
常用或特别的人类fasta参考基因组下载链接 - 知乎 (zhihu.com)

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 196,099评论 5 462
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 82,473评论 2 373
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 143,229评论 0 325
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,570评论 1 267
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,427评论 5 358
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,335评论 1 273
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,737评论 3 386
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,392评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,693评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,730评论 2 312
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,512评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,349评论 3 314
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,750评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,017评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,290评论 1 251
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,706评论 2 342
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,904评论 2 335

推荐阅读更多精彩内容