Somatic CNV 分析工具目前已经有多种了。目前常见有 GATK4 CNV、FACETS、Sequenza，ASCAT等。使用GISTIC工具可以寻找多个样品中，发生CNV的Fragments。

GISTIC默认输出all_lesions.conf_*.txt、amp_genes.conf_*.txt和del_genes.conf_*.txt等文件。这些文件可以被maftools读取并进行后续分析。然而，无数人被误导，以为导入这批数据后，得到的基因既是CNV基因。我自己在分析一批突变数据后，跟cBioportal上TCGA结果进行比较，最终发现先前对GISTIC输出的结果理解有误，并找到了GISTIC CNV基因汇总结果的正确打开方式。

cBioportal是这样展示CNA基因的

在cBioportal中，落入“all_lesions” 区间的基因，被视为“Driven CNA Genes”。例如在肺鳞癌的TCGA数据中：

注意基因右侧加圈的G，这些基因就说GISTIC的到的拷贝数变异显著区间内的基因。然而，若用这些区间内的基因并不能直接视为CNV变化的全部基因。例如，在该数据集中，肺鳞癌常见的PIK3CA扩增，就没有在CNV Peak中，但PIK3CA扩增在该组数据集中出现频率颇高。

cBioportal处理TCGA数据时，是这样判定CNV结果的：

For TCGA studies, the table in allthresholded.bygenes.txt (which is the part of the GISTIC output that is used to determine the copy-number status of each gene in each sample in cBioPortal) is obtained by applying both low- and high-level thresholds to to the gene copy levels of all the samples.

让GISTIC输出CNA基因列表

然而，GISTIC默认情况下并不会输出这个结果！我们需要添加额外的参数-savegene 1，让需要的结果输出出来，例如下面的代码：

gistic2 -b OutDir -seg Input.seg -refgene hg19.mat -conf 0.9 -savegene 1

上述命令运行完，在输出文件夹OutDir中，相对默认参数，会多出4个_data_by_genes.txt的文本。其中all_thresholded.by_genes.txt就是我们需要的结果。而cBioportal的说明也讲的很明确：

• -2 or Deep Deletion indicates a deep loss, possibly a homozygous deletion
• -1 or Shallow Deletion indicates a shallow loss, possibley a heterozygous deletion
• 0 is diploid
• 1 or Gain indicates a low-level gain (a few additional copies, often broad)
• 2 or Amplification indicate a high-level amplification (more copies, often focal)

而核查TCGA CNV的原始数据和cBioportal上展示的CNV结果，cBioportal仅仅考虑了2和-2两种情况。

将CNA基因汇总结果导入maftools

拿到GISTIC的拷贝数变异的数据后，使用下面的示例代码，将拷贝数汇总结果导入maftools：

library(maftools)

CNA = readr::read_delim("OutDir/all_thresholded.by_genes.txt","\t")
Genes = CNA$`Gene Symbol`
CNA = CNA[,-2];CNA = CNA[,-2] # GISTIC输出的CNV表中会包含GI号和基因所在的cytoband，这两列不需要
CNA = reshape2::melt(CNA,id=c('Gene Symbol'),value.name = 'CN')
colnames(CNA) = c('Gene','Sample_name','CN')
CNA = CNA[CNA$CN != 0,] 
CNA = CNA[CNA$CN != 1,] #  low-level gain 和 shallow loss的情况，酌情考虑
CNA = CNA[CNA$CN != -1,] # 我们的例子中仅考虑2与-2的情况
CN = rep('Amp',nrow(CNA))
CN[CNA$CN<0]='Del'
CNA$CN = CN

laml = read.maf(maf = "Input.maf",
                clinicalData = "Input.clinic.txt",
                cnTable = CNA)

oncoplot(maf = laml)

此时，导入maftools的突变结果将与cBioportal网站上一致了。最终得到的Oncoprint如下所示：

• cBioportal说明文档：https://docs.cbioportal.org/1.-general/faq