正常细胞中的DNA含量应该比较稳定。肿瘤细胞中DNA含量会增多，恶性肿瘤细胞表现为非整倍体异常增多。测定肿瘤细胞中DNA含量能直接反应细胞群的增值能力和恶性程度。肿瘤基因组研究通常要求肿瘤DNA样本中癌细胞的比例要大于80%，但是实际操作中有时候很难达到这样的要求。肿瘤样本中癌细胞总是混合一定未知比例的正常细胞，我们称肿瘤样本中癌细胞所占的比例为肿瘤纯度(Tumor Purity),称有染色体结构和数目异常导致的肿瘤样本中癌细胞的真正含量为倍性(Tumor ploidy)。估计肿瘤的纯度和倍性有利于癌症基因组进化和肿瘤内异质性研究。

Sequenza，一个R包，能够有效地估计肿瘤细胞的纯度和倍性，以及拷贝数的产生，杂合性的丢失LOH和突变频率谱。

原理是根据肿瘤样本基因组的拷贝数和体细胞突变的等位基因频率来计算纯度和倍性。

安装，需要安装两个工具，一个R package，一个python 包处理前期数据。

setRepositories(graphics = FALSE, ind = 1:6)
install.packages("sequenza")

python 工具：

sudo pip3 install sequenza-utils

Running sequenza

先准备三文件：

tumor bam
normal bam
reference fasta，这里用hg38

第一步：

1.1 一劳永逸

sequenza-utils gc_wiggle -w 50 --fasta Homo_sapiens_assembly38.fasta -o hg38.gc50Base.wig.gz

1.2 将T-N pair bam转换为seqz文件：

sequenza-utils bam2seqz -n Normal.bam -t Tumor.bam --fasta Homo_sapiens_assembly38.fasta -gc hg38.gc50Base.wig.gz -o sample.seqz.gz

1.3 seqz_binning：

sequenza-utils seqz_binning --seqz sample.seqz.gz -w 50 -o small.seqz.gz

第二步，在R中

三个功能：

2.1：chromosome.list 后面说，

test <- sequenza.extract("RA14_TU_CL1-P0-gDNA.seqz.gz", verbose = T, chromosome.list=c(paste0("chr", c(1:22, "X", "Y"))))

2.2：CP

CP <- sequenza.fit(test)

2.3：结果会生成在out.dir目录中。

sequenza.results(sequenza.extract = test,  cp.table = CP, sample.id = "Test",  out.dir="TEST")

3.画图：

3.1

cp.plot(CP)cp.plot.contours(CP, add = TRUE,   likThresh = c(0.999, 0.95),   col = c("lightsalmon", "red"), pch = 20)

3.2:Chromosome view

chromosome.view(mut.tab = test$mutations[[1]], baf.windows = test$BAF[[1]],
                ratio.windows = test$ratio[[1]], min.N.ratio = 1,
                segments = test$segments[[1]],
                main = test$chromosomes[1],
                cellularity = 0.89, ploidy = 1.9,
                avg.depth.ratio = 1)

2022-02-22

报错：

image.png

设置verbose=T，让它输出详细信息，看看是哪里报错了。

test <- sequenza.extract("small.seqz.gz", verbose = T)

image.png

看到是ChrM有错，官方的回答是：“Unfortunately chromosome M gives problems in the extraction step. Probably because of the lack of heterozygous positions.”

那就去掉chrM，再试：

test <- sequenza.extract("small.seqz.gz", verbose = T, chromosome.list=c(paste0("chr", c(1:22, "X", "Y"))))

其他参数：

parallel 设置并行线程数