1、研究背景
试剂污染和样品交叉污染是实验室中长期公认的一大问题，污染可能会致使实验结果的不许确。先前研究中曾在TCGA的RNA-seq数据中发现H-HPV18（Hela细胞系衍生的HPV-18病毒），所以本文做者但愿证明TCGA的RNA-seq样本被H-HPV18病毒污染，并肯定其污染来源。web
2、分析流程
3、结果解读
一、开发VirDetect：一种特异性病毒检测软件
为了在RNA-seq数据中检测病毒，做者开发了一个开源软件VirDetect，该软件使用STARV2.4 aligner进行RNA-seq reads的比对，将没有与人类基因组比对上的数据与病毒基因组数据进行比对（图1.A），做者在图1中对VirDetect的检测准确性进行了验证。算法
为了提升测序质量，避免资源浪费，做者将病毒基因组中的人类同源性区域和低复杂度区域进行了屏蔽（mask），做者取滑动窗口长度为75个核苷酸，作滑动窗口算法，以93%的核苷酸类似性做为同源片断的标准，将这些同源片断屏蔽（流程如图1.B所示）。做者对屏蔽效果进行验证，发如今对人类同源性片断和低复杂度片断同时屏蔽后，aligner不会再将reads比对到屏蔽片断上（图1.C）；当突变数<3时，中位敏感率（按读取到病毒基因组的reads比例计算）为99.6%；而当突变数>3时，中位敏感率呈线性降低，但图1.E显示中位PPV（阳性预测率，按读取到正确基因组的比例计算）值为97%，即便肿瘤突变负担较高，利用VirDetect将病毒reads读取到正确的基因组的几率仍然很高。数据库
二、利用VirDetect查找TCGA中的病毒污染
首先做者分析了北卡罗来纳大学对TCGA样本的测序结果，绘制了多种病毒在TCGA28种肿瘤中的表达率热图（图2）。图中的大部分结果和预期一致，如HBV（乙肝病毒）在肝癌中广泛存在，HPV16在头颈部鳞状细胞癌中广泛存在等等。但与预期不符的是：做者利用VirDetect在非子宫颈部肿瘤中检测到了HPV-18病毒（每种肿瘤样本的具体检测数值如图3.A所示），其中比较突出的有透明肾细胞癌（KIRC）、肺鳞状细胞癌（LUSC）等。做者认为能检测到这种数量级的病毒序列，最多是由试剂污染引发的，所以做者但愿在后续的研究中继续寻找污染的源头。express
图2.多种病毒在TCGA28种肿瘤中的表达率热图微信
先前有研究经过Hela细胞系中的23个特异性SNPs（单核苷酸多态性）证实了非子宫颈肿瘤的HPV-18是来源于Hela细胞的。做者对这一结果加以验证，图3.B结果说明除去CSEC（宫颈癌），和3个BLCA（膀胱癌）样本，其他全部17个非子宫颈癌样本的SNPs与Hela细胞的HPV-18基因组SNPs都彻底匹配，验证了非子宫颈肿瘤的HPV-18来源于Hela细胞系app
除了HPV-18外，在96个TCGA样本中还检测到了XMV43病毒（如图3.A下半部分所示），在HPV-18表达量最高的KIRC样本中，XMV43表达量一样最高（图3.C），两者的Spearman相关系数为0.44（p=0.006），同时做者在图3.D中统计了每种肿瘤样本中同时存在两种病毒的样本数，发现他们在BRCA、HNSC、KIRC、LUSC样本中均有明显的共表达。所以，做者认为TCGA样本的污染物中可能同时含有HPV-18和XMV43两种病毒的RNAide
图3.C-D HPV-18和XMV43的共表达状况ui
三、调查污染来源
为了调查形成样本污染的缘由，做者对UHRR（人类通用参考RNA，Universal Human Reference RNA）和大多数TCGA样本同时进行测序，并监测文库的构建（图4.A、C）。大部分UHRR样本的测序从2010年开始，当时的样本中并无检测到HPV-18或XMV43的序列。如图所示，发现两种病毒污染的高峰期发生在2011年4-6月左右，做者认为有可能发生污染的样本并不必定都是与其余UHRR（+）样本在同一天开始构建文库的，所以不能经过图4.A、C的结果直接判断病毒的来源是交叉污染。url
为了确保乳腺癌基因被加入到UHRR样本中，做者将实验室库存的MCF-7和ME16C两株乳腺癌细胞系-加入到了UHRR测序样本中（命名为UHRR+），在图4.B中做者利用箱线图对比了UHRR和UHRR+队列中两种病毒序列的计数，发现两者的HPV-18水平都很高，但UHRR+队列的XMV43水平明显高于UHRR，这代表添加到UHRR+中区别于UHRR的两个细胞株中的一株（MCF-7和ME16C）引入了XMV43病毒污染。spa
为了判断XMV43污染是由上述两株细胞中的哪一株引入的，做者利用实验室2013年准备的MCF-7和ME16C细胞株的RNA-seq数据检测两者的XMV43序列含量，发现MCF-7中的XMV43序列丰度（1%）高于ME16C（0.001%）；且整个XMV43基因组在MCF-7中的覆盖率大于10X，而只有40%的XMV43基因组在ME16C中的覆盖率大于10X（图5.A）。在MCF-7细胞株中，XMV43的序列同源性更高，且所有基因组覆盖率大于10X，所以TCGA样本中的XMV43污染更可能来源于MCF-7。
pBABE-puro hTERT载体是用来转导ME16C细胞株的，其中包含MMLV序列，它没有与人类基因组显著类似的序列和低复杂度区域，所以没有被VirDetect屏蔽。做者发现ME16C转录本覆盖了普罗霉素抗性基因和SV40启动子（两者均存在于pBABE-puro hTERT载体中）。正是该载体在ME16C细胞株转导中的使用，解释了在ME16C的RNA-seq数据中检测到XMV43的缘由（图5.B）
基于以上研究，做者给出了TCGA的RNA-seq样本受到污染的整个流程：一、XMV43病毒感染MCF-7细胞株；二、将MCF-7和ME16C细胞株的RNA添加到包含Hela细胞的UHRR中；三、UHRR与TCGA同时测序，致使了TCGA样本的污染。
图5.比较XMV43序列在两细胞株中的覆盖率
四、补充一种污染物：狂犬病病毒
在10种肿瘤的共计19个样本中检测到了狂犬病病毒，而且集中发生在2012.11-2013.3这段时间内（图6）。尽管病毒的表达量极低，可是在临近的一段时间内被检测到说明这一现象多是由该时间段的试剂污染引发的。
图6.狂犬病病毒在样本中的发现时间
小结
本篇文章中做者揭示了TCGA的RNA-seq数据库污染来源，为避免后续科研中出现数据偏差提供了帮助。在本项研究中，做者开发了能够在RNA-seq数据中检测病毒序列的开源软件-Virdetect，利用该软件检测RNA-seq数据中的病毒表达，证明了TCGA中污染的HPV-18来自Hela细胞，XMV43主要来自MCF-7细胞株，并由此获得了RNA-seq数据污染的流程：一、XMV43病毒感染MCF-7细胞株；二、将MCF-7和ME16C细胞株的RNA添加到包含Hela细胞的UHRR中；三、UHRR与TCGA同时测序，致使了TCGA样本的污染。
点击阅读原文，便可得到小编本次与你们分享的文章。本期的分享就到这里啦，下一期再见~
科研菌学术讨论群，在群内能够用本身的昵称，广告一概踢；其余公众号的宣传也不发，就算是要发，提早和小编商量和确认，否则也是一概踢哈。
欢迎添加小编微信↑↑↑
请你们加个人时候就备注好“学术讨论群”以及本身的“单位+专业+姓名”
编辑：麻辣烫
校审：糯米饭
本文分享自微信公众号 - 科研菌（keyanjun2020）。
若有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一块儿分享。