基因組發(fā)生重排而導(dǎo)致的,一般指長(zhǎng)度1 kb 以上的基因組片段的拷貝數(shù)增加或者減少, 主要表現(xiàn)為亞顯微水平的重復(fù)或者缺失具被。因此稱為“微”缺失/重復(fù)變異发乔。
在基因組水平上掸犬，CNV 不像單核苷酸多態(tài)性（SNPs）那樣高頻發(fā)生菜秦，但是它們對(duì)生物的表型廷区、功能、適應(yīng)性等有重要作用丰泊，且其產(chǎn)生的影響比 SNP 要更加強(qiáng)烈薯定。CNVs 是驅(qū)動(dòng)基因和基因組進(jìn)化的重要機(jī)制始绍。越來(lái)越多的研究表明瞳购，CNVs 通過(guò)各種分子機(jī)制（如基因劑量、基因融合）對(duì)家畜的孟德?tīng)栃誀詈蛿?shù)量性狀的表型變異具有重要影響亏推。如果CNV 存在于蛋白編碼區(qū)学赛，則改變蛋白功能，如果在調(diào)控區(qū)吞杭，則改變基因表達(dá)水平盏浇。
基因組上非等位的兩個(gè)高度同源的DNA序列在減數(shù)分裂或者有絲分裂的過(guò)程中發(fā)生錯(cuò)誤的配對(duì)，并發(fā)生序列交換芽狗，從而導(dǎo)致缺失绢掰、重復(fù)、倒位的出現(xiàn)〉尉ⅲ——就產(chǎn)生了拷貝數(shù)變異攻晒。
原文請(qǐng)看這兒https://zhuanlan.zhihu.com/p/166445484
CNV為基因組序列中長(zhǎng)度為 50 bp - 5 Mb 的重復(fù)和缺失變異，在群體中班挖，不同個(gè)體間重疊區(qū)域 CNV 的整合結(jié)果稱作 CNV 區(qū)域（CNV Region, CNVR）鲁捏。當(dāng)這些 CNVs 被發(fā)現(xiàn)在蛋白質(zhì)編碼基因（CNV 基因）或 miRNA（CNV-miRNA）的基因組區(qū)域時(shí)，該遺傳變異就可能作用于分子調(diào)控機(jī)制萧芙，并影響著生物的表型多樣性和對(duì)疾病的易感性给梅。
非等位同源重組（Non-allelic homologous recombination, NAHR）、非同源末端連接（Non-homologous end joining, NHEJ）双揪、復(fù)制叉停滯與模板交換（Fork stalling and template switching, FoSTeS）动羽、L1 介導(dǎo)的反轉(zhuǎn)錄轉(zhuǎn)座（L1-mediated retrotranspositio, LINE1）是基因組中形成重排的主要機(jī)制，同時(shí)也是大部分 CNV 產(chǎn)生的原因盟榴。
由于非同源染色體之間的兩個(gè)相似序列區(qū)域容易發(fā)生重組曹质，所以 NAHR常發(fā)生在減數(shù)分裂和有絲分裂過(guò)程。姐妹染色單體之間發(fā)生交叉的過(guò)程可能增加 DNA片段或丟失另一個(gè) DNA 片段擎场，從而導(dǎo)致染色體片段的重復(fù)羽德、缺失和倒位。
NHEJ 機(jī)制是細(xì)胞用來(lái)修復(fù)由電離輻射或活性氧物質(zhì)引起的 DNA 雙鏈斷裂（Double-strand breaks, DSBs）的生理形式迅办。
FoSTeS 是一種基于 DNA 復(fù)制的機(jī)制宅静，可以解釋復(fù)雜的基因組重排和 CNV。
L1 轉(zhuǎn)座通過(guò)逆轉(zhuǎn)錄和整合方式引發(fā) CNV 的產(chǎn)生站欺∫碳校基因組中倒位和易位等組合事件的鑒別在技術(shù)上仍然具有挑戰(zhàn)性，要完全鑒別其變異類型和成因的成本也很高矾策，相對(duì)而言磷账，CNV 更容易被鑒別，且可以為分析復(fù)雜變異提供強(qiáng)有力的遺傳學(xué)證據(jù)贾虽。
檢測(cè)手段：
芯片法和測(cè)序法
芯片法主要包括比較基因組雜交芯片（Array comparative genomic hybridization, a CGH）和 SNP 芯片（Singlenucleotide polymorphism arrays）逃糟。
DNA 測(cè)序法主要包括全基因組測(cè)序（Whole genome sequnecing, WGS）和單分子測(cè)序。
Paired-end mapping（PEM）方法蓬豁、Split read（SR）方法绰咽、Read depth（RD）方法和 De novo 組裝是檢測(cè)基因組拷貝數(shù)變異的主要策略。
RD 方法通過(guò)利用短讀取序列數(shù)據(jù)與參考基因組進(jìn)行比對(duì)地粪，標(biāo)準(zhǔn)化每個(gè)區(qū)域的讀取深度取募，其中低或零深度的區(qū)域被解釋為缺失，深度增加的區(qū)域被解釋為拷貝數(shù)的重復(fù)或擴(kuò)增蟆技。RD 方法相較于 PEM 和 SR 方法有著較高的檢出率和準(zhǔn)確性玩敏。
>現(xiàn)在基于二代測(cè)序數(shù)據(jù)分析 CNV的研究還存在著樣本少斗忌、測(cè)序深度較低、選用參考基因組質(zhì)量差異大等不足之處旺聚，這些因素均會(huì)對(duì) CNV 檢測(cè)結(jié)果產(chǎn)生不利影響
本文介紹的方法分別有：CNVCaller飞蹂、Matchclips2、Delly以及CNVnator翻屈、LUMPY
1陈哑、CNVCaller
CNVCaller軟件是由西北農(nóng)林科技大學(xué)姜雨老師團(tuán)隊(duì)開(kāi)發(fā)的一款專門用于分析拷貝數(shù)變異的軟件（CNVcaller: highly efficient and widely applicable software for detecting copy number variations in large populations），該軟件的假陽(yáng)性率相較于其他軟件更低伸眶、運(yùn)算速度更高惊窖，是一款本人強(qiáng)烈推薦的軟件。
詳情可以參考姜雨老師團(tuán)隊(duì)主導(dǎo)開(kāi)發(fā)的數(shù)據(jù)庫(kù)omicsDB (nwsuaf.edu.cn)
安裝參考：JiangYuLab/CNVcaller (github.com)
使用前準(zhǔn)備：
在Individual.Process.sh和CNV.Discovery.sh腳本文件中修改CNVcaller 的安裝路徑（該軟件的安裝路徑）厘贼，默認(rèn)是當(dāng)前工作目錄
image.png
即：export CNVcaller=/home/sll/miniconda3/CNVcaller‘
注意：該軟件在使用時(shí)每一步的輸入以及輸出文件都用絕對(duì)路徑
1）在當(dāng)前目錄創(chuàng)建referenceDB.windowsize文件界酒，構(gòu)建參考基因組數(shù)據(jù)庫(kù)
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] 滑動(dòng)窗口大小， >10x使用400-1000嘴秸，<10x使用1000-2000
-l the lower limit of GC content [default=0.2] 窗口可能含有的最低GC比例毁欣，低于該比例的窗口不會(huì)進(jìn)入后續(xù)計(jì)算。
-u the upper limit of GC content [default=0.7] 窗口可能含有的最高GC比例
-g the upper limit of gap content [default=0.5] 窗口可能含有的最高gap比例
作用：將參考基因組按用戶指定大性榔（-w）的滑動(dòng)窗口及一定的步長(zhǎng)（內(nèi)置是滑動(dòng)窗口大小的一半）分別統(tǒng)計(jì)基因組上每個(gè)窗口的 GC凭疮、repeat 及 gap 含量。
生成一個(gè)“referenceDB.800”的文件串述。
第一列执解，染色體名稱
第二列，窗口在對(duì)應(yīng)染色體上的序號(hào)
第三列纲酗，窗口實(shí)際起始位置
第四列衰腌，GC 含量
第五列，重復(fù)序列含量
第六列觅赊，gap 比例
2）計(jì)算每個(gè)窗口的絕對(duì)拷貝數(shù)
在這里有個(gè)問(wèn)題右蕊，我們只能使用他們提供的800bp窗口的dup文件來(lái)做，而他們提供的生成這個(gè)文件的腳本中吮螺，blasr軟件中的sawriter安裝不了饶囚，問(wèn)了軟件開(kāi)發(fā)作者，也說(shuō)這個(gè)軟件的安裝困難规脸。因此坯约，目前對(duì)我來(lái)說(shuō)自己生成不了(已解決熊咽，看后面的Dup文件的生成)莫鸭。
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 文件的標(biāo)頭信息，需要與 BAM 文件中 SM 標(biāo)簽保持一致（用戶可以通過(guò) 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 會(huì)根據(jù)給定的性染色體所有窗口讀段數(shù)的中位數(shù)與所有常染色體窗口讀段數(shù)的中位數(shù)比例來(lái)確定該個(gè)體的性別梨与。
運(yùn)行結(jié)束后堕花，三個(gè)默認(rèn)文件夾（RD_raw、RD_absolute粥鞋、RD_normalized）將會(huì)在當(dāng)前目錄下被創(chuàng)建缘挽，分別包含了每個(gè)樣本的全基因組所有窗口原始讀段數(shù)、通過(guò) link 文件合并后的讀段數(shù)呻粹，以及每個(gè)樣本經(jīng) GC 測(cè)序偏斜校正和標(biāo)準(zhǔn)化后的絕對(duì)拷貝數(shù)壕曼。標(biāo)準(zhǔn)化的文件名顯示了該個(gè)體基因組所有窗口的讀段數(shù)平均值、標(biāo)準(zhǔn)差與性別（1 為 XX 或 ZZ等浊，2 為 XY 或ZW）腮郊，其中平均值和標(biāo)準(zhǔn)差可用于樣本的質(zhì)控。絕對(duì)拷貝數(shù)為 1 時(shí)筹燕，表示為正常的拷貝數(shù)轧飞，即正常二倍體；0.5 表示雜合缺失撒踪；0 表示純合缺失过咬；1.5 表示雜合重復(fù)；2 表示純合重復(fù)制妄；絕對(duì)拷貝數(shù)超過(guò) 2 表示復(fù)雜的多次重復(fù)援奢。
3)拷貝數(shù)變異區(qū)域的確定（多樣本合并，單個(gè)樣本會(huì)報(bào)錯(cuò)）
通過(guò)綜合考慮絕對(duì)拷貝數(shù)的分布忍捡、變異的頻率及相鄰窗口的顯著相關(guān)性來(lái)初步確定 CNVR 的邊界（primaryCNVR）集漾。最后，將相鄰且拷貝數(shù)在群體中分布顯著相關(guān)的 CNVR進(jìn)一步合并得到最終的拷貝數(shù)變異檢測(cè)結(jié)果（mergedCNVR）砸脊。
cd RD_normalized
記得把第一步生成的referenceDB.800放進(jìn)去,輸入文件用絕對(duì)路徑＞咂！凌埂！
list 上一步輸出的RD_normalized文件夾中的文件驱显，一行一個(gè)，多個(gè)樣本瞳抓，需包含絕對(duì)路徑埃疫。
touch exclude_list
exclude_list 為空文件，提前創(chuàng)建
perl版本高于5.10.1時(shí)孩哑，修改
CNVcaller/bin/2.2.CNVRRedundancy.pl文件中的第48行
die "unexpected sample number!\n" unless length(@tmp_start_array) == length(@tmp);
為die "unexpected sample number!\n" unless length(scalar(@tmp_start_array)) == length(scalar(@tmp));
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 個(gè)體經(jīng)絕對(duì)拷貝數(shù)校正后的結(jié)果文件列表
-e 該列表中的樣本不被用于 CNVR 的檢測(cè)栓霜，但結(jié)果文件會(huì)記錄這些樣本的絕對(duì)拷貝數(shù)。exclude_list 為空文件代表所有個(gè)體將被用于CNVR檢測(cè)横蜒。
-f 當(dāng)一個(gè)窗口有超過(guò)該頻率的個(gè)體絕對(duì)拷貝數(shù)與正掣炻拷貝（“1”）顯著差異（雜合刪除或者
雜合復(fù)制）時(shí)销凑，就定義該窗口為候選拷貝數(shù)變異窗口。
-h 當(dāng)一個(gè)窗口有超過(guò)該頻數(shù)的個(gè)體絕對(duì)拷貝數(shù)與正辰龃叮拷貝（“1”）顯著差異（純合刪除或者
純合復(fù)制）時(shí)斗幼，就定義該窗口為候選拷貝數(shù)變異窗口。
##只需要滿足-f 與-h 中的任意一個(gè)即可
-r 在定義 CNVR 時(shí)抚垄，如果相鄰（沒(méi)有 overlap）候選拷貝數(shù)變異窗口的絕對(duì)拷貝數(shù)的相關(guān)系數(shù)高于該值將被合并蜕窿。
4)基因型判定
利用混合高斯模型將每個(gè)樣本的拷貝數(shù)歸入不同的基因型分類，并以 VCF 格式輸出呆馁，以方便后續(xù)通過(guò)全基因組關(guān)聯(lián)分析挖掘與重要經(jīng)濟(jì)性狀有關(guān)的拷貝數(shù)變異渠羞。
python /home/sll/miniconda3/CNVcaller/Genotype.py --cnvfile mergeCNVR --outprefix genotypeCNVR --nproc 24
--cnvfile CNVR 結(jié)果文件，包含全部樣本的拷貝數(shù)信息智哀，上一步的輸出結(jié) 果（mergeCNVR）次询。
--outprefix 輸出結(jié)果文件的前綴，默認(rèn)會(huì)輸出兩個(gè)文件瓷叫，其中屯吊，后綴為 tsv 的文件記錄了分類結(jié)果的基本統(tǒng)計(jì)信息，方便后續(xù)過(guò)濾低質(zhì)量的 CNVR摹菠；后綴為 vcf 的文件為常規(guī) VCF 格式盒卸。
--nproc 程序使用的進(jìn)程數(shù)，默認(rèn)為單進(jìn)程次氨，使用此參數(shù)可顯著減少程序運(yùn)行時(shí)間蔽介，但會(huì)增加內(nèi)存消耗。
可能會(huì)報(bào)錯(cuò)
image.png
修改Genotype.py腳本中的harabaz成harabasz即可煮寡。
sed -i "s/harabaz/harabasz/g" /home/sll/miniconda3/CNVcaller/Genotype.py
結(jié)果文件
genotypeCNVR.vcf文件內(nèi)容
CHROM: 所在染色體
POS：CNVR的起始坐標(biāo)位置
ID：CNVR的編號(hào)虹蓄，格式為：序列坐標(biāo):起始位置-終止位置
ALT: 變異類型，包括CN0幸撕、CN1薇组、CN2和CNH，分別代表0個(gè)拷貝坐儿、1個(gè)拷貝律胀、2個(gè)拷貝和超過(guò)2個(gè)的拷貝
INFO：包含CNVR的終止位置（END），CNVR的變異類型（SVTYPE）貌矿，基因分型的對(duì)數(shù)自然值（LOGLIKELIHOOD）和輪廓系數(shù)（SILHOUETTESCORE）
FORMAT：每個(gè)個(gè)體基因分型結(jié)果的輸出格式炭菌，GT和GP分別代表個(gè)體的基因分型結(jié)果和絕對(duì)拷貝數(shù)。
image.png
Dup文件的創(chuàng)建
blasr軟件安裝建議從github下載壓縮文件上傳
sawriter是目錄alignment/bin/下的sawritermc
創(chuàng)建自己所需要的dup文件:
- Split genome into short kmer sequences
python 0.1.Kmer_Generate.py [OPTIONS] FAFILE WINSIZE OUTFILE
<FAFILE> Reference sequence in FASTA format
<WINSIZE> The size of the window to use for CNV calling
<OUTFILE> Output kmer file in FASTA format
python /home/sll/miniconda3/CNVcaller/bin/0.1.Kmer_Generate.py /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna 400 kmer.fa
- Align the kmer FASTA (from step 1) to reference genome using blasr.
sawriter 創(chuàng)建sa索引文件
/home/sll/software/blasr-master/alignment/bin/sawritermc GCF_002263795.1_ARS-UCD1.2_genomic.fna.sa GCF_002263795.1_ARS-UCD1.2_genomic.fna
blasr 生成比對(duì)后的kmer.aln文件
blasr kmer.fa GCF_002263795.1_ARS-UCD1.2_genomic.fna --sa GCF_002263795.1_ARS-UCD1.2_genomic.fna.sa \
--out kmer.aln -m 5 \
--noSplitSubreads --minMatch 15 --maxMatch 20 \
--advanceHalf --advanceExactMatches 10 --fastMaxInterval \
--fastSDP --aggressiveIntervalCut --bestn 10
- Generate duplicated window record file.
python 0.2.Kmer_Link.py [OPTIONS] BLASR WINSIZE OUTFILE
<BLASR> blasr results (-m 5 format)
<WINSIZE> The size of the window to use for CNV calling
<OUTFILE> Output genome duplicated window record file
python /home/sll/miniconda3/CNVcaller/bin/0.2.Kmer_Link.py kmer.aln 1000 Bos_ARS1.2_window.link
shell循環(huán)(敬請(qǐng)膜拜大佬逛漫，也就是我黑低，哈哈哈)：
touch CNVCaller.sh
# 對(duì)基因組創(chuàng)建窗口文件(以后直接復(fù)制到操作目錄下用即可)
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
# 計(jì)算每個(gè)窗口的絕對(duì)拷貝數(shù)
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文件復(fù)制到RD_normalized目錄下
cp referenceDB.800 RD_normalized
# 進(jìn)入RD_normalized目錄
cd RD_normalized
# 將RD_normalized目錄下新生成的sex_1結(jié)尾的文件名，以絕對(duì)路徑的形式寫(xiě)入list.txt中
ls -R `pwd`/*sex_1 > list.txt
# 新建exclude_list文件
touch exclude_list
# 拷貝數(shù)變異區(qū)域的確定
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斷點(diǎn)計(jì)算方法
1)每個(gè)樣品分別進(jìn)行matchclips2計(jì)算：
/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
結(jié)果文件：
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勋眯、區(qū)域長(zhǎng)度：3' break point - 5' break point
6、UN下梢，堿基數(shù)量This number tells how many bases are repeated at the two break points
7客蹋、RD:n1;n2;n3:s
2)合并結(jié)果：
/home/software/matchclips2/cnvtable -L 10000 -cnvf cnvf.txt -O 0.5 -o overlap.txt
cnvf.txt為含有上一步結(jié)果文件名稱的txt文件，一列一個(gè)
-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循環(huán)
touch matchclips2.sh
##運(yùn)行
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
#將生成文件名列到一個(gè)文件中
find -name '*cnv' -exec basename {} \; > cnvf.txt
#進(jìn)行合并
/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讶坯；低版本的需要通過(guò)-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.
有群體時(shí)辆琅，建議將所有個(gè)體的.bcf文件合并后，再進(jìn)行g(shù)enotyping这刷。
運(yùn)行：詳情參考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轉(zhuǎn)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(基因分型婉烟，每個(gè)樣本都做)
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）
反轉(zhuǎn)（Inversion, INV）
染色體內(nèi)部易位（Intra-chromosomal Translocation, ITX）
染色體間易位（Inter-chromosomal Translocation, CTX)
4、CNVnator(拉跨）
檢測(cè)出的片段長(zhǎng)度很長(zhǎng)暇屋，假陽(yáng)性結(jié)果率較低似袁，還算比較準(zhǔn)確的一個(gè)軟件
1)提取mapping信息
home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -tree BC-448.sorted.addhead.markdup.bam -unique
提取比對(duì)后的reads，構(gòu)建樹(shù)咐刨，如果程序之前運(yùn)行失敗昙衅，已經(jīng)生成了一個(gè)root文件，在下次重新運(yùn)行時(shí)一定要?jiǎng)h除該root文件定鸟，如果不刪除绒尊，新的分析結(jié)果會(huì)追加到錯(cuò)誤的root文件中，影響后續(xù)分析仔粥。另外婴谱，可以添加 -chrom 函數(shù)指定染色體 如 -chrom 1 2 代表選擇1，2號(hào)染色體躯泰。
2)生成質(zhì)量分布圖HISTOGRAM
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -d -his 100
# 100指的是bin size谭羔，可以通過(guò)-eval參數(shù)進(jìn)行篩選，也可以根據(jù)經(jīng)驗(yàn)值進(jìn)行確定麦向，一般測(cè)序深度20-30x選取bin size大小100瘟裸，2-3x選取500，100x選取30诵竭。
# -d指定目錄话告，內(nèi)部存放給染色體的fasta文件兼搏，該參數(shù)指針對(duì)報(bào)錯(cuò)信息顯示不能解析基因組文件，此時(shí)需要指定參考基因組的位置沙郭，且各染色體需要拆分成單獨(dú)的fasta文件（目錄下有其文件也可以佛呻，只要有所有染色體的序列就好，軟件會(huì)自動(dòng)識(shí)別）
3)生成統(tǒng)計(jì)結(jié)果
/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在進(jìn)行變異檢測(cè)時(shí)病线，以提供的bin size對(duì)整個(gè)基因組進(jìn)行切割吓著，之后按照RD(read-depth)為基準(zhǔn)進(jìn)行cnv的檢測(cè)。
注意：從第二步開(kāi)始設(shè)置的bin size就不能變了
5)變異檢出
/home/software/CNVnator_v0.4.1/src/cnvnator -root BC-448.root -call 100 > BC-448-cnvout.txt
運(yùn)行之后輸出結(jié)果如下：
#第一列變異類型
#第二列位點(diǎn)信息
#第三列CNV大小
#第四列為標(biāo)準(zhǔn)化參數(shù)normalized_RD
#第5-8列為e-value值送挑，其中第五列越小绑莺，說(shuō)明結(jié)果越準(zhǔn)確
#第九列q0質(zhì)量值
6)轉(zhuǎn)為vcf格式
/home/software/CNVnator_v0.4.1/src/cnvnator2VCF.pl cnv.call.txt > cnv.vcf
5、LUMPY（挺牛的一個(gè)軟件）
每種算法都要其優(yōu)勢(shì)和不足之處惕耕，綜合運(yùn)用多種策略有助于提高檢測(cè)的靈敏度纺裁，lumpy就是這樣一款軟件，集合了read-pair司澎，split-read对扶，read-depth, 等多種策略來(lái)預(yù)測(cè)CNV
lumpy, svtools, svtyper安裝包在/home/sll/miniconda3/envs/python2.7/bin下
LUMPY的安裝使用python2.7的環(huán)境
同理，svtools和svtyper也使用python2.7的環(huán)境
安裝：conda install -c bioconda lumpy-sv
1. mapping
推薦采用bwa-mem算法將雙端序列比對(duì)到參考基因組上惭缰，為了加快運(yùn)行速度浪南，這里用samblaster軟件進(jìn)行markduplicate, 用法如下
samblaster --excludeDups \
--addMateTags \
--maxSplitCount 2 \
--minNonOverlap 20 \
samtools view -Sb - > sample.bam
2. extract discordant paired-end alignments
鏈接：http://www.reibang.com/p/1844870f3b25
discordant reads指的是R1和R2端比對(duì)之間的距離超過(guò)了期望的插入片段長(zhǎng)度或者比對(duì)到了不同鏈的reads,
samtools view -b -F 1294 \
sample.bam \
> sample.discordants.unsorted.bam
-F 1294用來(lái)提取不一致的聯(lián)配，用samtools flags 1294可以發(fā)現(xiàn)1294表示"PROPER_PAIR,UNMAP,MUNMAP,SECONDARY,DUP"漱受，帶上-F意味著以上這些標(biāo)記在我們篩選的聯(lián)配記錄中都不會(huì)出現(xiàn)络凿，也就意味著篩選的記錄要符合下面要求
不能是PROPER_PAIR: 就是比對(duì)工具認(rèn)為都正確比對(duì)到基因組上，在同一條染色體昂羡，在同一條鏈的情況絮记，常見(jiàn)的就是83,147和99,163
不能是UNMAP和MUNMAP，也就是配對(duì)的短讀至少有一個(gè)能夠比對(duì)到參考基因組上
也不能是SECONDARY虐先， 也就是他必須是主要聯(lián)配
光學(xué)重復(fù)怨愤，DUP, 就更加不能要了
3. extract split-reads alignments
split-reads指的是覆蓋了斷裂點(diǎn)的單端reads,這些reads根據(jù)斷裂點(diǎn)拆分成subreads后可以正確的比多到參考基因組上。在軟件的安裝目錄蛹批，自帶了一個(gè)名為extractSplitReads_BwaMem的腳本撰洗，用于提取split-reads, 用法如下
samtools view -h sample.bam \
| scripts/extractSplitReads_BwaMem -i stdin \
| samtools view -Sb - \
> sample.splitters.unsorted.bam
4. sort bams
軟件要求輸入的bam文件必須是排序之后的文件，所以對(duì)提取的兩個(gè)子bam進(jìn)行排序腐芍，用法如下
samtools sort \
sample.discordants.unsorted.bam \
sample.discordants
samtools sort \
sample.splitters.unsorted.bam \
sample.splitters
5. run lumpy
lumpyexpress是lumpy的一個(gè)封裝腳本差导，使用起來(lái)更加方便，基本用法如下
lumpyexpress \
-B sample.bam \
-S sample.splitters.bam \
-D sample.discordants.bam \
-o sample.vcf
6. genotype
檢測(cè)到的CNV, 可以用svtyper這個(gè)軟件預(yù)測(cè)在樣本中的分型結(jié)果猪勇，用法如下
svtyper軟件需要另外安裝
conda install -c bioconda svtyper
svtyper \
-B sample.bam \
-S sample.splitters.bam \
-i sample.vcf
> sample.gt.vcf
選擇清除與環(huán)境適應(yīng)性位點(diǎn)挖掘--Vst分析
Vst分析是類似于Fst的一個(gè)指標(biāo)设褐，用來(lái)衡量群體間每個(gè)CNVR差異大小的統(tǒng)計(jì)量，計(jì)算方法為Vst＝（Vt－Vs）/Vt，Vt表示所有樣本該區(qū)域拷貝數(shù)大小的方差助析，Vs表示兩個(gè)群體各自的方差根據(jù)各自群體大小加權(quán)之后的值犀被。Vst的值介于0-1之間，值越大表示群體間該區(qū)域拷貝數(shù)變異差異越大外冀，反之則越小寡键。
EXCEL計(jì)算.png
