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

[Xiaoxiao Hao]

Hi Daqiang,

I want to remove duplications before doing the mcall.

The bam I got from your another software LiBis, the command I used as follows:

LiBis -p 3 -n ${DIR_fastq}/IR_Csc1-r1_val_1.fq,${DIR_fastq}/IR_Csc1-r2_val_2.fq -r ${DIR_reference}/Homo_sapiens.GRCh37.73.dna.primary_assembly.fa 

1) remove duplications as follows:

We split bam files base on the tag using Bamtools, then use Picard removed duplications, at last merged all the bam files to one bam file using Picard MergeSamFiles.

2) I tried to use your MOABS mcall  get the CpG sites.

It always give me the following error information;

From the extension of file ./IR_Csc1-r1_val_1.markdup.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 ./IR_Csc1-r1_val_1.markdup.bam, the quality score format is Sanger format based at 33!

Protocol and read length are detected as WGBS and 144 bases for file ./IR_Csc1-r1_val_1.markdup.bam

For file ./IR_Csc1-r1_val_1.markdup.bam the number of all reads is 11991446 and the number of mapped reads is 11991446

Start processing file ./IR_Csc1-r1_val_1.markdup.bam on chrom 1

terminate called after throwing an instance of 'std::out_of_range'

  what():  basic_string::substr

/data/eg00/sgea/sge/ge2011.11/cell0/spool/EG01/job_scripts/193221: line 18:  6890 Aborted                 (core dumped) ../../moabs-v1.3.2.src.x86_64_Linux.data/bin/mcall -m ./IR_Csc1-r1_val_1.markdup.bam -r ${DIR_reference}/Homo_sapiens.GRCh37.73.dna.primary_assembly.fa

3) I checked my bam file, I do have ZS:Z field

[xhao@loginnode2 mcall-test-new]$ samtools view ../mcall-test/IR_Csc1-r1_val_1.markdup.bam | head -n 20 

E00526:319:H7Y7GCCX2:7:1208:19705:59710 73 1 10001 255 142M * 0 0 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAACCTACCCCAACCCTAACCCTAACCCTAACCCTAACCCTAACC FJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJF<AFJJJJAFJFJFFFJJ<<7FFJAA<AFFJJ7FJA7<--<7--7FJ--7-7A---7FF--77AAFF-AJJF-7FA7A---F PG:Z:MarkDuplicates NM:i:2 ZS:Z:++

E00526:319:H7Y7GCCX2:7:1123:7436:50656 147 1 13449 255 100M = 13451 -98 CCAATTTCACCAAAAATAAACCTCTTCCTAATAAACAACTACACCACTACCTAACACTATACCCTTCCTTTACTCTACCCACTAAAAACAATATTTATCA JJJ<FFFJJJJFFJJJFJFJJFJJFJJFJJFJJJJAFAFFFAAFFJFFJJFAJFJJJJJJJJJJJJJJJJFFFFFF7JJJFFJJAJJJJJJJJJJAFJJF PG:Z:MarkDuplicates.4 NM:i:1 ZS:Z:-+

E00526:319:H7Y7GCCX2:7:1123:7436:50656 99 1 13451 255 100M = 13449 98 AATTTCACCAAAAATAAACCTCTTCCTAATAAACAACTACACCACTACCTAACACTATACCCTTCCTTTACTCTACCCACTAAAAACAATATTTATCATA FJJAFAFFFJJAAAJF-FJAA<<<AFFFJFJJFJJAJF<F-FFFF7-<JJF<JJJJJF-FJJJJJJ7FJ-FFFF7JAFFJJJJJFJJJFJJ<JFAJJJFF PG:Z:MarkDuplicates.5 NM:i:1 ZS:Z:--

Any commands or suggestions will be really appreciated. 

Best,

Xiaoxiao

[ji zhang]

Hi Xiaoxiao,

I ran into the same problem. Have you found a solution? can you share?

Thanks!

[Owen Chen]

Dear xiaoxiao and 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.

For the issue of "throwing an instance of 'std::out_of_range'", it normally caused by having blank characters in the fasta d
efinition line of the reference file.
Check your reference file which is "Homo_sapiens.GRCh37.73.dna.primary_assembly.fa", if blanks do exist, remove it or change
it to "_". For example, change from ">chr1  AC:CM000663.2  gi:568336023" to ">chr1_AC:CM000663.2_gi:568336023" or ">chr1AC:
CM000663.2gi:568336023".
This issue actually already answered in the discussion group, you can find it in the link below.
"https://groups.google.com/g/moabs_msuite/c/UhzSqVV1a6M/m/sPQRNaFoBgAJ"
It may save you some time if you try search for the issue first.

br
Owen