【原】数据库里面很多资料被我们忽视掉了！

Virus expression detection reveals RNA-sequencing contamination in TCGA病毒表达检测揭示TCGA中的RNA测序污染

一、研究背景

试剂污染和样品交叉污染是实验室中长期公认的一大问题，污染可能会导致实验结果的不准确。先前研究中曾在TCGA的RNA-seq数据中发现H-HPV18（Hela细胞系衍生的HPV-18病毒），因此本文作者希望证实TCGA的RNA-seq样本被H-HPV18病毒污染，并确定其污染来源。

二、分析流程

三、结果解读

1、开发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读取到正确的基因组的概率仍然很高。

图1.VirDetect工作流程和准确性检测

2、利用VirDetect查找TCGA中的病毒污染

        首先作者分析了北卡罗来纳大学对TCGA样本的测序结果，绘制了多种病毒在TCGA28种肿瘤中的表达率热图（图2）。图中的大部分结果和预期一致，如HBV（乙肝病毒）在肝癌中普遍存在，HPV16在头颈部鳞状细胞癌中普遍存在等等。但与预期不符的是：作者利用VirDetect在非子宫颈部肿瘤中检测到了HPV-18病毒（每种肿瘤样本的具体检测数值如图3.A所示），其中比较突出的有透明肾细胞癌（KIRC）、肺鳞状细胞癌（LUSC）等。作者认为能检测到这种数量级的病毒序列，最可能是由试剂污染引起的，因此作者希望在后续的研究中继续寻找污染的源头。

图2.多种病毒在TCGA28种肿瘤中的表达率热图

图3.A HPV18和XMV43在TCGA肿瘤队列中的检测水平

        先前有研究通过Hela细胞系中的23个特异性SNPs（单核苷酸多态性）证明了非子宫颈肿瘤的HPV-18是来源于Hela细胞的。作者对这一结果加以验证，图3.B结果说明除去CSEC（宫颈癌），和3个BLCA（膀胱癌）样本，其余所有17个非子宫颈癌样本的SNPs与Hela细胞的HPV-18基因组SNPs都完全匹配，验证了非子宫颈肿瘤的HPV-18来源于Hela细胞系

• 除了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两种病毒的RNA

图3.B 非子宫颈癌样本与H-HPV18基因组SNPs的匹配程度

图3.C-D HPV-18和XMV43的共表达情况

3、调查污染来源

        为了调查造成样本污染的原因，作者对UHRR（人类通用参考RNA，Universal Human Reference RNA）和大多数TCGA样本同时进行测序，并监测文库的构建（图4.A、C）。大部分UHRR样本的测序从2010年开始，当时的样本中并没有检测到HPV-18或XMV43的序列。如图所示，发现两种病毒污染的高峰期发生在2011年4-6月左右，作者认为有可能发生污染的样本并不一定都是与其他UHRR（+）样本在同一天开始构建文库的，因此不能通过图4.A、C的结果直接判断病毒的来源是交叉污染。

• 为了确保乳腺癌基因被加入到UHRR样本中，作者将实验室库存的MCF-7和ME16C两株乳腺癌细胞系-加入到了UHRR测序样本中（命名为UHRR+），在图4.B中作者利用箱线图对比了UHRR和UHRR+队列中两种病毒序列的计数，发现二者的HPV-18水平都很高，但UHRR+队列的XMV43水平明显高于UHRR，这表明添加到UHRR+中区别于UHRR的两个细胞株中的一株（MCF-7和ME16C）引入了XMV43病毒污染。

图4.对UHRR/UHRR的测序及二者的病毒检测结果对比

• 为了判断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样本受到污染的整个流程：1、XMV43病毒感染MCF-7细胞株；2、将MCF-7和ME16C细胞株的RNA添加到包含Hela细胞的UHRR中；3、UHRR与TCGA同时测序，导致了TCGA样本的污染。

图5.比较XMV43序列在两细胞株中的覆盖率

4、补充一种污染物：狂犬病病毒

        在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数据污染的流程：1、XMV43病毒感染MCF-7细胞株；2、将MCF-7和ME16C细胞株的RNA添加到包含Hela细胞的UHRR中；3、UHRR与TCGA同时测序，导致了TCGA样本的污染。