肿瘤外显子数据处理系列教程（九）拷贝数变异分析（不同软件的比较）
写在前面
上一节我们用 GATK 的 CNV 流程分析了拷贝数变异，系列教程见：
GISTIC2：ftp://ftp.broadinstitute.org/pub/GISTIC2.0/README.txt
CNVkit：https://cnvkit.readthedocs.io/en/stable/pipeline.html
FACETS ：https://github.com/mskcc/facets
CNVkit
CNVkit 的用法比较简单，官方文档写得很详细，这里我不再说明。下面是我用到的脚本（在这里我比较关心的是最后拿到的 segment 文件，即 final.seg，该文件可以用直接载入到 igv 进行可视化，然后和上一节 GATK 的 CNV 结果进行比较）
## cnvkit.sh
GENOME=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
bed=/home/llwu/reference/index/human/CCDS/GRCh38.bed
cnvkit.py batch ../5.gatk/*[ABC]*bqsr.bam ../5.gatk/*techrep_2_bqsr.bam \
--normal ../5.gatk/*germline_bqsr.bam \
--targets ${bed} \
--fasta $GENOME \
--drop-low-coverage --scatter --diagram --method amplicon \
--output-reference my_reference.cnn --output-dir ./
cnvkit.py export seg *bqsr.cns -o gistic.segments
sed 's/_bqsr//' gistic.segments
awk '{print FILENAME"\t"$0}' *bqsr.cns | grep -v chromosome |sed 's/_bqsr.cns//g' |awk '{print $1"\t"$2"\t"$3"\t"$4"\t"$8"\t"$6}' >final.seg
下面两张图是 CNVkit 结果和上一节 GATK CNV 流程的结果，对比两个图，发现这两个软件拿到的结果基本一样，同样的 case4 技术重复的两个样本拷贝数缺失事件碎片化。
GISTIC2
GISTIC2 的安装比较复杂，需要安装 MATLAB，具体可以参考曾老师的博客：用 GISTIC 多个 segment 文件来找SCNA 变异 ：http://www.bio-info-trainee.com/1648.html
# 下载解压
mkdir GISTIC2
cd GISTIC2/
wget -c ftp://ftp.broadinstitute.org/pub/GISTIC2.0/GISTIC_2_0_23.tar.gz
tar zvxf GISTIC_2_0_23.tar.gz
# 安装MATLAB
cd MCR_Installer/
unzip MCRInstaller.zip
./install -mode silent -agreeToLicense yes -destinationFolder ~/software/GISTIC2/MATLAB_Compiler_Runtime/
# 添加环境变量
echo "export XAPPLRESDIR=/home/wsx/biosoft/GISTIC2/MATLAB_Compiler_Runtime/v83/X11/app-defaults:\$XAPPLRESDIR" >> ~/.bashrc
source ~/.bashrc
安装完成之后，新建一个目录，把上一步 cnvkit 拿到的 gistic.segments 文件拷贝到当前目录
mkdir SRP070662_gistic_results
cd SRP070662_gistic_results
cp /trainee/llwu/SRP070662/9.cnvkit/gistic.segments ./
然后把脚本 run_gistic_example 修改了一下，实际运行GISTIC脚本 my.gistic.sh 如下：
!/bin/sh
# run example GISTIC analysis
# output directory
echo --- creating output directory ---
basedir=`pwd`/SRP070662_gistic_results
mkdir -p $basedir
echo --- running GISTIC ---
# input file definitions
segfile=`pwd`/SRP070662_gistic_results/gistic.segments
markersfile=`pwd`/examplefiles/markersfile.txt
refgenefile=`pwd`/refgenefiles/hg38.UCSC.add_miR.160920.refgene.mat
alf=`pwd`/examplefiles/arraylistfile.txt
cnvfile=`pwd`/examplefiles/cnvfile.txt
# call script that sets MCR environment and calls GISTIC executable
./gistic2 -b $basedir -seg $segfile -refgene $refgenefile -genegistic 1 -smallmem 1 -broad 1 -brlen 0.5 -conf 0.90 -armpeel 1 -savegene 1 -gcm extreme
运行很快，十分钟左右就拿到结果，生成的文件有：
all_data_by_genes.txt
all_lesions.conf_90.txt
all_thresholded.by_genes.txt
amp_genes.conf_90.txt
amp_qplot.pdf
amp_qplot.png
broad_data_by_genes.txt
broad_significance_results.txt
broad_values_by_arm.txt
D.cap1.5.mat
del_genes.conf_90.txt
del_qplot.pdf
del_qplot.png
focal_dat.0.5.mat
focal_data_by_genes.txt
freqarms_vs_ngenes.pdf
gistic_inputs.mat
peak_regs.mat
perm_ads.mat
raw_copy_number.pdf
raw_copy_number.png
regions_track.conf_90.bed
sample_cutoffs.txt
sample_seg_counts.txt
scores.0.5.mat
scores.gistic
wide_peak_regs.mat
可视化结果生成的3张图，3 张图均未展示性染色体的拷贝数变异信息，图 1 垂直方向是基因组，水平方向是样本，如果倒过来，就和前面的 cnvkit 结果基本一致。图 2,3 表示的是拷贝数增加或者缺失的情况，水平方向下面是q值，绿线代表显着性阈值（q值 = 0.25）。图上面是 G-score，某个区域的 G-score 分数越高，该区域的发生拷贝数变异事件的可能性就越大，偶然性越低。可以在输出文件 amp_genes.conf_90.txt和del_genes.conf_90.txt查看每一个区域所包含的基因。
在这里，拷贝数缺失事件明显比拷贝数增加事件要多，可能是 case4 技术重复异常带来的影响，我试着删除掉 case4 技术重复的两个样本，然后对比一下结果（如下）：
facet
facet是一个 R 包，但不在 CRAN 中收录，安装的时候有点麻烦，在 Mac 和 Linux 上安装没问题，在 windows 上安装没成功。facet 的输入文件是处理后的 vcf 文件，在这里我用的是 GATK 的HaplotypeCaller工具生成的结果，并走了gvcf流程拿到的 merge.vcf，最后载入 R 中进行分析， R 语言分析部分的代码如下：
if(F){
install.packages('devtools')
devtools::install_github("mskcc/facets", build_vignettes = TRUE)
devtools::install_github("mskcc/pctGCdata")
}
options( stringsAsFactors = F )
library("facets")
# check if .Random.seed exists
seedexists <- exists(".Random.seed")
# save seed
if(seedexists)
oldSeed <- .Random.seed
# Alway use the same random seed
set.seed(0xfade)
# read the data
if(F){
library(vcfR)
vcf_file='merge.vcf'
### 直接读取群体gvcf文件即可
vcf <- read.vcfR( vcf_file, verbose = FALSE )
save(vcf,file = 'input_vcf.Rdata')
}
load(file = 'input_vcf.Rdata')
vcf@fix[1:4,1:4]
vcf@gt[1:14,1:4]
colnames(vcf@gt)
library(stringr)
tmp_f <- function(x){
#x=vcf@gt[,2]
gt=str_split(x,':',simplify = T)[,1]
ad=str_split(str_split(x,':',simplify = T)[,2],',',simplify = T)[,2]
dp=str_split(x,':',simplify = T)[,3]
return(data.frame(gt=gt,dp=as.numeric(dp),ad=as.numeric(ad)))
}
colnames(vcf@gt)
tms=colnames(vcf@gt)[!grepl('_germline',colnames(vcf@gt))][-1]
tmp <- lapply(tms, function(tm){
#tm=tms[1]
print(tm)
nm=paste0(strsplit(tm,'_')[[1]][1],'_germline')
print(nm)
if(nm %in% colnames(vcf@gt)){
snpm=cbind(vcf@fix[,1:2],
tmp_f(vcf@gt[,nm]),
tmp_f(vcf@gt[,tm]))
## only keep tumor or normal have mutations
kp=snpm[,3] %in% c('1/1','0/1') | snpm[,6] %in% c('1/1','0/1')
table(kp)
snpm=snpm[kp,]
## remove those show ./. positions
kp=snpm[,3] == './.' | snpm[,6]== './.'
print(table(!kp))
snpm=snpm[!kp,c(1,2,4,5,7,8)]
rcmat=snpm
rcmat$POS=as.numeric(rcmat$POS)
dim(rcmat)
rcmat=na.omit(rcmat)
colnames(rcmat)=c("Chromosome", "Position","NOR.DP","NOR.RD","TUM.DP","TUM.RD")
rcmat[,1]=gsub('chr','',rcmat$Chrom)
## fit segmentation tree
xx = preProcSample(rcmat)
## estimate allele specific copy numbers
oo=procSample(xx,cval=150)
oo$dipLogR
## EM fit version 1
fit=emcncf(oo)
tmp=fit$cncf
write.table(tmp,file = paste0('facets_cnv_',tm,'.seg.txt'),
row.names = F,col.names = F,quote = F)
head(fit$cncf)
fit$purity
fit$ploidy
png(paste0('facets_cnv_',tm,'.png'),res=150,width = 880,height = 880)
plotSample(x=oo,emfit=fit, sname=tm)
dev.off()
if(F){
fit2=emcncf2(oo)
head(fit2$cncf)
fit2$purity
fit2$ploidy
png(paste0('facets_cnv2_',tm,'.png'),res=150,width = 880,height = 880)
plotSample(x=oo,emfit=fit2, sname=tm)
dev.off()
}
return(c( fit$dipLogR,fit$ploidy,fit$purity))
}
})
tmp <- do.call(rbind,tmp)
rownames(tmp)=tms
colnames(tmp)=c('dipLogR', 'ploidy', 'purity')
write.csv(tmp,'ploidy_and_purity.csv')
我们就拿case1_biorep_A_techrep样本的结果进行展示，主要关心图的第三个 panel ，绘制每个 segment 的主要（黑色 = 2 ）和次要（红色 = 1）拷贝数。底部栏显示相关的细胞分数（cf），深蓝色表示高 cf，浅蓝色表示低 cf，米色表示正常。这个结果就和上面略有不同了，比如图中 19 号染色体在第三个 panel（第一、二个panel 好像还是正常的）中显示是拷贝数增加，但在上面几个软件的结果中 19 号染色体的拷贝数均没有太大的变化。毕竟 facet 的输入数据是 germline mutations，和前面几个软件都不太一样。
可能是软件基于不同类型的输入数据，算法不一样，facet 对 case4 技术重复的分析结果就比较正常了
系列教程写到这里，我们就肿瘤外显子数据分析了 Somatic SNV 和 CNV了，但是还有很多高级分析等着大家，我们下周再见吧。如果你发现教程存在错漏或者不妥的信息，请邮件联系我：Nickier0510@163.com
最后友情宣传一下生信技能树：
如果你看上面的教程有困难，也可以看看B站的WES视频哈