使用sequenza软件判定肿瘤纯度

软件官方主页是：http://www.cbs.dtu.dk/biotools/sequenza/ 算法描述及使用说明均可以在官网找到。软件结构介绍该软件包括两个部分:a python-based preprocessing toolan R package implementing the model fitting and visualization functionspython脚本部分可以分成2个步骤：First, it calculates the GC content in sliding windows from a genome reference file in FASTA format.Second, it processes the sequencing data from the tumor and normal specimens, which must be in the Pileup format, as output by SAMtools得到的结果就需要导入到R里面进行后续处理！！！R包的流程也可以分成3个步骤：first, sequenza. extract efficiently reads the input file into R, performs GC-content normalization of the tumor versus normal depth ratio, and performs allele-specific segmentation using the 'copynumber’ package.Second, **sequenza.fit** applies the model to infer cellularity and ploidy parameters and copy number profiles.Finally, sequenza.results returns the results of the estimation together with alternative solutions and visualization of the data and the model along the genome and the individual chromosomes.发表该软件的文章当时使用了10 个 ovarian serous carcinomas (OVCA) 和 20 个clear-cell renal cell carcinomas (KIRC) 肿瘤样本做例子。并且与SNP array data analyzed by allele-specific copy number analysis of tumors (ASCAT)检测结果的一致性，还跟两个流行的软件absCN-seq，ABSOLUTE做了对比分析说明自己开发的sequenza工具表现最优！软件使用方法python部分非常好安装：From Pypipip install sequenza-utilsFrom Sourcesgit clone https://bitbucket.org/sequenza_tools/sequenza-utilscd sequenza-utilspython setup.py testpython setup.py install安装完毕即可根据说明文档开始对自己的数据一步步处理：http://sequenza-utils.readthedocs.io/en/latest/ 老实说，这个是我见过最懒的说明书。后来发现我找错了，https://bitbucket.org/sequenza_tools/sequenza/ 其实有比较详尽的说明书。sequenza-utils -hif [  ! -f  hg38.gc50Base.txt.gz ]; thensequenza-utils  gc_wiggle -f $reference  -w 50 -o - | gzip > hg38.gc50Base.txt.gzfisequenza-utils bam2seqz -gc hg38.gc50Base.txt.gz \-F $reference -n $normal_bam -t  $tumor_bam | gzip > ${sample}.seqz.gzsequenza-utils seqz_binning -w 50 -s ${sample}.seqz.gz | gzip > ${sample}_small.seqz.gz得到的文件如下，取决于自己实验的wes数据量182M May  1 17:47 hg38.gc50Base.txt.gz54K May  1 19:36 S12_1_small.seqz.gz39K May  1 23:07 S12_13_small.seqz.gz1.6M May  2 01:03 S12_2_small.seqz.gz9.0M May  2 03:28 S12_22_small.seqz.gz其实可以直接用varscan结果，不走自己的python上游流程利用varscan的somatic mutation calling的结果，导入R，如下；snp <- read.table("varscan.snp", header = TRUE, sep = "\t")cnv <- read.table("varscan.copynumber", header = TRUE, sep = "\t")seqz.data <- VarScan2seqz(varscan.somatic = snp, varscan.copynumber = cnv)write.table(seqz.data, "my.sample.seqz", col.names = TRUE, row.names = FALSE, sep = "\t")制作得到的my.sample.seqz文件即可进入R流程！R 语言部分也是非常好安装，只是一个R包而已Reference manual:sequenza.pdfVignettes:Sequenza user guide这个时候的说明书详尽程度尚可！library(sequenza)#In the package is provided a small seqz filedata.file <-  system.file("data", "example.seqz.txt.gz", package = "sequenza")data.file# 根据前面介绍的，三部曲# sequenza.extract: process seqz data, normalization and segmentationtest <- sequenza.extract(data.file, verbose = FALSE)#sequenza.fit: run grid-search approach to estimate cellularity and ploidyCP <- sequenza.fit(test)#sequenza.results: write files and plots using suggested or selected solutionsequenza.results(sequenza.extract = test,cp.table = CP, sample.id = "Test",out.dir="TEST")## 然后是理解输出结果文件及各种成图真的是不能太方便了，所有的输出结果尽在眼前！！！不过，CP那个步骤可能会比较耗时。结果文件非常多，需要仔细理解

当然了， 我感兴趣的其实只有："cellularity"    "ploidy.estimate"   "ploidy.mean.cn"0.89    2   1.762051146082690.89    2   1.762051146082690.89    2   1.76205114608269一些炫目的图：

可以看到这个软件顺便找了CNV，所以：生信技能树GATK4系列教程GATK4的gvcf流程你以为的可能不是你以为的新鲜出炉的GATK4培训教材全套PPT，赶快下载学习吧曾老湿最新私已：GATK4实战教程GATK4的CNV流程-hg38