2021-05-25 picard和beagle软件使用过程中遇到的各种问题

声明:因为本人使用过程中遇到了很多原来人碰得到的问题,放防止以后遗忘,所以做一下记录,因为前总结的很到位,所以直接引用到自己日常记录,供日后查阅。

主要引用为:【原创文章】辣椒BSA分析对数据进行质控、去重、索引等预处理的几个姿势 (https://www.baishujun.com/archives/7191.html)!

去除PCR扩增重复序列

由于构建文库的时候进行了pcr扩增,然后进行测序,这么测序出来的数据存在很多重复序列,在这里需要去除,然后建立索引,那么这个预处理工作就结束了。

遇到的问题最多的也是在这里。

按照之前的资料,对于PE测序数据一般用picard,效果好,SE可以用samtools或picard去重。

所以我首先用picard的MarkDuplicates进行了去重,遇到N个错误!

1,临时目录不够

错误提示:

java运行问题,原因/tmp空间太小造成的,有两个解决办法。

(1),给picard设置特定的tmp目录 (原作者推荐

picard MarkDuplicates I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt TMP_DIR=/www/tmp

本人项目中也是用特定目录解决, 在picard运行命令之前,加一行命令:

export _JAVA_OPTIONS=-Djava.io.tmpdir=./tmp

picard -Xmx32g MarkDuplicates I=S99.new1.bam O=S99.rmdup.bam CREATE_INDEX=ture REMOVE_DUPLICATES=ture M=S99.marked_dup_metrics.txt

(2),为了一劳永逸的解决这个问题,我直接挂载了一个1t的分区给/tmp(原作者)


给/tmp,/home和根目录/分别挂载1t,给数据目录/www挂载20t,还留几t机动调配。

#####此处还有一点项目中的查询结果针对 java.lang.OutOfMemoryError: GC overhead limit exceeded

如果在垃圾收集中花费太多的时间太少,GC会抛出这个异常。GC上的CPU时间占用了98%,只有不到2%的堆被恢复。

此功能旨在防止应用程序长时间运行,而由于堆太小,进行很少或没有进度。

您可以使用命令行选项关闭此功能  -XX:-UseGCOverheadLimit (我在使用beagle尝试使用过这个但是没有解决问题,其他原因)。

2,java堆大小设置过低

错误提示:

面对“Java heap space”的问题,主要是由于我的测序数据大,占用内存大,但是java程序默认的堆大小分配不够,所以出现这个问题,我们可以给当前运行的java程序设定一个内存值,从而覆盖默认设置。

picard -Xmx40g MarkDuplicates I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt (Xmx40g,设置JVM最大可用内存为40G.)

3,“Insert size out of range”错误

关于这个错误,从字面意思上来看是插入大小超过了范围,并且我发现基本都是在检测重复完成之后才出现,这个错误应当和程序读取reads的大小或者java写入磁盘文件大小相关方面的限制有关,但是没有找到具体原因,从错误提示也可以看出,这是在确认sam数据时出现的问题,那么我们可以采用宽松一点的审计策略,就可以绕过这个“Insert size out of range”错误。

picard -Xmx40g MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt

我用的是conda环境,所以上面的命令是我运行picard进行序列去重成功的命令。

如果不想使用conda的picard,那么我们也可以直接调用java包,在picard命令中已经包含了当前picard的jar包的地址,我们可以修改命令如下:

java -Xmx40g -jar /www/soft/miniconda2/envs/pepper/share/picard-2.20.3-0/picard.jar MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt

使用了VALIDATION_STRINGENCY=LENIENT参数以后

我们可以看到大量的提示

 其它去重软件

经过这一轮之后,整理发现,还有好几个很优秀的去重软件,

1,Sambamba

地址:http://lomereiter.github.io/sambamba/

原文描述如下:

Marks (by default) or removes duplicate reads. For determining whether a read is a duplicate or not, the same criteria as in Picard are used.

这是我推荐的替代picard的软件

她筛选reads是否重复的标准和picard是一样的(For determining whether a read is a duplicate or not, the same criteria as in Picard are used)。

关键是这个软件支持多线程,速度就快多了,她也支持是标记或者删除重复序列,并且自动建立bai索引文件

该软件还有其它的sort,index,merge等功能,详细地址为,http://lomereiter.github.io/sambamba/docs/sambamba-markdup.html#OPTIONS,

需要代理才能打开。

sambamba markdup -p -t 16 S347.merged.sorted.viewed.bam S347.sorted.viewed.markdupba.bam

-p显示过程,-t线程数16.

对于一个46.70GB的bam文件进行去重,picard耗时78.11分钟,而sambamba仅仅用时33分钟,快了一半吧。

但是我用picard去重标记重复后的文件为48.23G,用Sambamba去重后的文件为46.27g,这是否预示着picard标记了更多重复?

运行sambamsa软件,我遇到一个Too many open files问题

sambamba-markdup: Cannot open file `/tmp/sambamba-pid87983-markdup-xrmt/PairedEndsInfozkeh29' in mode `w+' (Too many open files)

这是由于我们的系统限制导致的,解决方法有两个。

(1),临时办法

首先

运行命令

ulimit -n

查看线程限制,默认一般为1024,我们修改为10240.

运行命令

ulimit -n 10240

(2),上面的方法重启后就失效了,如果想重启后生效,需要修改/etc/security/limits.conf文件,root权限。

在文件末尾添加

* soft nofile 10240

* hard nofile 10240

*代表所有该系统的用户

这个10240值不能大于/proc/sys/fs/nr_open文件中的数值,否则将无法正常登录系统。

我的nr_open值为1048576

2,samtools

samtools也有标记或者删除重复序列的功能,命令分别为

samtools rmdup

删除重复序列,但是软件还是推荐标记重复序列

samtools markdup -@ 8 S347.merged.sorted.viewed.bam S347.sorted.viewed.markdupba.bam

通过资料查询,发现这个工具对于bam有一定要求:bam文件以染色体位置排序,MC tag需要用fixmate -m来提供,fixmate -m操作的bam文件必须是利用read name来排序的。

我放弃了,没有测试。

3,samblaster

地址:https://github.com/GregoryFaust/samblaster

命令示例:

samblaster -i samp.disc.sam -o samp.split.sam

没有多线程,操作的是sam文件,没有测试。

4,biobambam2

地址:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4075596/

https://github.com/gt1/biobambam2

bammarkduplicates功能可以去重,没测试

5,bamUtil

https://github.com/statgen/bamUtil

有兴趣的试试吧



生成索引

1,用Sambamba去重时,会自动生成bai索引文件。

但是当我使用命令

sambamba index -p S347.sorted.viewed.markdup.bam

却显示错误

Segmentation fault (core dumped)

2,picard MarkDuplicates去重时加入CREATE_INDEX=true参数,也可以建立索引。

单独使用picard BuildBamIndex命令

picard BuildBamIndex -Xmx64g I=S352.sorted.viewed.markdup.bam

对于picard MarkDuplicates生成的bam提示“非bam”错误

picard.sam.BuildBamIndex done. Elapsed time: 0.00 minutes.

Runtime.totalMemory()=2147483648

To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp

Exception in thread “main” htsjdk.samtools.SAMException: Input file must be bam file, not sam file.

这个真是搞不懂了

用flagstat检查了一下,好像bam文件正常啊

samtools flagstat -@ 30 S346.sorted.viewed.markdup.bam >S346.sorted.markdup.stat

对于sambamba markdup生成的bam,不说不是bam文件了,但是提示“ Insert size out of range”错误

picard.sam.BuildBamIndex done. Elapsed t ime: 29.70 minutes.

Runtime.totalMemory()=5200936960

To get help, see http://broadinstitute.github.io/picard/index.html#Gett ingHelp

Exception in thread “main” htsjdk.samtools.SAMFormatException: SAM vali dation error: ERROR: Record 704691907, Read name A00174:10:HCCJKDSXX:4: 2426:2564:33301, Insert size out of range

at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.ja va:455)

3,使用samtools index建立索引

samtools index S347.sorted.viewed.markdup.bam

却提示错误

Region 536874585..536874735 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6

samtools index: failed to create index for “S347.sorted.viewed.markdup.bam”: Numerical result out of range

查询了一下资料,原始BAM索引格式(BAI)具有固定的bin大小系统,并且最大参考序列长度为512Mbp,使其不适合大染色体。

辣椒的参考基因组cm334,zunla超过3个G,远远大于536,870,912,bai格式索引不适合辣椒基因组分析!

辣椒只能使用csi格式索引文件,然后使用bcftools来call snp了

samtools index -c S347.sorted.viewed.markdup.bam

把变异检测前对数据进行预处理的命令总结一下:

#数据指控

fastp -q 15 -u 40 -M 0 -n 5 -l 80 -w 16 -i  S347_jikangda-A_Green-beifen-1_AHCCJKDSXX_S347_L004_R1_001.fastq.gz -I  S347_jikangda-A_Green-beifen-1_AHCCJKDSXX_S347_L004_R2_001.fastq.gz -o S347.clean_R1.fastq.gz -O S347.clean_R2.fastq.gz -h report/S347-report.html

#建立参考基因组索引

bwa index Capsicum.annuum.L_Zunla-1_Release_2.0.fasta -p zunla

#与参考基因组比对并转换成bam文件

bwa mem -k 32 -M -t 16 -R '@RG\tID:S347\tPL:illumina\tLB:S347\tSM:S347'  zunla  S347.clean_R1.fastq.gz  S347.clean_R2.fastq.gz | samtools view -S -b - > S347.bam

#对bam文件排序

samtools sort -@ 10 -m 80G -O bam -o S347.sorted.viewed.bam S347.bam

#merge两个lanes

samtools merge S347.merged.sorted.viewed.bam S347.sorted.viewed.bam S411.sorted.viewed.bam (此处samtools 其实可以多线程

#去除PCR重复序列

picard -Xmx40g  MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt

#构建bam文件的csi索引

samtools index -c S347.sorted.viewed.markdup.bam

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

推荐阅读更多精彩内容