工具安装及试用总结:对WES数据做germline的CNV calling(非癌症的疾病研究)(未完。。。)

发现好多CNV calling 工具都好古早。。。安装和试用时关于版本的问题调试比较多。。。所以想把自己遇到的报错贴出来,方便后人debug

1、 EXCAVATOR2

原理及软件介绍:使用EXCAVATOR2检测WES的CNV
https://mp.weixin.qq.com/s/WcbCXq9Y7FGtvZXS7-HCEA
上面这个链接基本简单地做了说明,下面就简单记录一下我自己在安装和使用上遇到的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

安装时遇到的问题:

R, SAMtools, Perl基本都是服务器上早就装好的,版本一般都不低,所以没什么问题。但是在装Hmisc这个R包的时候:
我用的R-3.5装的时候,给我报错 latticeExtra 这个包 not available,说这个 latticeExtra 需要的版本更高。。。我换了R-3.6装,在装 latticeExtra 的时候说 jpeg 的包有问题。。。逐个试了很多R包安装方法,都不行。

解决方法:

由于我们是实验室用一个服务器,最开始默认的R和配套的lib可能比较旧,或者有各种只有root权限才能修改的东西。所以这里可以重新用conda 建一个环境,下载一个靠谱的R,重头安装一遍你需要的包。
或者,组里有同学的R和lib可以完成这个包的安装,就直接引用到你自己的环境变量里吧!哈哈哈哈哈哈

alias R=‘谢谢大哥的R路径’ 
export R_LIBS_SITE="谢谢大哥的R lib路径:$R_LIBS_SITE"

运行的必要条件:

首先就是一定要记得在这个软件的解压后路径下运行命令!
因为这些perl是直接读取你运行命令的这个位置,然后会在后面用这个路径的字符串编辑一些新的路径,
所以一定要在这个软件的解压后路径下运行命令!(其实就是我自己被蠢到过。。。)

运行时遇到的问题:

在运行第一步TargetPerla.pl的时候,遇到如下报错:

~/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

主要问题是加载不到这个:libpng12.so.0 (libpng15它都不认的,好像只认这一个版本。。。)

解决方法:

单独下了这个lib文件的相关文件,直接放到原来的lib路径下,行不通。。。
conda install libpng=1.2就OK啦!如果这个conda的lib路径原来不在环境变量里,新加进去就OK了:

export LD_LIBRARY_PATH=确定装有libpng12.so.0的路径:$LD_LIBRARY_PATH

2、CoNIFER

使用环境:python2.7

准备probes文件

只认chr1-22&XY,要把chrM去掉,否则如下报错,

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'

运行中的可能错误1

conifer.py的第564行有个“f._has_Index()”,随着pysam包的版本不同,该命令的写法不同,可以都试一下
https://sourceforge.net/p/conifer/discussion/general/thread/d2fbc181/?limit=25
可以通过conda list先确定一下你的pysam版本,然后修改到对应的。
搞pysam的时候我被玄学到了,先是镜像的问题,只装上了0.6,上述任何版本的修改都没用。。。。
修改镜像之后装了最新的0.16,失败。。。
随手装了0.9,使用的时候import失败。。。
然后换成0.8,安装和import成功,改成“f.has_Index()”可以成功运行

运行中的可能错误2

再有是关于tables,由于conifer实在是太古早了,其中的语法都是tables2.0的版本,会有如下各种报错。。。

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``

全网看了一圈,已经下载不到tables2.0了,所以还是用的现成装的tables3.5,边改边test,根据这个网页对应把旧的语法改成新的就行。http://www.pytables.org/MIGRATING_TO_3.x.html?highlight=creategroup 后面的call步骤也是要记得修改新版语法。

运行中的可能错误3

analyse有一个可能的报错“IndexError: boolean index did not match indexed array along dimension 0; dimension is 24661 but corresponding boolean dimension is 24660”,可能是因为numpy版本的问题,解决方法是
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

运行中的可能错误4

plotcalls画图的部分提醒补安装了matplotlib,出现了报错:

/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测试数据跑通

需要额外关注和计算的参数

analyse中的--svd参数,官网教程给了说明,应该根据你的样本数和数据方差去选择合适的svd数,具体可以看文献理解

3.XHMM

官方说明:http://atgu.mgh.harvard.edu/xhmm/tutorial.shtml
(曾经打开过,并下载到了安装包“statgen-xhmm-998f7c405974.zip”,然而在正式要用的时候,我翻不翻墙都没能再打开这个链接。。。)
好在文章有protocol:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4065038/
【发现可以用的tutorial了:[https://statgen.bitbucket.io/xhmm/tutorial.html]https://statgen.bitbucket.io/xhmm/tutorial.html)】
但是在make编译的时候有一大堆报错,在网上兜了一圈,好像都说make容易有问题,在bioconda上有现成的packagehttps://anaconda.org/bioconda/xhmm
conda install -c bioconda xhmm一句话搞定安装~

首先需要GATK的DepthOfCoverage来计算一个覆盖深度的值,但是这个工具是属于GATK3的,GATK4从4.1.6版本才开始重新复原这个tool,而我自己手头是4.1.3的版本,所以就用conda新建了一个环境,专门装了一个gatk3.8conda create -n gatk3 -c bioconda gatk
安装后需要注册一下,解决操作参考:https://zhuanlan.zhihu.com/p/129858566
由于我们实验室有几个服务器,配置略有不同,在某个服务器中运行DepthOfCoverage的过程中发现了如下报错:

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文件,用了注册时的jar不大行,所以去找了3.8.1的jar进行替换,尝试成功,可以正常运行。
gatk3.8及以前的版本可以在google云上找到:https://console.cloud.google.com/storage/browser/gatk-software/package-archive/gatk;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

之后按tutorial的步骤一步步操作即可。

文章的protocol中说关于filter的具体参数可以用后面作图的protocol把一些值的范围都找出来,但是这里暂时受限于Plink/Seq的locdb参考数据下载不下来,这个网址打不开http://atgu.mgh.harvard.edu/plinkseq/resources.shtml
如果有可以下载的路径求分享!
但是这个步骤是optional的,所以最后使用的时候我选择了跳过。。。用tutorial的参考值进行的后续分析()

4、CANOES

使用说明可以在这里找到https://github.com/ShenLab/CANOES
其中需要的软件工具里,GATK的GCContentByInterval,也是个在GATK4(至少4.1.0.3)里没有的,所以上面创建的gatk3.8环境又有用了!
软件原理是基于每批次WES的背景值进行分析,所以要按上机批次进行分析,最好30个以上一批,20个以上也行(这里的用法我之前理解错了,感谢DXR师姐提点!)
该软件工具只能处理常染色体,需要分析X染色体的话可以参考下面这篇文章的方法部分,通过修改R包加入了X染色体上的分析:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7418612/
canoes.reads.txt文件的第一列处理完默认是chr1,chr2,...,chr21,chr22但是在后面的分析中,只认,1,2,...,21,22,所以中间可以对这个文件进行一下预处理。
在 # call CNVs with the Viterbi algorithm#这步中发现Viterbi 这个function中的viterbi.pointers[i, ] <- apply(temp.matrix, 2, which.max)会报错,发现是viterbi.matrix[i, ] <- apply(temp.matrix, 2, max)中,如果temp.matrix中偶尔会有个别NaN值,会被认为是最大值,从而后面的全部变成了NaN。这里想到的解决方法是在449行加上一句temp.matrix[is.na(temp.matrix)] <- (-Inf)直接把Na替换成负无穷,这样就不会被误认为是max了。后面可以正常运行了。

5、CODEX2

说明文档
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值可根据中间的碎石图选择范围,软件推荐1-10。总体运行顺畅没有显著问题

6、HMZDelFinder

R包和sample:https://github.com/BCM-Lupskilab/HMZDelFinder/
自己的WES数据可利用官方给的函数calcRPKMsFromBAMs来制作(4是函数中apply家族需要使用的核数)
calcRPKMsFromBAMs(bedFile, bamdir, sampleNames, rpkmDir,4)
该函数运行过程中容易出现一个问题,调用data.table的fread()读取bed文件时的报错:

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

尝试修改了文件换行符、data.table的版本数,都没能解决。了解了fread()的功能主要是快速读取大文件,私以为此步骤的耗时速度并非关键因素,所以将官方R文件中的第98行的bed <- fread(bedFile)修改为bed <- read.table(bedFile,header=F)。解决问题~
后面读取vcf和RPKM文件的步骤里都用到了fread(),在我自己的R里都会有问题,所以都改成了read.table() + as.data.table()进行表格的格式转换。反正这里经历了比较痛苦的逐行debug过程。。。

7、ExomeDepth

该R包运行问题较少,多批次多样本可进行循环处理,具体操作可参考://www.greatytc.com/p/a650a9d9a861

8、CONTRA

上一条引用的博主小姐姐用过这个软件,这里直接引用一下://www.greatytc.com/p/f23cc2c4b45d
软件论文:https://academic.oup.com/bioinformatics/article/28/10/1307/212453
说明文档:http://contra-cnv.sourceforge.net/
软件文章的讨论部分说到该软件对数据要求不高,也没提到批次效应的问题,所以应该可以把所有样本都丢进去做。
蓝鹅,我连软件都没下载下来。。。https://sourceforge.net/projects/contra-cnv/files/CONTRA.V2.0/

后面的Flag:

ADTEx

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,837评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,551评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,417评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,448评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,524评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,554评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,569评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,316评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,766评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,077评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,240评论 1 343
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,912评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,560评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,176评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,425评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,114评论 2 366
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,114评论 2 352

推荐阅读更多精彩内容