I am trying to build a mappability mask with Heng Li's SNPable program https://lh3lh3.users.sourceforge.net/snpable.shtml . The instructions there are pretty brief and do not detail all the necessary steps, which makes it challenging for introductory bioinformaticians. Suppose the reference genome is genome.fa, copy-pasting the given instructions, they are:

Extract all overlapping k-mer subsequences as read sequences:

splitfa genome.fa 35 | split -l 20000000  

Align all reads to the genome with BWA

bwa aln -R 1000000 -O 3 -E 3 genome.index xxaa > xxaa.sai

Suppose the unsorted BWA alignment results are xx??.sam.gz. Generate rawMask with:

gzip -dc xx??.sam.gz | gen_raw_mask.pl > rawMask_35.fa  

Generate the final mask:

gen_mask -l 35 -r 0.5 rawMask_35.fa > mask_35_50.fa   

In step 2, genome.index comes from bwa index. I believe there are a few undetailed steps to do between steps 2 and 3. When I did step 1, it wrote a series of files xaa xab xac ... xlk xlm, then step 2 aligned these and wrote the files xaa.sai xab.sai xac.sai ... xlk.sai xlm.sai. I used bwa samse to generate sam files. For each sai file do:

bwa samse genome.index ${prefix}.sai > ${prefix}.sam

I then merged all of these sam files into one big sam file:

samtools merge merged.sam ${allsamfiles}

and then processed steps 3 and 4 with that merged sam file:

gzip -dc merged.sam | gen_raw_mask.pl > rawMask_35.fa  
gen_mask -l 35 -r 0.5 rawMask_35.fa > mask_35_50.fa   

However, the problems are that this final mask_35_50.fa is not contiguous, in other words there are lots of different sequences for each chromosome, e.g.

cat  mask_35_50.fa | grep ">chr4" | grep -v random | grep -v Un | wc -l # remove scaffolds and contigs with grep -v random and grep -v Un

returns 166, not 1 as I would have hoped.

What is going on here? It seems like I need to sort the merged.sam file, but i tried this in various ways and things got even worse. Using

samtools sort -@ 4 -m 2G merged.sam -o merged_sorted.sam 

then doing steps 3 and 4 with merged_sorted.sam , the final mask was even more fragmented. I also tried sorting by readname

samtools sort -@ 4 -m 2G -n merged.sam -o merged_sortedreadname.sam 

but that also did not work.

Does anybody know what I am doing wrong?