Control-FREEC检测拷贝数变异

概述

Control-FREEC是一款认可度较高的检测染色体片段拷贝数变异的软件

官方文献为

• Control-FREEC: a tool for assessing copy number and allelic content using next generation sequencing data. V. Boeva, T. Popova, K. Bleakley, P. Chiche, I. Janoueix-Lerosey, O. Delattre and E. Barillot. Bioinformatics, 2012, 28(3):423-5. PMID: 22155870.

  CNA detection part of Control-FREEC (simply FREEC)

• Control-free calling of copy number alterations in deep-sequencing data using GC-content normalization. V. Boeva, A. Zinovyev, K. Bleakley, J.-P. Vert, I. Janoueix-Lerosey, O. Delattre and E. Barillot. Bioinformatics, 2011, 27(2):268-9. PMID: 21081509.
  LOH detection part of Control-FREEC

官方manual网址为

http://boevalab.inf.ethz.ch/FREEC/

原理

（待有时间补充）

安装

软件作者将源码释放在了github，所以可以直接clone

gitclone https://github.com/BoevaLab/FREEC.git

然后进入FREEC文件直接使用，可视化的时候需要安装R包rtracklayer，依赖XML，可以用conda安装，命令为

conda install -c r r-xml=3.98_1.5

此外还需要下载软件自带数据库mapping标签http://xfer.curie.fr/get/vyIi4w8EONl/out100m2_hg38.zip

SNP数据http://xfer.curie.fr/get/zVAhyLdWmfo/hg19_snp131.SingleDiNucl.1based.txt

配置文件

WGS

[general]
BedGraphOutput=TRUE
bedtools=/usr/bin/bedtools
breakPointThreshold=0.9
breakPointType=2
chrFiles=hg19chromFa#分染色体的参考基因组，可以用gatk配套的
chrLenFile=Data/chr_hg19.len#染色体长度数据，可以用软件自带的Data下面数据
coefficientOfVariation = 0.01
contamination=0
contaminationAdjustment=FALSE
degree=3
forceGCcontentNormalization=0
GCcontentProfile = output/GC_profile.cnp
gemMappabilityFile =Data/out100m2_hg19.gem#bam比对质量的编码对应质量值的表,软件自带
intercept=1
minCNAlength=1
minMappabilityPerWindow=0.85
minExpectedGC=0.35
maxExpectedGC=0.55
minimalSubclonePresence=20
maxThreads=8
noisyData=FALSE
outputDir=output/
ploidy=2
printNA=TRUE
readCountThreshold=10
sambamba=sambamba#提前安装sambaba
SambambaThreads=8
samtools=samtools#提前安装samtools，若没有加入环境变量，则需要写成绝对路径
sex=XX#性别，男为XY，女为XX，在X染色体del或dup的判定时起作用
step=10000#WGS的话步长设为10000
telocentromeric=50000#WGS默认
uniqueMatch=TRUE
window=1000000#与step一起起作用，wes和区间捕获时需要设置更小

[sample]
mateFile=sample.HQ.bam#样本bam
#mateCopyNumberFile = chr_19.noDup0.pileup.gz_sample.cpn
inputFormat=BAM#SAM, BAM, pileup, bowtie, eland, arachne, psl (BLAT), BED, Eland. All formated except BAM are also accepted GZipped.
mateOrientation = FR
#miniPileup

[control]
#无

[BAF]
SNPfile = testChr19/hg19_snp131.SingleDiNucl.1based.txt#下载软件自带
minimalCoveragePerPosition = 1
makePileup=testChr19/hg19_snp131.SingleDiNucl.1based.bed#下载软件自带的正常人数据
fastaFile=database/hg19/refseq/ucsc.hg19.fasta
minimalQualityPerPosition=0
shiftInQuality=33

[target]
#captureRegions=如果是WGS，捕获区间不需要设置

WES

$cat config.txt 
[general]
BedGraphOutput=TRUE
bedtools=/usr/bin/bedtools
breakPointThreshold=0.9
breakPointType=2
chrFiles=hg19chromFa
chrLenFile=~/Software/Control_FREEC/Data/chr_hg19.len
coefficientOfVariation = 0.01
contamination=0
contaminationAdjustment=FALSE
degree=1
forceGCcontentNormalization=1
GCcontentProfile = output/GC_profile.cnp
gemMappabilityFile =~/Software/Control_FREEC/Data/out100m2_hg19.gem
intercept=1
minCNAlength=3
minMappabilityPerWindow=0.85
minExpectedGC=0.35
maxExpectedGC=0.55
minimalSubclonePresence=20
maxThreads=8
noisyData=TRUE
outputDir=output/
ploidy=2
printNA=FALSE
readCountThreshold=50
sambamba=~/Software/Sambamba/sambamba
SambambaThreads=8
samtools=~/miniconda2/bin/samtools
sex=XX
#step=10000 当样品是区间捕获panel或者全外显子捕获时，step不设置
#telocentromeric=50000 同上
uniqueMatch=TRUE
window=0 #全外或者panel数据设为0，因为需要提供bed，不需要像全基因组那样平均地设置窗口

[sample]
mateFile=2_Alignment/result/sample.HQ.bam
#mateCopyNumberFile = chr_19.noDup0.pileup.gz_sample.cpn
inputFormat=BAM#SAM, BAM, pileup, bowtie, eland, arachne, psl (BLAT), BED, Eland. All formated except BAM are also accepted GZipped.
mateOrientation = FR
#miniPileup

[control]
#临床样本无对照

[BAF]
SNPfile =testChr19/hg19_snp131.SingleDiNucl.1based.txt
minimalCoveragePerPosition = 1
makePileup=testChr19/hg19_snp131.SingleDiNucl.1based.bed
fastaFile=~/database/hg19/refseq/ucsc.hg19.fasta
minimalQualityPerPosition=0
shiftInQuality=33

[target]

captureRegions=~/database/hg19/target_region/SureSelect_Human_All_Exon_V5_Regions.bed

参数解释：

[general]

chrLenFile指定参考物种染色体长度的文件，共三列，第一列为编号，第二列为染色体名字，第四列为染色体长度。需要注意的是，软件只会分析在该文件中出现的染色体区域。http://bioinfo-out.curie.fr/projects/freec/src/hg19.len

samtools指定samtools路径/usr/local/bin/samtools

ploidy指定参考物种染色体组的个数，通常我们都是分析人的CNV,人是二倍体生物，这个参数的值就是2。

breakPointThreshold官方推荐的取值范围是0.01到0.08，数值越小，预测到的CNV越多。

freec通过分析某一区域的测序深度来预测CNV, 对于全基因组数据，根据滑动窗口模型进行分析，window参数指定窗口的大小，step指定步长；对于全外显子数据，计算测序深度时按照exon区域进行计算，所以window设置为0。

当不提供对照样本时，必须设置chrFiles和GCcontentProfile两个参数。

chrFiles参数的值为一个目录，该目录下时每条染色体的fasta格式的序列。

GCcontentProfile参数的值为一个文件，记录了染色体上固定窗口区域内的GC含量，可以用gccount软件生成。

[sample]项：

mateFile参数指定待分析样本的输入文件，通常都是bam格式的，也支持SAM, pileup等其他格式；

inputFormat指定输入文件的格式;

mateOrientation指定测序方向，对于单端测序的数据，对应的值为0；对于illumina 双端测序的reads, 对应的值为FR。如果输入的bam文件是排序之后的bam文件，需要将该参数的值设为0。

[control]

设置对照样本的输入文件，和sample的设置是一样的。

[BAF]

SNPfile指定已知SNP位点的文件，格式为VCF.

wget https://xfer.curie.fr/get/nil/n5jk4aipGAr/hg19_snp142.SingleDiNucl.1based.txt.gz
gunzip  hg19_snp142.SingleDiNucl.1based.txt.gz

[target]

captureRegions参数的值是bed格式的文件，指定捕获的目的区域，共3列，第一列染色体名字，第二列和第三列分别为区域的起始和终止位置

wget https://xfer.curie.fr/get/nil/n5jk4aipGAr/hg19_snp142.SingleDiNucl.1based.txt.gz
gunzip -c hg19_snp142.SingleDiNucl.1based.txt.gz | awk {'printf  ("%s\t%s\t%s\n", $1,$2-1,$2)'} >hg19_snp142.SingleDiNucl.1based.bed

运行脚本

配置config.txt，然后依次运行以下命令。

#统计样本GC含量
~/Software/Control_FREEC/gccount/gccount -conf config.txt
#Freec跑cnv结果
freec -conf config.txt
#对结果的区间文件进行注释
bedtools intersect -a output/sample.HQ.bam_CNVs -b Database/hg19/all_gene.bed -wa   -wb|  bedtools groupby -i - -g 1-7 -c 11 -o collapse >output/SV_result.xls
#对数据结果可视化，第一步，计算显著性
cat ~/Software/Control_FREEC/FREEC-11.5/scripts/assess_significance.R | R --slave --args sample.HQ.bam_CNVs sample.HQ.bam_ratio.txt 
第二步，绘制曼哈顿图
cat ~/Software/Control_FREEC/FREEC-11.5/scripts/makeGraph.R | R --slave --args 2 sample.HQ.bam_ratio.txt sample.HQ.bam_BAF.txt

结果解析

（待补充）

可视化结果

若要取单独的染色体作图，如下图类似，可以通过修改makeGraph.R里的作图函数op <- par(mfrow = c(5,5))改为c(1,1)，设计每页的小图数为1，可将所有染色体单独画出。

如图，10XWGS数据分析结果，可以明显看到染色体21出现了异常的三拷贝，并且在整条染色体都有体现。而其他染色体显示正常。通过细胞实验室辅助，验证以上结论

参考：

http://boevalab.inf.ethz.ch/FREEC/tutorial.html

h.R里的作图函数op <- par(mfrow = c(5,5))改为c(1,1)，设计每页的小图数为1，可将所有染色体单独画出。

[外链图片转存中…(img-A5ZQf8Sr-1593660328999)]

如图，10XWGS数据分析结果，可以明显看到染色体21出现了异常的三拷贝，并且在整条染色体都有体现。而其他染色体显示正常。通过细胞实验室辅助，验证以上结论

参考：

http://boevalab.inf.ethz.ch/FREEC/tutorial.html

http://blog.sciencenet.cn/blog-2609994-1169205.html