基因组结构变异 (structural variation) 检测

检测SVs的意义：

SVs对基因组的影响比起SNP更大，一旦发生往往会给生命体带来重大影响，比如导致出生缺陷、癌症等；

有研究发现，基因组上的SVs比起SNP而言，更能代表人类群体的多样性特征；

稀有且相同的一些结构性变异往往和疾病(包括癌症)的发生相互关联甚至还是其直接的致病诱因。

目前已有的检测SVs的策略：

Read Pair，一般称为Pair-End Mapping，简称RP或者PEM；

Split Read，简称SR；

Read Depth，简称RD，也有人将其称为RC——Read Count的意思，它与Read Depth是同一回事，顾名思义都是利用read覆盖情况来检测变异的方法；

序列从头组装(de novo Assembly， 简称AS)的方法。

(1) 这个插入片段长度的分布是RP方法进行变异检测的一个关键信息

如果插入片段长度有异常，即明显偏离分布中心的片段长度，则它可能是隐含的变异信号！

举个例子，如果我们发现它这个计算出来的插入片段长度与正态分布的中心相比大了200bp(假设这个200bp已经大于3个标准差了)，那么就意味着参考基因组比read1和read2所在的片段要长200bp，通过类似这样的方式，我们就可以发现read1和read2所在的序列片段相比与参考基因组而言发生了200bp的删除(Deletion)。

(2) 通过比对read1和read2之间的序列位置关系，还能够发现更多非线性的序列变异

比如，序列倒置(Inversion)，因为，按照PE的测序原理，read1和read2与参考基因组相比对，正好是一正一负，要么是read1比上正链，read2比上负链，要么是反过来，而且read1和read2都应处于同一个染色体上，如果不是这种现象，那么就很可能是序列的非线性结构性变异所致

RP方法存在的缺陷：

RP对于大长度Deletion(通常是大于1kbp)比较敏感，准确性也高，而50bp-200bp的这个范围内的变异是它的一个检测暗区。

对于Deletion的检测，由于要求插入片段长度的变化要具有统计意义上的显著性，所以它所能检测到的片段长度就会受插入片段长度的标准差(SD)所影响。简单来说，就是越大的序列删除越偏离正常的长度中心，才越容易被检测到

它所能检测的Insertion序列，长度无法超过插入片段的长度

如果这个Insertion序列很长——举个极端的例子——整个插入片段都是Insertion序列，那么你会发现read1和read2根本就不会比对上基因组，它在基因组上一点信号都没有，你甚至都不知道有这个序列的存在。另外，Insertion的检测精度也同时受限于插入片段长度的标准差。

PEM方法共有两种策略来鉴定SVs/CNVs：

clustering approach：使用预先定义的距离(predefined distance)来检测异常(discordant)的PE reads，那这里就产生了一个问题：如何得到合适的预定义距离？

model-based approach：对全基因范围的PE read的插入长度进行概率分布建模，然后基于统计检验得到异常的PE reads；

使用PEM方法检测CNV的工具列表

Tool

URL

Language

Input

Comments

BreakDancer

Perl, C++

Alignment files

Predicting nsertions, deletions, inversions, inter- and intra-chromosomal translocations

PEMer

Perl, Python

FASTA

Using simulation-based error models to call SVs

VariationHunter

C

DIVETa

Detecting insertions, deletions and inversions

commonLAW

C++

Alignment files

Aligning multiple samples simultaneously to gain accurate SVs using maximum parsimony model

GASV

Java

BAM

A geometric approach for classification and comparison of structural variants

Spanner

N/A

N/A

N/A

Using PEM to detect tandem duplications

SR算法的核心也是对非正常PE比对数据的利用

对比上面提到的RP方法：

RP中的非正常比对，通常是read1和read2在距离或者位置关系上存在着不正常的情形，而它的一对PE read都是能够“无伤”地进行比对的；

但SR一般是指这两条PE的read，有一条能够正常比对上参考基因组，但是另一条却不行的情形。

这个时候，比对软件(比如BWA)会尝试把这条没能够正常比上基因组的read在插入片段长度的波动范围内，使用更加宽松的Smith-Waterman局部比对方法，尝试搜索这条read最终可能比对得上的位置。如果这条read有一部分能够比上，那么BWA会对其进行软切除——soft-clip(CIGAR序列中包含S的那些read，图9中也有这类比对情况)，标记能够比上的子序列(但不能比上的序列还是留在原read中，这也是软切除的含义)。

SR的一个优势在于，它所检测到的SVs断点能精确到单个碱基，但是也和大多数的RP方法一样，无法解决复杂结构性变异的情形。

使用Split-Read方法检测CNV的工具列表

Tool

URL

Language

Input

Comments

AGE

C++

FASTA

A dynamic-programming algorithm using optimal alignments with gap excision to detect breakpoints

Pindel

C++

BAM /FASTQ

Using a pattern growth approach to identify breakpoints of various SVs

SLOPE

C++

SAM/FASTQ/MAQb

Locating SVs from targeted sequencing data

SRiC

N/A

N/A

BLAT output

CalibratingSV calling using realistic error models

其实CNV实质上是序列Deletion或Duplication，是可以归类于Deletion和Insertion这个大的分类的，只是由于它的发生有着其独特的特点，而且往往还比较长，所以也就习惯了独立区分。

Read Depth方法一般用于CNVs的检测，其检测原理为：

基因组区域的Read Depth与其拷贝数相关

全基因组测序(WGS)得到的覆盖深度呈现出来的是一个泊松分布——因为基因组上任意一个位点被测到的几率都是很低的——是一个小概率事件，在很大量的测序read条件下，其覆盖就会呈现一个泊松分布，如下图。

拷贝数增加会使得该区域的Read Depth高于期望值，而拷贝数缺失使得该区域的Read Depth低于高于期望值

目前有两种利用Read depth信息检测CNV的策略：

一种是，通过检测样本在参考基因组上read的深度分布情况来发现CNV，这类适用于单样本，也是用的比较多的一个方法；

另一种则是通过识别并比较两个样本在基因组上存在丢失和重复倍增的区域，以此来获得彼此相对的CNV，适用于case-control模型的样本，或者肿瘤样本Somatic CNV的发现，这有点像CGH芯片的方法，因此又被称作CNV-seq，具体信息见 2.1. CNV-seq。

CNVnator使用的是第一种策略，同时也广泛地被用于检测大的CNV，当然还有很多冷门的软件，这里就不再列举了；CNV-seq使用的则是第二种策略。

使用Read-Depth方法检测CNV的工具列表

Tool

URL

Language

Input

Comments

SegSeqa

Matlab

Aligned read positions

Detecting CNV breakpoints using massively parallel sequence data

CNV-seqa

Perl, R

Aligned read positions

Identifying CNVs using the difference of observed copy number ratios

RDXplorerb

Python, Shell

BAM

Detecting CNVs through event-wise testing algorithm on normalized read depth of coverage

BIC-seqa

Perl, R

BAM

Using the Bayesian information criterion to detect CNVs based on uniquely mapped reads

CNAsega

R

BAM

Using flowcell-to-flowcell variability in cancer and control samples to reduce false positives

cn.MOPSb

R

BAM/read count matrices

Modelling of read depths across samples at each genomic position using mixture Poisson model

ReadDepth

R

BED files

Using breakpoints to increase the resolution of CNV detection from low-coverage reads

rSW-seqa

C

Aligned read positions

Identifying CNVs by comparing matched tumor and control sample

CNVnator

C++

BAM

Using mean-shift approach and performing multiple-bandwidth partitioning and GC correction

CNVnorma

R

Aligned read positions

Identifying contamination level with normal cells

CMDSb

C, R

Aligned read positions

Discovering CNVs from multiple samples

mrCaNaVar

C

SAM

A tool to detect large segmental duplications and insertions

CNVeM

N/A

N/A

N/A

Predicting CNV breakpoints in base-pair resolution

cnvHMM

C

Consensus sequence from SAMtools

Using HMM to detect CNV

其实从上面看下来，SVs检测最大的难点实际上是read太短导致的。就因为read太短，我们不能够在比对的时候横跨基因组重复区域；就因为read太短，很多大的Insertion序列根本就没能够看到信息；就因为read太短，比对才那么纠结，我们才需要用各种数学模型来 猜测这个变异到底应该是什么等等。

但从理论上来讲，三代测序和de novo assembly 的方法应该要算是基因组结构性变异检测上最有效的方法，它们都能够检测所有类型的结构性变异。

最早的CNV检测实验方法，是采用染色体核型分析 (karyotyping)和FISH(荧光原味杂交)

在2003年，出现了arrayCGH方法和SNP芯片方法，这两种全基因范围的高通量检测方法，它们的检测原理是：

依据拷贝数与杂交信号的强度正相关，对特定的待检测基因组区段设计特异性杂交探针，将配对的Test Genome和Reference Genome分别用不同颜色的荧光标记，然后与芯片进行杂交，比较两种荧光信号的强弱

比如SNP6.0芯片就是其中的一款代表产品，TCGA里面主要是通过Affymetrix SNP6.0 array这款芯片来测拷贝数变异，值得注意的是，并不是只有TCGA利用了SNP6.这个芯片数据，著名的CCLE计划也对一千多细胞系处理了SNP6.0芯片，数据也是可以下载的。

对SNP6.0的拷贝数芯片来说，通常是用PICNIC等软件处理原始数据，就可以得到的segment记录文件，每个样本一个结果，下面是示例结果：

Chromosome Start End Num_Probes Segment_Mean

1 61735 1510801 226 -0.0397

1 1627918 1672603 17 -0.92

1 1687587 16153497 8176 0.0077

1 16153536 16153925 5 -2.7441

1 16154201 16155010 4 -0.8711

1 16165661 72768498 34630 0.0048

1 72768916 72811148 46 -1.7394

1 72811904 95674710 14901 0.0026

1 95676511 95676518 2 -1.6636

表明了某条染色体的某个区域内，SNP6.0芯片设计了多少个探针，芯片结果的拷贝数值是多少(这个区域的拷贝数用Segment_Mean)。通常二倍体的Segment_Mean值为0，可以用-0.2和0.2来作为该区域是否缺失或者扩增，也有人选择0.4作为阈值。

但是这些基于芯片杂交的检测方法存在一些缺点：

存在杂交信号噪声；

只能对已知的变异进行检测，无法发现新型的变异和罕见变异；

低分辨率；

有限的基因组覆盖度，有些区域因为一些原因没有设计检测探针，因此也就覆盖不到这些区域；

CNV-seq：通过低倍全基因组测序检测CNV技术，2009年由Xie等开发提出

原理：

CNV-seq检测原理派生自aCGH(比较基因组杂交芯片)，但不针对芯片数据，而是使用NGS测序数据进行检测。将等量的待测样本DNA和正常样本DNA分别进行建库进行NGS测序后，与reference sequence进行比对，通过比较待测样本和正常样本每个滑动窗口内reads数目的多少(待测样本reads数/正常对照样本reads数)来确定待测样本每个位点的拷贝数。

其实简单来说，就是前面提到的Read Depth方法

它可以检出全基因组水平的大片段的CNV

检测效果：

采用双盲实验设计，对染色体结构异常的样本进行检测，滑动窗口大小为20kb的情况下，CNV-seq和高密度SNP-array对已知致病CNVs都能达到100%的检出。

与中等密度SNP-array相比，CNV-seq表现更优。

基于PCR-free的WGS文库构建方法，GC偏好性性更小，数据波动更小，表现更稳定。因为没有捕获设计(与WES比较)，reads在全基因组范围覆盖均一性和随机性相比捕获测序的方法更好

CNV-seq与芯片平台灵敏度比较

从上图可以看到CNV-seq的检测灵敏度明显要高于芯片平台，而且以mate-pair方式建库的NGS检测方法的检出率更高

SCNA: somatic copy number alterations，一般指的是癌症中tumor tissue中的拷贝数变异

SCNA检测存在的难度：

进行tumor tissue的bulk取样测序，组织中的细胞异质性高，信号比较弱

可能受到正常的germline细胞的污染

来自父本和母本的reads混在一起

注：

目前还不懂用GATK4检测CNV流程中每一步的原理，只是简单地将生物技能树中的内容复制过来，而且近期也没有对这个分析流程进行深究的打算，仅当做一个知识的补充——毕竟GATK4并不是目前用于call SNV的主流工具

1. 首先制作外显子坐标记录文件

# bed to intervals_list

$ cat exon_probe.hg38.gene.bed|awk '{print "chr"$0}' >hg38.chr.bed

$ java -jar ~/biosoft/picardtools/2.9.2/picard.jar BedToIntervalList \

I=hg38.chr.bed \

O=list.interval_list \

SD=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict

# Preprocess Intervals

$ gatk PreprocessIntervals \

-L list.interval_list \

--sequence-dictionary /home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict \

--reference /home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta \

--padding 250 \

--bin-length 0 \

--interval-merging-rule OVERLAPPING_ONLY \

--output targets.preprocessed.interval.list

2. 然后把bam文件转为外显子reads数

for i in /PATH/TO/*.bam

do

j=$(basename "$i" _recal.bam)

echo $j

## step1 : CollectReadCounts

gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" CollectReadCounts \

-I $i \

-L $tl \

-R $GENOME \

--format HDF5 \

--interval-merging-rule OVERLAPPING_ONLY \

--output $j.clean_counts.hdf5

## step2 : CollectAllelicCounts，这一步非常耗时，而且占空间

gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" CollectAllelicCounts \

-I $i \

-L $tl \

-R $GENOME \

-O $j.allelicCounts.tsv

done

3. 接着合并所有的normal样本的数据创建 cnvponM.pon.hdf5

$ gatk --java-options "-Xmx20g" CreateReadCountPanelOfNormals \

--minimum-interval-median-percentile 5.0 \

--output cnvponM.pon.hdf5 \

--input counts/OSCC_01_N.clean_counts.hdf5 \

--input counts/OSCC_04_N.clean_counts.hdf5

值得注意的是这个cnvponM.pon.hdf5文件， h5py文件是存放两类对象的容器，数据集(dataset)和组(group)，dataset类似数组类的数据集合，和numpy的数组差不多。group是像文件夹一样的容器，它好比python中的字典，有键(key)和值(value)。group中可以存放dataset或者其他的group。”键”就是组成员的名称，”值”就是组成员对象本身(组或者数据集)。

5. 最后走真正的CNV流程

for i in ./counts/*

do

j=$(basename "$i" .clean_counts.hdf5)

gatk --java-options "-Xmx20g" DenoiseReadCounts \

-I $i \

--count-panel-of-normals q \

--standardized-copy-ratios $j.clean.standardizedCR.tsv \

--denoised-copy-ratios $j.clean.denoisedCR.tsv

done

for i in denoisedCR/*

do

j=$(basename "$i" .clean.denoisedCR.tsv)

## ModelSegments的时候有两个策略，是否利用CollectAllelicCounts的结果

gatk --java-options "-Xmx20g" ModelSegments \

--denoised-copy-ratios $i \

--output segments \

--output-prefix $j

gatk --java-options "-Xmx20g" CallCopyRatioSegments \

-I segments/$j.cr.seg \

-O segments/$j.clean.called.seg

CNVnator的CNV检测原理：

切分bins，定量bins的RD信号值；

首先将整个基因组切分成连续而不重叠的bins，每个bins的大小相同。然后计算每个bin的RD信号强度，由于在之前的研究中发现GC含量会影响RD信号的定量

因此根据GC含量对RD信号进行了校正，公式如下：

分割不同CN的片段 (segments)。从全基因组范围的bins的RD信号分布图中，找出不同的CN (Copy Number) 的区域，确定相邻CN区域的breakpoint

该过程用到了图像处理领域常用的一个算法：mean shift 算法，该算法在下文补充部分有进一步的说明，点 这里 查看

对于图中的每一个点(每一个点表示一个bin)获得它的mean-shift向量——指向与它最相似的那些bins(要么朝左要么朝右)；

根据这些bins的mean-shift向量的朝向，确定相邻CN区域的breakpoint——当相邻两个bins的mean-shift向量的朝向为背靠背形式时，breakpoint位于这两个bins之间；

安装CNVnator软件是一种神级坑的体验，它底层需要依赖ROOT这个计算包，所以在安装CNVnator之前你需要将ROOT安装上，并且设置好ROOT的环境变量：

1. 安装配置和启动ROOT

# 1.1 安装ROOT可以从ROOT官网下载源码自己编译，但是如果自己编译八成又会遇到各种稀奇古怪的报错，友情提示还是直接下载官方提供的二进制包吧

$ wget -c https://root.cern.ch/download/root_v6.14.06.Linux-centos7-x86_64-gcc4.8.tar.gz

$ tar zxvf root_v6.14.06.Linux-centos7-x86_64-gcc4.8.tar.gz

# 1.2 然后编辑你的.bashrc或者.bash_profile

export ROOTSYS=/PATH/TO/root

export LD_LIBRARY_PATH=$ROOTSYS/lib:$LD_LIBRARY_PATH

# 1.3 启动ROOT

source root/bin/thisroot.sh

2. 安装CNVnator

# 由于官方提供编译过的二进制安装包，所以只能自己硬着头皮下载源码编译了，伤心

# 2.1 下载CNVnator

$ wget -c https://github.com/abyzovlab/CNVnator/releases/download/v0.3.3/CNVnator_v0.3.3.zip

$ unzip CNVnator_v0.3.3.zip

# 2.2 编译CNVnator源码包中夹带的samtools，这一步很可能报错！我死在了这一步

$ cd src/samtools

$ make

# 2.3 编译CNVnator，这一步也很可能报错！！

$ cd ../

$ make

上面讲了这么多怎么安装CNVnator，然后依据自己的亲身经历告诫各位，强烈不推荐采用上面的方法自己安装，这样很可能你会经历一个从入门到放弃的一个辛酸的过程

强烈推荐用conda安装

# 用以下命令寻找可以通过conda安装CNVnator的镜像

$ anaconda search -t conda cnvnator

# 通过上一步的搜索找到一个符号要求的镜像https://conda.anaconda.org/pwwang

$ conda install --channel https://conda.anaconda.org/pwwang cnvnator

# 或者用以下的简写形式

$ conda install -c pwwang cnvnator

提取mapping信息

$ ./cnvnator [-genome name] -root out.root [-chrom name1 ...] -tree [file1.bam ...]

Example:

cnvnator -root Sample1.root -tree Sample1.sorted.bam -unique

注意：这一步非常耗内存，可以通过指定-chrom参数来只对指点的染色体提取mapping信息，来降低内存压力，但是染色体名必须与SAM/BAM文件的Header信息中的染色体名一致

@HD VN:1.4 GO:none SO:coordinate

@SQ SN:1 LN:249250621

@SQ SN:2 LN:243199373

@SQ SN:3 LN:198022430

生成质量分布图HISTOGRAM

$ ./cnvnator [-genome name] -root file.root [-chrom name1 ...] -his bin_size [-d dir]

Example:

cnvnator -root Sample1.root -his 100 -d hg38/Homo_sapiens_assembly38.fasta

这一步操作内存消耗不高，但也可以通过指定-chrom参数来只对指点的染色体进行操作

需要提供每条染色体的fasta序列，用-d参数指定

生成统计结果

$ ./cnvnator -root file.root [-chrom name1 ...] -stat bin_size

Example:

cnvnator -root Sample1.root -stat 100

RD信息分割partipition

$ ./cnvnator -root file.root [-chrom name1 ...] -partition bin_size [-ngc]

Example:

cnvnator -root Sample1.root -partition 100

可以通过设置-ngc参数关闭GC校正RD信息功能

变异检出

$ ./cnvnator -root file.root [-chrom name1 ...] -call bin_size [-ngc]

Example:

cnvnator -root Sample1.root -call 100 > Sample1.cnvnator.vcf

运行结果默认以标准输出输出STDOUT

输出形式如下：

CNV_type coordinates CNV_size normalized_RD e-val1 e-val2 e-val3 e-val4 q0

normalized_RD -- normalized to 1.

e-val1 -- is calculated using t-test statistics.

e-val2 -- is from the probability of RD values within the region to be in

the tails of a gaussian distribution describing frequencies of RD values in bins.

e-val3 -- same as e-val1 but for the middle of CNV

e-val4 -- same as e-val2 but for the middle of CNV

q0 -- fraction of reads mapped with q0 quality

CNVkit是为control-case配对样本的杂交捕获测序(如WES)的CNV鉴定而开发的

它的检测过程如下图

CNVkit主要通过batch工具完成这个复杂的多步计算，且对于不同的使用情景，有以下几种用法：

# From baits and tumor/normal BAMs

cnvkit.py batch *Tumor.bam --normal *Normal.bam \

--targets my_baits.bed --annotate refFlat.txt \

--fasta hg19.fasta --access data/access-5kb-mappable.hg19.bed \

--output-reference my_reference.cnn --output-dir results/ \

--diagram --scatter

# Reusing a reference for additional samples

cnvkit.py batch *Tumor.bam -r Reference.cnn -d results/

# Reusing targets and antitargets to build a new reference, but no analysis

cnvkit.py batch -n *Normal.bam --output-reference new_reference.cnn \

-t my_targets.bed -a my_antitargets.bed \

-f hg19.fasta -g data/access-5kb-mappable.hg19.bed

其实batch工具是CNVkit内部多个工具命令的组合打包，完整的执行过程如下：

cnvkit.py access hg19.fa -o access.hg19.bed

cnvkit.py autobin *.bam -t baits.bed -g access.hg19.bed [--annotate refFlat.txt --short-names]

# For each sample...

cnvkit.py coverage Sample.bam baits.target.bed -o Sample.targetcoverage.cnn

cnvkit.py coverage Sample.bam baits.antitarget.bed -o Sample.antitargetcoverage.cnn

# With all normal samples...

cnvkit.py reference *Normal.{,anti}targetcoverage.cnn --fasta hg19.fa -o my_reference.cnn

# For each tumor sample...

cnvkit.py fix Sample.targetcoverage.cnn Sample.antitargetcoverage.cnn my_reference.cnn -o Sample.cnr

cnvkit.py segment Sample.cnr -o Sample.cns

# Optionally, with --scatter and --diagram

cnvkit.py scatter Sample.cnr -s Sample.cns -o Sample-scatter.pdf

cnvkit.py diagram Sample.cnr -s Sample.cns -o Sample-diagram.pdf

(1)虽然CNVkit是为捕获测序开发的但是也可以用于WGS的分析，可以将WGS理解为对全基因组的捕获测序

可以通过--method wgs参数，设置为WGS模式；

一旦设为WGS模式，整个基因组就成了sequencing-accessible regions(通过--access指定)，那么对于捕获测序需要提供的目标区域BED文件(通过-t参数指定)，此时会使用--access指定的区域，因此就不需要再提供了；

需要提供注释文件，通过--annotate参数指定；

$ cnvkit.py batch Sample1.bam Sample2.bam -n Control1.bam Control2.bam \

-m wgs -f hg19.fasta --annotate refFlat.txt

为了提高WGS检测CNV的准确性，开发者提供了以下几点建议：

(2)当没有配对样本时，如何解决？

官方提供的解决方案是，通过在--normal/-n选项后不添加任何参数(即不知道BAM文件)来构造一个 “flat” reference：在这个人为构造的reference中，所有的bins的coverage相同

$ cnvkit.py batch *Tumor.bam -n -t my_baits.bed -f hg19.fasta \

--access data/access-5kb-mappable.hg19.bed \

--output-reference my_flat_reference.cnn -d example2/

(3)若要处理的情况是WGS且为非配对样本，即是上面提到的两个情况的组合，此时怎么解决？

$ cnvkit.py batch Sample.bam -n -m wgs -f hg19.fasta --annotate refFlat.txt

$ lumpyexpress \

-B Sample1.sorted.bam \

-S Sample1.discordants.sorted.bam \

-D Sample1.splitters.sorted.bam \

-o Sample1.lumpu.sv.vcf

# SV检测

$ delly call \

-g hg38/Homo_sapiens_assembly38.fasta \

-o Sample1.delly.sv.bcf \

-n Sample1.sorted.bam

# 过滤结果

$ delly filter \

-f germline \

-p \

-q 20 Sample1.delly.sv.bcf \

-o Sample1.delly.sv.filter.bcf

大多数的reads都能在基因组中找到它唯一的mapping位置，但是如果这条reads是来源于重复序列区域，那么它就会得到multiple mapping的位置，且每条mapping的结果的分值都相同，此时对于这些multiple mapping的reads如何处置呢？

目前有两种主要的处理办法：

一种处理方式是在mapping之前，将reference中的重复序列过滤(mask)；

另一种方法是，从这些得分相同的multiple mapping位置中随机挑选一个，最为这条reads最后的归属；

对于paired-end的测序数据，还可以利用额外的信息来辅助multiple mapping情况的处置：

利用PE reads mapping的插入片段距离和双端的mapping方向，将那些距离或方向存在异常的mapping位置丢弃，以此来排除最不可能的mappning情况

当参考基因组中存在多个相同或相似的拷贝区域时，会对CNV calling带来一定的误导

例如，在参考基因组(单倍体)中存在两个完全相同的片段重复，分别记为A和B，而在检测的样本中

只存在A拷贝，由于人是二倍体，样本会存在两份A拷贝，那么来着样本的A区域的reads在mapping过程中会随机均匀地分布在参考基因组的A区域和B区域(分值相同的multiple mapping随机指定一个mapping位置)，则A区域和B区域的Read Depth为实际的一半，因此这两个基因组区域都会被鉴定为Deletion

为了解决这种情况，CNVnator提出了一种解决方案：

对于那些multiple mapping的reads，会被赋予一个mapping quality，为0，这些reads被称为q0 reads

统计每个called CNV segments的q0 reads的比例，然后统计这些比例的频数分布，得到下面的直方图：

若一个CNV片段的q0比例大于50%，则暗示着它有可能是在参考基因组中存在高度相似的多个基因组区域，此时的CNV calling的结果可能不准确，需要进行过滤

假设我们有一堆点，和一个小的圆形窗口，我们要完成的任务就是将这个窗口移动到最大灰度密度处(也就是点最多的地方)。如下图所示：

初始窗口是蓝色的C1，它的圆心为蓝色方框的C1_o，而窗口中所有点质心却是C1_r，很明显圆心和点的质心没有重合。所以移动圆心C1_o到质心C1_r，这样我们就得到了一个新的窗口。这时又可以找到新的窗口内所以点的质心，大多数情况下还是不重合，所以重复上面的操作直到：将新窗口的圆心和它所包含的点的质心重合，这样我们的窗口会落在像素值(和)最大的地方。如上图C2是窗口的最后位置，它包含的像素点最多。

每一次迭代中需要移动圆心到质心，这个移动的方向称为Mean Shift向量

那这个Mean Shift向量怎么计算呢？

对于给定d维空间 Rd 中的n个样本点xi,i=1,2,…,n在xd点的Mean Shift向量的基本形式定义为：

其中，Sh是一个半径为h的高维球区域：Sh(x)={y:(y−x)T(y−x)≤h2}，k表示在这n个样本点中有k个落入球Sh中。

直观上来看，这k个样本点在x处的偏移向量即为：对落入Sh区域中的k个样本点相对于点x的偏移向量求和然后取平均值；

几何解释为：如果样本点xi服从一个概率密度函数为f(x)的分布，由于非零的概率密度函数的梯度指向概率密度增加最大的方向，因此从平均上来说，Sh区域内的样本点更多的落在沿着概率密度梯度的方向。因此，Mean Shift向量Mh(x)应该指向概率密度梯度的方向

参考资料：

(4) Zhou B, Ho SS, Zhang X, et al.Whole-genome sequencing analysis of CNV using low-coverage and paired-end strategies is efficient and outperforms array-based CNV analysis Journal of Medical Genetics 2018;55:735-743.

(5) Xie C , Tammi M T . CNV-seq, a new method to detect copy number variation using high-throughput sequencing[J]. BMC Bioinformatics, 2009, 10(1):80-0.

(6) Zhao M, Wang Q, Wang Q, Jia P, Zhao Z. Computational tools for copy number variation (CNV) detection using next-generation sequencing data: features and perspectives. BMC Bioinformatics. 2013;14 Suppl 11(Suppl 11):S1.

(13) Abyzov A, Urban AE, Snyder M, Gerstein M. CNVnator: an approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing. Genome Res. 2011;21(6):974-84.