MOABS: error information: XR:Z or ZR:Z:, or ZS:Z: field not found, please enable -R option in bsmap

[ji zhang]

Dear Professor Sun:

          You have published a software called "MOABS'', which is a very good job I think. Nowadays, I use it to analyze DNA methylation data.

         I successfully installed moabs and it runs very smoothly with the test run. I’m trying to run mcall on a bismark.sam alignment file, but bumped into the following error:

        $ mcall -m ./bismark/ NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam  -r ../../../UCSC/hg38_nochr.fa -a 0 --outputDir ./MOABS/

Options are saved in file run.config and printed here:
cytosineMinScore=20
excludedFlag=0
fullMode=0
keepTemp=0
mappedFiles=./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam
minFragSize=0
minMMFragSize=0
nextBaseMinScore=3
outputDir=./MOABS/
processPEOverlapSeq=1
qualityScoreBase=0
reference=../../../UCSC/hg38_nochr.fa
reportCHX=X
reportCpX=G
requiredFlag=0
skipRandomChrom=1
statsOnly=0
threads=1
trimRRBSEndRepairSeq=2
trimWGBSEndRepairPE1Seq=3
trimWGBSEndRepairPE2Seq=3
Program started
From the extension of file ./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam, program is parsing file according to BAM foramt
XR:Z or ZR:Z:, or ZS:Z: field not found, please enable -R option in bsmap
For file ./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam, the quality score format is Sanger format based at 33!
Protocol and read length are detected as WGBS and 140 bases for file ./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam

For file ./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam the number of all reads is 759969592 and the number of mapped reads is 759969592
Start processing file /share/pub/zhangj/Project/NE2/WGBS/bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam on chrom 1
Segmentation fault

    I checked my bam file, I do have XR:Z field

$samtools view ./bismark/NE2WGBS-trimmed5-pair1_bismark_bt2_pe.sort.bam | head

E00514:540:H22MMCCX2:6:1107:11718:54594_1:N:0:GAATTCGT+GCCTCTAT 83 1 9988 8 97M = 9992 -101 GGAGAGAGGGAGGTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ NM:i:13 MD:Z:0N0N0N0N0N0N0N0N0N0N0N0N0N84 XM:Z:................................................................................................. XR:Z:CT XG:Z:GA
E00514:540:H22MMCCX2:6:2114:26819:52660_1:N:0:GAATTCGT+GCCTCTAT 83 1 9989 8 91M1D25M = 10029 -156 GGAGGATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTACCCTAACCCTAACCCTAACCCTAA FJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJFJFJJJJJJJJJJJFJJJJJJJJFJJJJJJJFJJJJJJJJJJJJ NM:i:13 MD:Z:0N0N0N0N0N0N0N0N0N0N0N0N79^A25 XM:Z:.................................................................................................................... XR:Z:CT XG:Z:GA

Do you have any idea how to fix it?  Any commands or suggestions will be really appreciated.

Best regards,

Ji Zhang

------------------

graduate  student of biomedical engineering,

Wenzhou Medical University, China

[Owen Chen]

Dear Ji Zhang,

For the issue of "XR:Z or ZR:Z:, or ZS:Z: field not found, please enable -R option in bsmap", in addition to the "ZS:Z" fiel
d, you need "XR:Z or ZR:Z" field in the bam file too. 

You can set "XR:Z" field in the bam file by enable -R option in bsmap which is the suggestion moabs gave.      

You need to modify the main program moabs to achieve this. 

Open the master script moabs located in path "${MOABS_ROOT}/bin/moabs" using any editor you like, in the line 570 of the scr
ipt, you should see something like this 'doSys("$mapper -a $mates[0] -b $mates[1] -o $o $others > moabs.bsmap.$o.log 2>&1");
', change it to 'doSys("$mapper -R -a $mates[0] -b $mates[1] -o $o $others > moabs.bsmap.$o.log 2>&1");'.
You need also change line 575 which locate in another if arm from  'doSys("$mapper -a $mates[0]              -o $o $others >
moabs.bsmap.$o.log 2>&1");' to 'doSys("$mapper -R -a $mates[0]              -o $o $others > moabs.bsmap.$o.log 2>&1");'. 

Rerun it the issue should disappear.

br

Owen