發(fā)現(xiàn)好多CNV calling 工具都好古早钝凶。。唁影。安裝和試用時(shí)關(guān)于版本的問(wèn)題調(diào)試比較多耕陷。。据沈。所以想把自己遇到的報(bào)錯(cuò)貼出來(lái)哟沫，方便后人debug

1、 EXCAVATOR2

原理及軟件介紹：使用EXCAVATOR2檢測(cè)WES的CNV
https://mp.weixin.qq.com/s/WcbCXq9Y7FGtvZXS7-HCEA
上面這個(gè)鏈接基本簡(jiǎn)單地做了說(shuō)明锌介，下面就簡(jiǎn)單記錄一下我自己在安裝和使用上遇到的bug和解決方法吧~

安裝的必要條件：

EXCAVATOR2 was conceived for running on 64-bit UNIX desktop machines with at least 4 CPUs and 4 GB RAM.
In order to work properly EXCAVATOR2 needs R (version≥2.14.0) and the Hmisc library (R package), SAMtools(version≥0.1.17),andPerl(version≥5.8.8)tobecorrectlyinstalledonyoursystem

安裝時(shí)遇到的問(wèn)題：

R, SAMtools, Perl基本都是服務(wù)器上早就裝好的嗜诀，版本一般都不低，所以沒(méi)什么問(wèn)題孔祸。但是在裝Hmisc這個(gè)R包的時(shí)候：
我用的R-3.5裝的時(shí)候隆敢，給我報(bào)錯(cuò) latticeExtra 這個(gè)包 not available，說(shuō)這個(gè) latticeExtra 需要的版本更高崔慧。拂蝎。。我換了R-3.6裝惶室，在裝 latticeExtra 的時(shí)候說(shuō) jpeg 的包有問(wèn)題温自。。皇钞。逐個(gè)試了很多R包安裝方法悼泌，都不行。

解決方法：

由于我們是實(shí)驗(yàn)室用一個(gè)服務(wù)器鹅士，最開(kāi)始默認(rèn)的R和配套的lib可能比較舊券躁，或者有各種只有root權(quán)限才能修改的東西。所以這里可以重新用conda 建一個(gè)環(huán)境掉盅，下載一個(gè)靠譜的R也拜，重頭安裝一遍你需要的包。
或者趾痘，組里有同學(xué)的R和lib可以完成這個(gè)包的安裝慢哈，就直接引用到你自己的環(huán)境變量里吧！哈哈哈哈哈哈

alias R=‘謝謝大哥的R路徑’ 
export R_LIBS_SITE="謝謝大哥的R lib路徑:$R_LIBS_SITE"

運(yùn)行的必要條件：

首先就是一定要記得在這個(gè)軟件的解壓后路徑下運(yùn)行命令永票！
因?yàn)檫@些perl是直接讀取你運(yùn)行命令的這個(gè)位置卵贱，然后會(huì)在后面用這個(gè)路徑的字符串編輯一些新的路徑滥沫，
所以一定要在這個(gè)軟件的解壓后路徑下運(yùn)行命令！（其實(shí)就是我自己被蠢到過(guò)键俱。兰绣。。）

運(yùn)行時(shí)遇到的問(wèn)題：

在運(yùn)行第一步TargetPerla.pl的時(shí)候编振，遇到如下報(bào)錯(cuò)：

~/software/EXCAVATOR2_Package_v1.1.2/lib/OtherLibrary/bigWigAverageOverBed: error while loading shared libraries: libpng12.so.0: cannot open shared object file: No such file or directory
Error in file(file, "rt") : cannot open the connection
Calls: read.table -> file
In addition: Warning message:
In file(file, "rt") :
cannot open file '~/EXCAVATOR2_Package_v1.1.2/data/targets/hg19/AJTK_w10000/MAP/Mapout.txt': No such file or directory
Execution halted

主要問(wèn)題是加載不到這個(gè)：libpng12.so.0 （libpng15它都不認(rèn)的缀辩，好像只認(rèn)這一個(gè)版本。踪央。臀玄。）

解決方法：

單獨(dú)下了這個(gè)lib文件的相關(guān)文件，直接放到原來(lái)的lib路徑下畅蹂，行不通健无。。液斜。
用conda install libpng=1.2就OK啦累贤！如果這個(gè)conda的lib路徑原來(lái)不在環(huán)境變量里，新加進(jìn)去就OK了：

export LD_LIBRARY_PATH=確定裝有l(wèi)ibpng12.so.0的路徑:$LD_LIBRARY_PATH

2旗唁、CoNIFER

使用環(huán)境：python2.7

準(zhǔn)備probes文件

只認(rèn)chr1-22&XY畦浓，要把chrM去掉，否則如下報(bào)錯(cuò)检疫，

Traceback (most recent call last):
  File "conifer.py", line 682, in <module>
    args.func(args)
  File "conifer.py", line 545, in CF_bam2RPKM
    probes = cf.loadProbeList(probe_fn)
  File "/picb/dermatogenomics/chenjieyi/software/conifer_v0.2.2/conifer_functions.py", line 96, in loadProbeList
    probes.append({'probeID': probeID, 'chr':chrStr2Int(row['chr']),'start':int(row['start']),'stop':int(row['stop']), 'name':row['name']})
  File "/picb/dermatogenomics/chenjieyi/software/conifer_v0.2.2/conifer_functions.py", line 58, in chrStr2Int
    return int(chr)
ValueError: invalid literal for int() with base 10: 'M'

運(yùn)行中的可能錯(cuò)誤1

conifer.py的第564行有個(gè)“f._has_Index()”，隨著pysam包的版本不同祷嘶，該命令的寫(xiě)法不同屎媳，可以都試一下
https://sourceforge.net/p/conifer/discussion/general/thread/d2fbc181/?limit=25
可以通過(guò)conda list先確定一下你的pysam版本，然后修改到對(duì)應(yīng)的论巍。
搞pysam的時(shí)候我被玄學(xué)到了烛谊，先是鏡像的問(wèn)題，只裝上了0.6嘉汰，上述任何版本的修改都沒(méi)用丹禀。。鞋怀。双泪。
修改鏡像之后裝了最新的0.16，失敗密似。焙矛。。
隨手裝了0.9残腌，使用的時(shí)候import失敗村斟。贫导。。
然后換成0.8蟆盹，安裝和import成功孩灯，改成“f.has_Index()”可以成功運(yùn)行

運(yùn)行中的可能錯(cuò)誤2

再有是關(guān)于tables，由于conifer實(shí)在是太古早了逾滥，其中的語(yǔ)法都是tables2.0的版本钱反，會(huì)有如下各種報(bào)錯(cuò)。匣距。面哥。

balabala Error 'openFile'
AttributeError: 'File' object has no attribute 'createGroup'
AttributeError: 'File' object has no attribute 'createTable'
tables.exceptions.NoSuchNodeError: group ``/`` does not have a child named ``_f_getChild``

全網(wǎng)看了一圈，已經(jīng)下載不到tables2.0了毅待，所以還是用的現(xiàn)成裝的tables3.5尚卫，邊改邊test，根據(jù)這個(gè)網(wǎng)頁(yè)對(duì)應(yīng)把舊的語(yǔ)法改成新的就行尸红。http://www.pytables.org/MIGRATING_TO_3.x.html?highlight=creategroup 后面的call步驟也是要記得修改新版語(yǔ)法吱涉。

運(yùn)行中的可能錯(cuò)誤3

analyse有一個(gè)可能的報(bào)錯(cuò)“IndexError: boolean index did not match indexed array along dimension 0; dimension is 24661 but corresponding boolean dimension is 24660”，可能是因?yàn)閚umpy版本的問(wèn)題外里，解決方法是
The error is in line 142 of conifer.py, instead of:

rpkm = RPKM_data[start_probeID:stop_probeID,:]

it should be:

rpkm = RPKM_data[start_probeID-1:stop_probeID,:]

參考：https://github.com/UBC-Stat-ML/conifer/issues/26

運(yùn)行中的可能錯(cuò)誤4

plotcalls畫(huà)圖的部分提醒補(bǔ)安裝了matplotlib怎爵，出現(xiàn)了報(bào)錯(cuò)：

/data/dermatogenomics4/software/anaconda3/envs/py27/lib/python2.7/site-packages/matplotlib/pyplot.py:522: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)

解決方法：在conifer.py 的 line 460 附近，在import matplotlib后面一行加上matplotlib.rcParams.update({'figure.max_open_warning': 0})

OK_{測(cè)試數(shù)據(jù)跑通}

需要額外關(guān)注和計(jì)算的參數(shù)

analyse中的--svd參數(shù)盅蝗，官網(wǎng)教程給了說(shuō)明鳖链，應(yīng)該根據(jù)你的樣本數(shù)和數(shù)據(jù)方差去選擇合適的svd數(shù)，具體可以看文獻(xiàn)理解

3.XHMM

官方說(shuō)明：http://atgu.mgh.harvard.edu/xhmm/tutorial.shtml
（曾經(jīng)打開(kāi)過(guò)墩莫，并下載到了安裝包“statgen-xhmm-998f7c405974.zip”芙委，然而在正式要用的時(shí)候，我翻不翻墻都沒(méi)能再打開(kāi)這個(gè)鏈接狂秦。灌侣。。）
好在文章有protocol：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4065038/
【發(fā)現(xiàn)可以用的tutorial了：[https://statgen.bitbucket.io/xhmm/tutorial.html]（https://statgen.bitbucket.io/xhmm/tutorial.html)】
但是在make編譯的時(shí)候有一大堆報(bào)錯(cuò)裂问，在網(wǎng)上兜了一圈侧啼，好像都說(shuō)make容易有問(wèn)題，在bioconda上有現(xiàn)成的packagehttps://anaconda.org/bioconda/xhmm
堪簿，conda install -c bioconda xhmm一句話搞定安裝~

首先需要GATK的DepthOfCoverage來(lái)計(jì)算一個(gè)覆蓋深度的值痊乾，但是這個(gè)工具是屬于GATK3的，GATK4從4.1.6版本才開(kāi)始重新復(fù)原這個(gè)tool戴甩，而我自己手頭是4.1.3的版本符喝，所以就用conda新建了一個(gè)環(huán)境，專門(mén)裝了一個(gè)gatk3.8conda create -n gatk3 -c bioconda gatk
安裝后需要注冊(cè)一下甜孤，解決操作參考：https://zhuanlan.zhihu.com/p/129858566
由于我們實(shí)驗(yàn)室有幾個(gè)服務(wù)器协饲，配置略有不同畏腕，在某個(gè)服務(wù)器中運(yùn)行DepthOfCoverage的過(guò)程中發(fā)現(xiàn)了如下報(bào)錯(cuò)：

ERROR StatusLogger Unable to create class org.apache.logging.log4j.core.impl.Log4jContextFactory specified in jar:file:/[conda env]/opt/gatk-3.8/GenomeAnalysisTK.jar!/META-INF/log4j-provider.properties
ERROR StatusLogger Log4j2 could not find a logging implementation. Please add log4j-core to the classpath. Using SimpleLogger to log to the console...

查了一圈是要想辦法替換conda env中的jar文件，用了注冊(cè)時(shí)的jar不大行茉稠，所以去找了3.8.1的jar進(jìn)行替換描馅，嘗試成功，可以正常運(yùn)行而线。
gatk3.8及以前的版本可以在google云上找到：https://console.cloud.google.com/storage/browser/gatk-software/package-archive/gatk;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

之后按tutorial的步驟一步步操作即可铭污。

文章的protocol中說(shuō)關(guān)于filter的具體參數(shù)可以用后面作圖的protocol把一些值的范圍都找出來(lái)，但是這里暫時(shí)受限于Plink/Seq的locdb參考數(shù)據(jù)下載不下來(lái)膀篮，這個(gè)網(wǎng)址打不開(kāi)http://atgu.mgh.harvard.edu/plinkseq/resources.shtml
如果有可以下載的路徑求分享嘹狞！
但是這個(gè)步驟是optional的，所以最后使用的時(shí)候我選擇了跳過(guò)誓竿。磅网。。用tutorial的參考值進(jìn)行的后續(xù)分析（）

4筷屡、CANOES

使用說(shuō)明可以在這里找到https://github.com/ShenLab/CANOES
其中需要的軟件工具里涧偷，GATK的GCContentByInterval，也是個(gè)在GATK4（至少4.1.0.3）里沒(méi)有的毙死，所以上面創(chuàng)建的gatk3.8環(huán)境又有用了燎潮！
軟件原理是基于每批次WES的背景值進(jìn)行分析，所以要按上機(jī)批次進(jìn)行分析扼倘，最好30個(gè)以上一批确封，20個(gè)以上也行（這里的用法我之前理解錯(cuò)了，感謝DXR師姐提點(diǎn)０π俊）
該軟件工具只能處理常染色體隅肥，需要分析X染色體的話可以參考下面這篇文章的方法部分，通過(guò)修改R包加入了X染色體上的分析：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7418612/
canoes.reads.txt文件的第一列處理完默認(rèn)是chr1,chr2,...,chr21,chr22但是在后面的分析中袄简，只認(rèn)，1,2,...,21,22泛啸，所以中間可以對(duì)這個(gè)文件進(jìn)行一下預(yù)處理绿语。
在 # call CNVs with the Viterbi algorithm#這步中發(fā)現(xiàn)Viterbi 這個(gè)function中的viterbi.pointers[i, ] <- apply(temp.matrix, 2, which.max)會(huì)報(bào)錯(cuò)，發(fā)現(xiàn)是viterbi.matrix[i, ] <- apply(temp.matrix, 2, max)中候址，如果temp.matrix中偶爾會(huì)有個(gè)別NaN值吕粹，會(huì)被認(rèn)為是最大值，從而后面的全部變成了NaN岗仑。這里想到的解決方法是在449行加上一句temp.matrix[is.na(temp.matrix)] <- (-Inf)直接把Na替換成負(fù)無(wú)窮匹耕，這樣就不會(huì)被誤認(rèn)為是max了。后面可以正常運(yùn)行了荠雕。

5稳其、CODEX2

說(shuō)明文檔
http://htmlpreview.github.io/?https://github.com/yuchaojiang/CODEX2/blob/master/demo/CODEX2.html
https://github.com/yuchaojiang/CODEX2
R<3.5
source("https://bioconductor.org/biocLite.R")
R≥3.5
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install()
K值可根據(jù)中間的碎石圖選擇范圍驶赏，軟件推薦1-10〖染希總體運(yùn)行順暢沒(méi)有顯著問(wèn)題

6煤傍、HMZDelFinder

R包和sample：https://github.com/BCM-Lupskilab/HMZDelFinder/
自己的WES數(shù)據(jù)可利用官方給的函數(shù)calcRPKMsFromBAMs來(lái)制作（4是函數(shù)中apply家族需要使用的核數(shù)）
calcRPKMsFromBAMs(bedFile, bamdir, sampleNames, rpkmDir,4)
該函數(shù)運(yùn)行過(guò)程中容易出現(xiàn)一個(gè)問(wèn)題，調(diào)用data.table的fread()讀取bed文件時(shí)的報(bào)錯(cuò)：

Error in fread(bedFile) : 
  Internal error: invalid head position. jump=0, headPos=0x7f60a986513f, thisJumpStart=0x7f60a975e000, sof=0x7f60a975e000

嘗試修改了文件換行符嘱蛋、data.table的版本數(shù)蚯姆，都沒(méi)能解決。了解了fread()的功能主要是快速讀取大文件洒敏，私以為此步驟的耗時(shí)速度并非關(guān)鍵因素龄恋，所以將官方R文件中的第98行的bed <- fread(bedFile)修改為bed <- read.table(bedFile,header=F)。解決問(wèn)題~
后面讀取vcf和RPKM文件的步驟里都用到了fread()凶伙，在我自己的R里都會(huì)有問(wèn)題郭毕，所以都改成了read.table() + as.data.table()進(jìn)行表格的格式轉(zhuǎn)換。反正這里經(jīng)歷了比較痛苦的逐行debug過(guò)程镊靴。铣卡。。

7偏竟、ExomeDepth

該R包運(yùn)行問(wèn)題較少煮落，多批次多樣本可進(jìn)行循環(huán)處理，具體操作可參考：http://www.reibang.com/p/a650a9d9a861

8踊谋、CONTRA

上一條引用的博主小姐姐用過(guò)這個(gè)軟件蝉仇，這里直接引用一下：http://www.reibang.com/p/f23cc2c4b45d
軟件論文：https://academic.oup.com/bioinformatics/article/28/10/1307/212453
說(shuō)明文檔：http://contra-cnv.sourceforge.net/
軟件文章的討論部分說(shuō)到該軟件對(duì)數(shù)據(jù)要求不高，也沒(méi)提到批次效應(yīng)的問(wèn)題殖蚕，所以應(yīng)該可以把所有樣本都丟進(jìn)去做轿衔。
藍(lán)鵝，我連軟件都沒(méi)下載下來(lái)睦疫。害驹。。https://sourceforge.net/projects/contra-cnv/files/CONTRA.V2.0/

后面的Flag：

ADTEx