群体遗传变异——拷贝数变异分析
基因组发生重排而导致的,一般指长度1 kb 以上的基因组片段的拷贝数增加或者减少, 主要表现为亚显微水平的重复或者缺失。因此称为“微”缺失/重复变异。
在基因组水平上，CNV 不像单核苷酸多态性（SNPs）那样高频发生，但是它们对生物的表型、功能、适应性等有重要作用，且其产生的影响比 SNP 要更加强烈。CNVs 是驱动基因和基因组进化的重要机制。越来越多的研究表明，CNVs 通过各种分子机制（如基因剂量、基因融合）对家畜的孟德尔性状和数量性状的表型变异具有重要影响。如果CNV 存在于蛋白编码区，则改变蛋白功能，如果在调控区，则改变基因表达水平。
基因组上非等位的两个高度同源的DNA序列在减数分裂或者有丝分裂的过程中发生错误的配对，并发生序列交换，从而导致缺失、重复、倒位的出现。——就产生了拷贝数变异。
原文请看这儿https://zhuanlan.zhihu.com/p/166445484
CNV为基因组序列中长度为 50 bp - 5 Mb 的重复和缺失变异，在群体中，不同个体间重叠区域 CNV 的整合结果称作 CNV 区域（CNV Region, CNVR）。当这些 CNVs 被发现在蛋白质编码基因（CNV 基因）或 miRNA（CNV-miRNA）的基因组区域时，该遗传变异就可能作用于分子调控机制，并影响着生物的表型多样性和对疾病的易感性。
非等位同源重组（Non-allelic homologous recombination, NAHR）、非同源末端连接（Non-homologous end joining, NHEJ）、复制叉停滞与模板交换（Fork stalling and template switching, FoSTeS）、L1 介导的反转录转座（L1-mediated retrotranspositio, LINE1）是基因组中形成重排的主要机制，同时也是大部分 CNV 产生的原因。
由于非同源染色体之间的两个相似序列区域容易发生重组，所以 NAHR常发生在减数分裂和有丝分裂过程。姐妹染色单体之间发生交叉的过程可能增加 DNA片段或丢失另一个 DNA 片段，从而导致染色体片段的重复、缺失和倒位。
NHEJ 机制是细胞用来修复由电离辐射或活性氧物质引起的 DNA 双链断裂（Double-strand breaks, DSBs）的生理形式。
FoSTeS 是一种基于 DNA 复制的机制，可以解释复杂的基因组重排和 CNV。
L1 转座通过逆转录和整合方式引发 CNV 的产生。基因组中倒位和易位等组合事件的鉴别在技术上仍然具有挑战性，要完全鉴别其变异类型和成因的成本也很高，相对而言，CNV 更容易被鉴别，且可以为分析复杂变异提供强有力的遗传学证据。
image.png
检测手段：
芯片法和测序法
芯片法主要包括比较基因组杂交芯片（Array comparative genomic hybridization, a CGH）和 SNP 芯片（Singlenucleotide polymorphism arrays）。
DNA 测序法主要包括全基因组测序（Whole genome sequnecing, WGS）和单分子测序。
Paired-end mapping（PEM）方法、Split read（SR）方法、Read depth（RD）方法和 De novo 组装是检测基因组拷贝数变异的主要策略。
RD 方法通过利用短读取序列数据与参考基因组进行比对，标准化每个区域的读取深度，其中低或零深度的区域被解释为缺失，深度增加的区域被解释为拷贝数的重复或扩增。RD 方法相较于 PEM 和 SR 方法有着较高的检出率和准确性。
image.png
>现在基于二代测序数据分析 CNV的研究还存在着样本少、测序深度较低、选用参考基因组质量差异大等不足之处，这些因素均会对 CNV 检测结果产生不利影响
本文介绍的方法分别有：CNVCaller、Matchclips2、Delly以及CNVnator
1、CNVCaller
CNVCaller软件是由西北农林科技大学姜雨老师团队开发的一款专门用于分析拷贝数变异的软件（CNVcaller: highly efficient and widely applicable software for detecting copy number variations in large populations），该软件的假阳性率相较于其他软件更低、运算速度更高，是一款本人强烈推荐的软件。
详情可以参考姜雨老师团队主导开发的数据库omicsDB (nwsuaf.edu.cn)
安装参考：JiangYuLab/CNVcaller (github.com)
使用前准备：
在Individual.Process.sh和CNV.Discovery.sh脚本文件中修改CNVcaller 的安装路径（该软件的安装路径），默认是当前工作目录
image.png
即：export CNVcaller=/home/sll/miniconda3/CNVcaller‘
注意：该软件在使用时每一步的输入以及输出文件都用绝对路径
1）在当前目录创建referenceDB.windowsize文件，构建参考基因组数据库
perl /home/sll/miniconda3/CNVcaller/bin/CNVReferenceDB.pl /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna -w 800
-w the window size (bp) for all samples [default=800] 滑动窗口大小， >10x使用400-1000，<10x使用1000-2000
-l the lower limit of GC content [default=0.2] 窗口可能含有的最低GC比例，低于该比例的窗口不会进入后续计算。
-u the upper limit of GC content [default=0.7] 窗口可能含有的最高GC比例
-g the upper limit of gap content [default=0.5] 窗口可能含有的最高gap比例
作用：将参考基因组按用户指定大小（-w）的滑动窗口及一定的步长（内置是滑动窗口大小的一半）分别统计基因组上每个窗口的 GC、repeat 及 gap 含量。
生成一个“referenceDB.800”的文件。
第一列，染色体名称
第二列，窗口在对应染色体上的序号
第三列，窗口实际起始位置
第四列，GC 含量
第五列，重复序列含量
第六列，gap 比例
2）计算每个窗口的绝对拷贝数
在这里有个问题，我们只能使用他们提供的800bp窗口的dup文件来做，而他们提供的生成这个文件的脚本中，blasr软件中的sawriter安装不了，问了软件开发作者，也说这个软件的安装困难。因此，目前对我来说自己生成不了。
bash Individual.Process.sh -b /home/sll/2022-DNA-cattle/CNV/YLHN1.sorted.addhead.markdup.bam -h YLHN1 -d Btau5.0.1_800_link -s none
-h BAM 文件的标头信息，需要与 BAM 文件中 SM 标签保持一致（用户可以通过 samtools view -H BAM 查看）。
-d 校正所需要的dup文件。查看物种dup文件[JiangYuLab/CNVcaller (github.com)](https://github.com/JiangYuLab/CNVcaller/tree/master#generate-your-own-duplicated-window-record-file)。
-s 性染色体名称。CNVcaller 会根据给定的性染色体所有窗口读段数的中位数与所有常染色体窗口读段数的中位数比例来确定该个体的性别。
运行结束后，三个默认文件夹（RD_raw、RD_absolute、RD_normalized）将会在当前目录下被创建，分别包含了每个样本的全基因组所有窗口原始读段数、通过 link 文件合并后的读段数，以及每个样本经 GC 测序偏斜校正和标准化后的绝对拷贝数。标准化的文件名显示了该个体基因组所有窗口的读段数平均值、标准差与性别（1 为 XX 或 ZZ，2 为 XY 或ZW），其中平均值和标准差可用于样本的质控。绝对拷贝数为 1 时，表示为正常的拷贝数，即正常二倍体；0.5 表示杂合缺失；0 表示纯合缺失；1.5 表示杂合重复；2 表示纯合重复；绝对拷贝数超过 2 表示复杂的多次重复。
3)拷贝数变异区域的确定（多样本合并，单个样本会报错）
通过综合考虑绝对拷贝数的分布、变异的频率及相邻窗口的显著相关性来初步确定 CNVR 的边界（primaryCNVR）。最后，将相邻且拷贝数在群体中分布显著相关的 CNVR进一步合并得到最终的拷贝数变异检测结果（mergedCNVR）。
cd RD_normalized
记得把第一步生成的referenceDB.800放进去,输入文件用绝对路径！！！
list 上一步输出的RD_normalized文件夹中的文件，一行一个，多个样本，需包含绝对路径。
touch exclude_list
exclude_list 为空文件，提前创建
bash /home/sll/miniconda3/CNVcaller/CNV.Discovery.sh -l /home/sll/2022-DNA-cattle/CNV/RD_normalized/list.txt -e /home/sll/2022-DNA-cattle/CNV/RD_normalized/exclude_list -f 0.1 -h 1 -r 0.1 -p primaryCNVR -m mergeCNVR
-l 个体经绝对拷贝数校正后的结果文件列表
-e 该列表中的样本不被用于 CNVR 的检测，但结果文件会记录这些样本的绝对拷贝数。exclude_list 为空文件代表所有个体将被用于CNVR检测。
-f 当一个窗口有超过该频率的个体绝对拷贝数与正常拷贝（“1”）显著差异（杂合删除或者
杂合复制）时，就定义该窗口为候选拷贝数变异窗口。
-h 当一个窗口有超过该频数的个体绝对拷贝数与正常拷贝（“1”）显著差异（纯合删除或者
纯合复制）时，就定义该窗口为候选拷贝数变异窗口。
##只需要满足-f 与-h 中的任意一个即可
-r 在定义 CNVR 时，如果相邻（没有 overlap）候选拷贝数变异窗口的绝对拷贝数的相关系数高于该值将被合并。
4)基因型判定
利用混合高斯模型将每个样本的拷贝数归入不同的基因型分类，并以 VCF 格式输出，以方便后续通过全基因组关联分析挖掘与重要经济性状有关的拷贝数变异。
python /home/sll/miniconda3/CNVcaller/Genotype.py --cnvfile mergeCNVR --outprefix genotypeCNVR --nproc 24
--cnvfile CNVR 结果文件，包含全部样本的拷贝数信息，上一步的输出结 果（mergeCNVR）。
--outprefix 输出结果文件的前缀，默认会输出两个文件，其中，后缀为 tsv 的文件记录了分类结果的基本统计信息，方便后续过滤低质量的 CNVR；后缀为 vcf 的文件为常规 VCF 格式。
--nproc 程序使用的进程数，默认为单进程，使用此参数可显著减少程序运行时间，但会增加内存消耗。
可能会报错
image.png
修改Genotype.py脚本中的harabaz成harabasz即可。
sed -i "s/harabaz/harabasz/g" /home/sll/miniconda3/CNVcaller/Genotype.py
结果文件
genotypeCNVR.vcf文件内容
CHROM: 所在染色体
POS：CNVR的起始坐标位置
ID：CNVR的编号，格式为：序列坐标:起始位置-终止位置
ALT: 变异类型，包括CN0、CN1、CN2和CNH，分别代表0个拷贝、1个拷贝、2个拷贝和超过2个的拷贝
INFO：包含CNVR的终止位置（END），CNVR的变异类型（SVTYPE），基因分型的对数自然值（LOGLIKELIHOOD）和轮廓系数（SILHOUETTESCORE）
FORMAT：每个个体基因分型结果的输出格式，GT和GP分别代表个体的基因分型结果和绝对拷贝数。
image.png
shell循环(敬请膜拜大佬，也就是我，哈哈哈)：
touch CNVCaller.sh
# 对基因组创建窗口文件(以后直接复制到操作目录下用即可)
perl /home/sll/miniconda3/CNVcaller/bin/CNVReferenceDB.pl /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna -w 800
# 计算每个窗口的绝对拷贝数
ls *bam|cut -d"." -f 1 | sort -u | while read id; do bash /home/sll/miniconda3/CNVcaller/Individual.Process.sh -b `pwd`/${id}.sorted.addhead.markdup.bam -h ${id} -d /home/sll/miniconda3/CNVcaller/Btau5.0.1_800_link -s none; done
# 将referenceDB.800文件复制到RD_normalized目录下
cp referenceDB.800 RD_normalized
# 进入RD_normalized目录
cd RD_normalized
# 将RD_normalized目录下新生成的sex_1结尾的文件名，以绝对路径的形式写入list.txt中
ls -R `pwd`/*sex_1 > list.txt
# 新建exclude_list文件
touch exclude_list
# 拷贝数变异区域的确定
bash /home/sll/miniconda3/CNVcaller/CNV.Discovery.sh -l `pwd`/list.txt -e `pwd`/exclude_list -f 0.1 -h 1 -r 0.1 -p primaryCNVR -m mergeCNVR
# 基因型判定
python /home/sll/miniconda3/CNVcaller/Genotype.py --cnvfile mergeCNVR --outprefix genotypeCNVR --nproc 8
nohup bash CNVCaller.sh &
2、Matchclips2（刘正喜师姐的脚本，还得是师姐）
基于long soft clips的CNV断点计算方法
1)每个样品分别进行matchclips2计算：
/home/software/matchclips2/matchclips -t 4 -f /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna -b BC-448.sorted.addhead.markdup.bam -o CNV.BC-448.txt
结果文件：
1、CHROM
2、起始位置
3、中止位置
4、DEL或DUP
For DEL, the 5' break point is BEGIN and 3' is END,
For DUP, the 5' break point is END and 3' is BEGIN.、
5、区域长度：3' break point - 5' break point
6、UN，碱基数量This number tells how many bases are repeated at the two break points
7、RD:n1;n2;n3:s
2)合并结果：
/home/software/matchclips2/cnvtable -L 10000 -cnvf cnvf.txt -O 0.5 -o overlap.txt
cnvf.txt为含有上一步结果文件名称的txt文件，一列一个
-cnv STR file names separated by white spaces
-i STR IDs for the files separated by white spaces
-cnvf STR a file with list of file names and(or) ids
-l INT minimum length to include, INT=3
-L INT maximum length to include, INT=10000000
-t STR only process STR type of CNVs
-R STR subset cnvs in region STR
-chr 1-22XY are treated as chr[1-22XY]
-O FLOAT minimum reciprocal overlap ratio [0.0, 1.0], FLOAT=0.5
-o STR outputfile, STR=STDOUT
shell循环
touch matchclips2.sh
##运行
ls *bam|cut -d"." -f 1 | sort -u | while read id; do /home/software/matchclips2/matchclips -t 4 -f /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna -b ${id}.sorted.addhead.markdup.bam -o ${id}.cnv; done
#将生成文件名列到一个文件中
find -name '*cnv' -exec basename {} \; > cnvf.txt
#进行合并
/home/software/matchclips2/cnvtable -L 10000 -cnvf cnvf.txt -O 0.5 -o overlap.txt
nohup bash matchclips2.sh &
3、Delly
Delly v0.7.8 可以一次性call 所有类型的SVs；低版本的需要通过-t指定SV类型
输入文件
1、bam文件：Bam files need to be sorted, indexed and ideally duplicate marked. If multiple libraries are present for a single sample these need to be merged in a single bam file with unique ReadGroup tags.
2、.excl: a file include information of telomere and centromere regions and also all unplaced contigs, those regions will be removed from calling.
有群体时，建议将所有个体的.bcf文件合并后，再进行genotyping。
运行：详情参考https://github.com/yhwu/matchclips2
1)call sv
delly call -g /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna BC-448.sorted.addhead.markdup.bam > BC-448.cnv.vcf
2）bcf转vcf格式
bcftools view $i.target.bcf >$i.target.vcf
3)Merge SV sites into a unified site list(合并)
delly merge -o sites.bcf s1.bcf s2.bcf ... sN.bcf
4)Genotype this merged SV site list across all samples. This can be run in parallel for each sample(基因分型，每个样本都做)
delly call -g $DB/target-ref.fa -v sites.bcf -o s1.geno.bcf s1.bam
5）Merge all genotyped samples to get a single VCF/BCF using bcftools merge
bcftools merge -m id -O b -o merged.bcf s1.geno.bcf s2.geno.bcf ... sN.geno.bcf
插入（Insertion, INS）
缺失（Deletion, DEL）
反转（Inversion, INV）
染色体内部易位（Intra-chromosomal Translocation, ITX）
染色体间易位（Inter-chromosomal Translocation, CTX)
4、CNVnator(拉跨）
检测出的片段长度很长，假阳性结果率较低，还算比较准确的一个软件
1)提取mapping信息
home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -tree BC-448.sorted.addhead.markdup.bam -unique
提取比对后的reads，构建树，如果程序之前运行失败，已经生成了一个root文件，在下次重新运行时一定要删除该root文件，如果不删除，新的分析结果会追加到错误的root文件中，影响后续分析。另外，可以添加 -chrom 函数指定染色体 如 -chrom 1 2 代表选择1，2号染色体。
2)生成质量分布图HISTOGRAM
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -d -his 100
# 100指的是bin size，可以通过-eval参数进行筛选，也可以根据经验值进行确定，一般测序深度20-30x选取bin size大小100，2-3x选取500，100x选取30。
# -d指定目录，内部存放给染色体的fasta文件，该参数指针对报错信息显示不能解析基因组文件，此时需要指定参考基因组的位置，且各染色体需要拆分成单独的fasta文件（目录下有其文件也可以，只要有所有染色体的序列就好，软件会自动识别）
3)生成统计结果
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -stat 100
4)RD（reads depth）信息分割partipition
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -partition 100
##cnvnator在进行变异检测时，以提供的bin size对整个基因组进行切割，之后按照RD(read-depth)为基准进行cnv的检测。
注意：从第二步开始设置的bin size就不能变了
5)变异检出
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -call 100 > BC-448-cnvout.txt
运行之后输出结果如下：
#第一列变异类型
#第二列位点信息
#第三列CNV大小
#第四列为标准化参数normalized_RD
#第5-8列为e-value值，其中第五列越小，说明结果越准确
#第九列q0质量值
6)转为vcf格式
/home/software/CNVnator_v0.4.1/src/cnvnator2VCF.pl cnv.call.txt > cnv.vcf
选择清除与环境适应性位点挖掘--Vst分析
Vst分析是类似于Fst的一个指标，用来衡量群体间每个CNVR差异大小的统计量，计算方法为Vst＝（Vt－Vs）/Vt，Vt表示所有样本该区域拷贝数大小的方差，Vs表示两个群体各自的方差根据各自群体大小加权之后的值。Vst的值介于0-1之间，值越大表示群体间该区域拷贝数变异差异越大，反之则越小。
EXCEL计算.png
