CNVcaller遇到的坑

使用CNVcaller进行CNV检测

• 1.下载CNVcaller

cd /public/jychu/soft
git clone git://github.com/JiangYuLab/CNVcaller.git
如果下载不了，可以下载压缩文件再上传至服务器
unzip CNVcaller-master.zip         #解压
mv CNVcaller-master CNVcaller             #重命名
在运行前，需要基于 CNVcaller 的安装环境修改两个 shell 程序中 CNVcaller 的安装路径，默认是当前工作目录

• 2.构建参考基因组数据库
  作用：将参考基因组按用户指定大小（-w）的滑动窗口及一定的步长（内置是滑动窗口 大小的一半）分别统计基因组上每个窗口的 GC、repeat 及 gap 含量。

cp /public/jychu/refs/Gallus_gallus.GRCg6a.dna.toplevel.fa  /public/jychu/chicken_body_size/CNV/
perl  ~/soft/CNVcaller/bin/CNVReferenceDB.pl  Gallus_gallus.GRCg6a.dna.toplevel.fa -w 1000     #滑动窗口1000bp步长500bp
输出文件：上述运行完成会在工作目录下生成一个名为“referenceDB.1000”的文件，其 中的“1000”为滑动窗口的大小，文件的每列信息如下： 
第一列，染色体名称 
第二列，窗口在对应染色体上的序号 
第三列，窗口实际起始位置 
第四列，GC 含量 
第五列，重复序列含量 
第六列，gap 比例

• 3.计算每个窗口的绝对拷贝数
  作用：计算每个个体基因组所有窗口的绝对拷贝数。该过程主要包括三步：1. 解析每个 样本的 BAM（BWA 比对生成）文件，并统计每个窗口的读段数；2. 对高相似度的（≥97%） 窗口的读段数进行合并；3. 对每个窗口合并后的读段数进行 GC 偏好性校正并将其除以所 有窗口校正后读段数的中位数，以获得该窗口的绝对拷贝数。在这一步，每个样本运行需要 约 500 MB 内存，所以一个节点可以同时提交多个任务。
  输入：上一步输出的数据库文件“ReferenceDB.400”和 BAM 文件。默认数据库文件放 在当前目录下.输出：运行结束后，三个默认文件夹（RD_raw、RD_absolute、RD_normalized）将会在当前 目录下被创建，分别包含了每个样本的全基因组所有窗口原始读段数、通过 link 文件合并后 的读段数，以及每个样本经 GC 测序偏斜校正和标准化后的绝对拷贝数。标准化的文件名显 示了该个体基因组所有窗口的读段数平均值、标准差与性别（1 为 XX 或 ZZ，2 为 XY 或 ZW），其中平均值和标准差可用于样本的质控。绝对拷贝数为 1 时，表示为正常的拷贝数， 即正常二倍体；0.5 表示杂合缺失；0 表示纯合缺失；1.5 表示杂合重复；2 表示纯合重复； 绝对拷贝数超过 2 表示复杂的多次重复。

命令：bash Individual.Process.sh -b <bam> -h <header> -d <dup> -s <sex_chromosome>
鸡的dup文件下载： wget http://animal.nwsuaf.edu.cn/code/source/download/CNVcaller/database/duplicated_window_record_file/chicken.tar.gz（提供的为5.0版本）
tar -xzvf chicken.tar.gz  #解压文件
cp  ~/soft/CNVcaller/Individual.Process.sh ./
samtools view -H BC01.sox5.bam  #查看文件的标签
gatk AddOrReplaceReadGroups -I BC01.sox5.bam -O BC01.bam -LB WGS -PL illumina -PU BC01 -SM BC01 (如果未加标签)
ls *.bam|sed "s/\./\t/g"|cut -f1>sample.tmp
cat sample.tmp|while read id
do
arr=(${id})
bam=${arr[0]}
bash Individual.Process.sh -b /public/jychu/chicken_body_size/CNV/${bam}.sox5.bam -h $bam -d chickendup/Gallus5_1000_link -s none   &（-s Z）
done

如果物种没有提供的dup文件，可以按下面的步骤制作

Step 1: Split genome into short kmer sequences.
python ~/soft/CNVcaller/bin/0.1.Kmer_Generate.py Gallus_gallus.GRCg6a.dna.toplevel.fa 1000 chicken6_1000.fa
Step 2: Align the kmer FASTA (from step 1) to reference genome using blasr.
conda create -n blasr
conda activate blasr
conda install blasr
或者编译安装
git clone https://github.com/mchaisso/blasr.git  （或者下载压缩文件） 
blasr chicken6_1000.fa Gallus_gallus.GRCg6a.dna.toplevel.fa --sa Gallus_gallus.GRCg6a.dna.toplevel.fa.sa --out chicken6_1000.aln -m 5 --noSplitSubreads --minMatch 15 --maxMatch 20 --advanceHalf --advanceExactMatches 10 --fastMaxInterval --fastSDP --aggressiveIntervalCut --bestn 10
Gallus_gallus.GRCg6a.dna.toplevel.fa.sa 文件需要用balsr中的sawriter去创建
#sawriter的下载和使用（需要用root权限） 这是我弄了好几天才能编译成功的一种方式
git clone git://github.com/samtools/htslib.git    #安装blasr前需要安装htlib
cd htslib
https://github.com/samtools/htscodecs.git
autoheader     # If using configure, generate the header template...
autoconf       # ...and configure script (or use autoreconf to do both)
./configure    # Optional, needed for choosing optional functionality
make
make install
git clone git://github.com/mchaisso/blasr.git
export HDF5INCLUDEDIR=/usr/include/hdf
export HDF5LIBDIR=/usr/lib/hdf
cd blasr/
make
cd alignment/bin
/home/cjy/soft/blasr/alignment/bin/sawritermc Gallus_gallus.GRCg6a.dna.toplevel.fa   #利用sawriter对基因组进行index                         
我有直接编译好的sawriter可执行文件，如果你们安装不了，可以直接私信我要
Step 3: Generate duplicated window record file.
python ~/soft/CNVcaller/bin/0.2.Kmer_Link.py chicken6_1000.aln 1000 chicken1000.link

• 4.拷贝数变异区域的确定
  作用：通过综合考虑绝对拷贝数的分布、变异的频率及相邻窗口的显著相关性来初步确定 CNVR 的边界（primaryCNVR）。最后，将相邻且拷贝数在群体中分布显著相关的 CNVR
  进一步合并得到最终的拷贝数变异检测结果（mergedCNVR）。

cd /public/jychu/chicken_body_size/CNV/RD_normalized
ls `pwd`/*sex_1 >list
touch exclude_list  #新建一个空文件
cd ..
cp  ~/soft/CNVcaller/CNV.Discovery.sh ./
CNVcaller=~/soft/CNVcaller
bash CNV.Discovery.sh -l RD_normalized/list -e RD_normalized/exclude_list -f 0.1 -h 3 -r 0.5 -p primaryCNVR -m mergeCNVR
参数详解：
Required arguments: 
-l|--RDFileList          individual normalized read depth file list
-e|--excludedFileList    list of samples exclude from CNVR detection
-f|--frequency           minimum frequency of gain/loss individuals for candidate CNV window definition
                         [recommend 0.1]
-h|--homozygous          number of homozygous gain/loss individuals for candidate CNV window definition 
                         [recommend 3]
-r|--pearsonCorrelation  minimum of Pearson’s correlation coefficient between the two adjacent non-overlapping windows
                         0.5 for sample size (0, 30]
                         0.4 for sample size (30, 50]
                         0.3 for sample size (50, 100]
                         0.2 for sample size (100, 200]
                         0.15 for sample size (200, 500]
                         0.1 for sample size (500,+∞)
-p|--primaryCNVR         primary CNVR result 
-m|--mergedCNVR          merged CNVR result

• 5.基因型判定
  作用：利用混合高斯模型将每个样本的拷贝数归入不同的基因型分类，并以 VCF 格式输出， 以方便后续通过全基因组关联分析挖掘与重要经济性状有关的拷贝数变异。

cp  ~/soft/CNVcaller/Genotype.py ./
python Genotype.py --cnvfile mergeCNVR --outprefix Genotype --nproc 24
报错：ModuleNotFoundError: No module named '_sysconfigdata_x86_64_conda_linux_gnu'
find / -name '_sysconfigdata_x86_64*
/publicchu/miniconda3/pkgs/python-3.7.3-h0371630_0b/python3.7/_sysconfigdata_x86_64_conda_cos6_linux_gnu.py ~icken_body_size/CNV/
mv _sysconfigdata_x86_64_conda_cos6_linux_gnu.py _sysconfigdata_x86_64_conda_linux_gnu.py
报错：from sklearn.metrics import calinski_harabaz_score
ImportError: cannot import name 'calinski_harabaz_score' from 'sklearn.metrics'
修改calinski_harabaz_score 为calinski_harabasz_score
sed -i "s/harabaz/harabasz/g" Genotype.py

输出：genotypeCNVR.vcf 和genotypeCNVR.tsv。二者均包含每个个体的基因分型结果，以 及每个 CNVR 基因分型的似然值和轮廓系数，区别在于 genotypeCNVR.vcf 为 VCF 格式， genotypeCNVR.tsv 为 tab 分割的表格格式。此外，genotypeCNVR.tsv 还包含更多的字段用于 存放分类结果的统计信息，如每一类的频数及对应的绝对拷贝数的平均值与标准差。 VCF 各字段含义如下： CHROM: CNVR 所在的序列名称。 POS：CNVR 的起始坐标位置。 ID：CNVR 的编号，格式为：序列标号:起始位置-终止位置。 ALT：变异类型，包括 CN0、CN1、CN2 和 CNH，分别代表零个拷贝、一个拷贝、两个拷 贝和超过两个拷贝。 INFO：包含 CNVR 的终止位置（END），CNVR 的变异类型（SVTYPE），基因分型的对数 似然值（LOGLIKELIHOOD）和轮廓系数（SILHOUETTESCORE）。 FORMAT：每个个体基因分型结果的输出格式，GT 和 CP 分别代表个体的基因分型结果和 绝对拷贝数。
参数详解： --cnvfile CNVR 结果文件，包含全部样本的拷贝数信息，由上一步的 CNV.Discovery.sh 得到。 --outprefix 输出结果文件的前缀，默认会输出两个文件，其中，后缀为 tsv 的文件记录了分 类结果的基本统计信息，方便后续过滤低质量的 CNVR；后缀为 vcf 的文件为常规 VCF 格 式，各字段具体含义见 VCF 注释部分。 可选参数：--merge为了得到更多的双等位 CNVR 用于后续分析，可以使用--merge 选项。使用后，会增加一个后缀为 _merge.vcf 的 输 出 文 件 ， 类 似 .vcf 文 件 ， 区 别 在 于 _merge.vcf 中把所有重复算作一种变异类型