这篇文章的肿瘤外显子数据分析值得批评一下

研究背景

• Lenvatinib 是一种多靶点酪氨酸激酶抑制剂（TKI），可有效靶向血管内皮生长因子受体（VEGFR），并已证明治疗晚期和不可切除的 HCC （uHCC） 的有效结果。

• PD1或PD-L1的免疫检查点阻断疗法已被广泛用于治疗多种类型的人类癌症，具有持久的反应和更高的生存率，包括黑色素瘤和HCC。

文章包括多种组学分析结果，我们重点关注的是基因组学的分析结果。

研究方法

• 患者和样本：人群队列或样本比较复杂，总共有4个队列：Main/ New/ Healthy/ Tissue。

• Main：共有 51 例 HCC 患者接受了 Lenvatinib 和抗 PD1 抗体的联合治疗，26 例为应答者，取样 26 例治疗前血浆，10 例治疗后血浆样本。有 25 名无应答，取样 25 份治疗前血浆样本和 6 份治疗后血浆样本。

• New：13例 uHCC 联合治疗前（基线）血浆样本，8 名有应答和 5 名无应答。

• Healthy：另外有 15 人是健康对照的血浆样本。

• 40-HCC ：采集了40份HCC患者组织样本。该队列包括 20 名无应答（手术切除后再免疫治疗）和 20 名（治疗后再进行手术切除）有应答。这些样本进行了whole genome mutation sequencing 基因组突变测序，Western blot，和蛋白质组学分析。我们重点关注这里的基因组突变测序结果。

• 40-HCC队列的测序策略：WES，使用QuarXeq Human All Exon Probes 3.0 和 QuarHyb One Reagent Kit 试剂盒进行外显子捕获和测序文库构建。使用华大的 DNBSEQ-T7测序仪，双端测序，读长为150 bp，测序深度是 145x~359x

• 数据处理：这里注意一下，方法部分没有指出如何 call somatic 突变，只是简单提到了 GATK ，没有提到对照样本。

研究结果

这里仅仅关心 40-HCC 队列的WES 结果：

• 突变结果：对 uHCC 患者的 40 个癌症组织进行了全外显子组测序 （WES），包括 20 个无应答和 20 个应答：

• fig3：几乎是用 maftools 一键出图。该数据集中 top 突变基因包括PSPC1、ZNF806、KLF18、FER1L5、MUC16、NEFH、PRG4和TP53（图3C）。TP53在52%的样本中发生突变。图3E显示出了应答组的 TMB 显著高于无应答者。且不关心该分析结果是否正确。首先 fig3 这张图，一般来说，文章正文只会放fig3C 和 E，其他子图都无需放到文章正文。

• 差异突变基因组：在应答和无应答两组比较差异突变，fig S2 显示21个基因在应答组中显著富集（p < 0.05），包括ATM、DNAH14、OBSCN、PRAG1、USP9Y、ZNF724、ZNF93、OTOGL、TASOR2、ANAPC1、BAZ2B、CEP350、CTNNA3、FBN2、KIF13B、KMT2B、RIF1、SDK2、UNC79、WNK2和ZNF680

文章评价

文章的 40 HCC 队列 WES 数据，有几个疑点：

• 1. PSPC1 基因在所有样本中均有突变 100%，这是比较罕见的。
  

• 2. 在方法部分没有指出是如何 call 突变或者如何call somatic突变，只是简单提了一下 GATK 和 ANNOVAR 注释，且没有提到对照样本。
  

• 3. top 突变基因除 TP53 外没有其他明显基因，类似 KRAS 、NRAS 、PIK3CA、PTEN、BRAF、ERBB2 等癌症明星突变基因都没有看到。这里可以和 TCGA HCC 队列的突变图谱简单比较一下：Comprehensive and Integrative Genomic Characterization of Hepatocellular Carcinoma
  

• 4. 这 40 HCC 队列采样时间点不一样，无应答（手术切除后再免疫治疗）和 有应答（治疗后再进行手术切除）。但是分析出来的结果，有应当患者（联合治疗后接受手术切除），得到的TMB，比无应当患者（联合治疗前接受手术切除） 还要高

• 5. WES 数据突变分析部分过于简单，fig3 和 fig S2 的子图几乎都是一键出图。在文章中出现这么简单的分析并不常见。

• 6. fig3D 中，应答和无应答患者的 TMB 都较高，和 TCGA 肝癌的 TMB 相比高了很多。

结果重现

以上两个大图都是该文章的 WES 突变图谱，如果想要重现，可以下载文件附件进行重现。这次重现数据来自于文章附件：https://www.cell.com/cms/10.1016/j.celrep.2024.113877/attachment/13e44277-8b84-4b92-8aaf-4eb427b3f6c4/mmc2.xlsx

临床信息

附件中的临床信息，是整理过后的临床三线表，无法获取每一个患者对应的临床信息如 sex/HBV_DNA/HBsAg/OS/Recurrence/Stage，只能获取到患者免疫治疗是否有应答的分组信息

# 清空环境并载入R包
rm




    
(list = ls())
library(maftools)
library(stringr)
library(ggpubr)
library(tidyr)
library(data.table)
library(pheatmap)
library(ggrepel)
library(ggsci)
library(ggplot2)
library(VennDiagram)
library(ggVennDiagram)
clinical = readxl::read_xlsx("mmc2.xlsx", sheet = 7,skip = 3




    
)[,c(1,2)]
clinical = as.data.frame(clinical)
colnames(clinical)[2] = "Response"
write.table(clinical, file = "clinical.txt",sep = "\t",quote = F,row.names = F)

突变结果

文章正文的 fig3 是突变结果，其中：（A）40-HCC 队列中鉴定出的不同类别突变的条形图和箱线图。（B）Ti 和 Tv 图，说明了 HCC 中 6 种 SNV 的分布。（C）Oncoplot 显示 40-HCC 队列的体细胞图谱 （D）40-HCC 队列的有无应答的患者的 TMB 与所有TCGA 队列的比较。（E）40-HCC 队列中有无应答的患者之间的 TMB 比较 （F）突变是否同时或互斥的显著性。（G）40-HCC 队列中 MUC16 和 NEFH 基因共生突变的条形图分析。

annovar = readxl::read_xlsx("mmc2.xlsx", sheet = 4,skip = 2)[,-2]
write.




    
table(annovar, file = "annovar.txt",sep = "\t",quote = F,row.names = F)
annovar_maf = annovarToMaf(annovar = "annovar.txt", 
                           refBuild = 'hg19', 
                           tsbCol = 'Sample code', 
                           table = 'refGene')
annovar_maf = read.maf(annovar_maf, clinicalData = "clinical.txt")
# fig3A
plotmafSummary(maf = annovar_maf)

# fig3B
maf.titv = titv(maf = annovar_maf, plot = FALSE, useSyn =




    
 TRUE)
plotTiTv(res = maf.titv)

# fig3C
oncoplot(maf = annovar_maf,
         clinicalFeatures = "Response",
         sortByAnnotation = T,
         showTumorSampleBarcodes = T,
         barcode_mar = 5,
         gene_mar = 7,
         top = 29
        )

# fig3D
maf_re = subsetMaf(maf = annovar_maf,tsb = clinical[clinical$Response == "responder",




    
1])
maf_no = subsetMaf(maf = annovar_maf,tsb = clinical[clinical$Response != "responder",1])
tcgaCompare(maf = maf_re)
tcgaCompare(maf = maf_no)

# fig3E
tmb = tmb(annovar_maf,captureSize = 50,logScale = T)
tmb = merge(tmb,clinical)
ggplot(data=tmb,mapping = aes(x=Response, y=total,color=Response)) + 
geom_boxplot() + 
stat_boxplot




    
(geom = "errorbar",width=0.15)+
#geom_dotplot(binaxis='y', stackdir='center', dotsize=1,aes(fill = Response))+
geom_jitter(aes(shape=Response), position=position_jitter(0.1),size=5)+
theme_bw()+  
geom_signif(comparisons = list(c("responder", "Nonresponder")),test = t.test,color="black")+
theme(text = element_text(size=18))

# fig3F
somaticInteractions(maf = annovar_maf, top = 25, pvalue = c(0.05, 0.01))

# fig3G




    
png(file = "fig3G.png",width = 900,height = 300)
oncoplot(maf = annovar_maf,
         clinicalFeatures = "Response",
         sortByAnnotation = T,
         showTumorSampleBarcodes = F,
         genes = c("MUC16","NEFH"),
         drawColBar = F,
        )
dev.off()

而附图中的 fig S2 也是突变信息，其中：（A）有无应答突变基因的差异及显著性 （B）有无应答突变基因的比例 （C）有无应答突变基因 ATM 和 OBSCN 的位点分布 （D）有无应答突变基因所在通路的比较 （E）有无应答 TP53 通路的比较

# fig S2A
pt.vs.rt <- mafCompare(m1 = maf_no, m2 = maf_re,




    
 m1Name = 'nonresponder', m2Name = 'responder', minMut = 5)
forestPlot(mafCompareRes = pt.vs.rt, pVal = 0.05)

# fig S2B
genes = c("ATM","DNAH14","OBSCN","PRAG1","USP9Y",
"ZNF724","ZNF93","OTOGL","TASOR2","ANAPC1",
"BAZ2B","CEP350","CTNNA3","FBN2","KIF13B",
"KMT2B","RIF1","SDK2","UNC79","WNK2","ZNF680")
coBarplot(m1 = maf_no, m2 = maf_re, 
          m1Name = 'nonresponder', m2Name =




    
 'responder',
          genes = genes)

# fig S3C
lollipopPlot2(m1 = maf_no, m2 = maf_re, 
              gene = "ATM", 
              AACol1 = "aaChange", AACol2 = "aaChange", 
              m1_label = 'all', m2_label = 'all',labPosAngle = 90,
              m1_name = 'nonresponder', m2_name = 'responder',
              roundedRect = T,showDomainLabel = F)
lollipopPlot2(m1 = maf_no, m2 = maf_re, 
              gene = "OBSCN", 
              AACol1 = "aaChange", AACol2 = "aaChange", 
              m1_label = 'all', m2_label = 'all',labPosAngle = 90,
              m1_name = 'nonresponder', m2_name = 'responder',
              roundedRect = T,showDomainLabel = F)

# fig S4D
pathways