align Falcon-unzip contigs to itself using minimap

Question

[olechnwin]

Hello All,

I've run Falcon unzip and got the primary contigs and haplotigs of my denovo assembled genome.

I am trying to visualize the alignment between the primary contigs (cns_p_ctg.fasta) and haplotigs (cns_h_ctg.fasta). I tried to concatenate primary contigs and haplotigs as the reference genome (cns_p_h_ctg.fasta) then align each primary contigs and haplotigs to this reference genome using minimap.

Basically this is what I tried to do:

concat cns_p_ctg.fasta cns_h_ctg.fasta > cns_p_h_ctg.fasta
minimap2 -d cns_p_h_ctg.mmi cns_p_h_ctg.fasta
minimap2 -a cns_p_h_ctg.mmi cns_h_ctg.fasta > cns_h_ctg_aligned_cns_p_h_ctg.sam
samtools view -bS cns_h_ctg_aligned_cns_p_h_ctg.sam > cns_h_ctg_aligned_cns_p_h_ctg.aligned.bam

However, at the end no bam or sam files were generated. These are the log:

+ /opt/minimap2-2.12_x64-linux/minimap2 -d cns_p_h_ctg.mmi cns_p_h_ctg.fasta
[M::mm_idx_gen::86.704*1.79] collected minimizers
[M::mm_idx_gen::105.810*2.01] sorted minimizers
[M::main::114.081*1.93] loaded/built the index for 13543 target sequence(s)
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 13543
[M::mm_idx_stat::114.990*1.92] distinct minimizers: 99850267 (23.08% are singletons); average occurrences: 7.594; average spacing: 5.308
[M::mm_idx_gen::123.520*1.89] collected minimizers
[M::mm_idx_gen::125.275*1.91] sorted minimizers
[M::main::126.829*1.90] loaded/built the index for 3528 target sequence(s)
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 3528
[M::mm_idx_stat::127.356*1.89] distinct minimizers: 29020765 (67.65% are singletons); average occurrences: 2.047; average spacing: 5.310
[M::main] Version: 2.12-r827
[M::main] CMD: /opt/minimap2-2.12_x64-linux/minimap2 -d cns_p_h_ctg.mmi cns_p_h_ctg.fasta
[M::main] Real time: 127.812 sec; CPU: 241.706 sec
+ /opt/minimap2-2.12_x64-linux/minimap2 -a cns_p_h_ctg.mmi cns_h_ctg.fasta
[WARNING] For a multi-part index, no @SQ lines will be outputted. Please use --split-prefix.
[M::main::9.474*0.99] loaded/built the index for 13543 target sequence(s)
[M::mm_mapopt_update::10.820*0.99] mid_occ = 842
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 13543
[M::mm_idx_stat::11.655*0.99] distinct minimizers: 99850267 (23.08% are singletons); average occurrences: 7.594; average spacing: 5.308
[M::worker_pipeline::285.765*2.78] mapped 3652 sequences
[M::worker_pipeline::547.295*2.88] mapped 3964 sequences
[M::worker_pipeline::991.046*2.92] mapped 4739 sequences
[M::worker_pipeline::1849.678*2.91] mapped 1121 sequences
[M::main::1852.085*2.90] loaded/built the index for 3528 target sequence(s)
[M::mm_mapopt_update::1852.085*2.90] mid_occ = 842
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 3528
[M::mm_idx_stat::1852.564*2.90] distinct minimizers: 29020765 (67.65% are singletons); average occurrences: 2.047; average spacing: 5.310
[M::worker_pipeline::2310.411*2.92] mapped 3652 sequences
[M::worker_pipeline::2753.968*2.94] mapped 3964 sequences
[M::worker_pipeline::3307.200*2.94] mapped 4739 sequences
[M::worker_pipeline::4549.630*2.90] mapped 1121 sequences
[M::main] Version: 2.12-r827
[M::main] CMD: /opt/minimap2-2.12_x64-linux/minimap2 -a cns_p_h_ctg.mmi cns_h_ctg.fasta
[M::main] Real time: 4549.653 sec; CPU: 13212.292 sec
+ samtools view -bS cns_h_ctg_aligned_cns_p_h_ctg.sam
[E::sam_parse1] missing SAM header
[W::sam_read1] Parse error at line 2
[main_samview] truncated file.

I cannot figure out how to align each individual contigs to a concatenated contigs as "reference genome". Tried split-prefix but was not sure what to input.

Thank you in advance for your help.

UPDATE:

Thanks to the link provided below by ATpoint.

The problem is solved by manually adding sam header. Basically just change the samtools line above to this one (adding -T reference.fasta) :

samtools view -b -T cns_p_h_ctg.fasta cns_h_ctg_aligned_cns_p_h_ctg.sam > cns_h_ctg_aligned_cns_p_h_ctg.aligned.bam

Answer 1

[ATpoint]

Maybe see https://github.com/lh3/minimap2/issues/171

[olechnwin]

Thank you so much for referring me to that issue. That fixed my problem. I wonder how to accept your comment as answer.

[ATpoint]

Glad it worked.

[olechnwin]

This is unrelated but I remember yesterday there was only a button to upvote. There was no button to accept. What did you do?

[ATpoint]

I moved my comment to answer. You can only accept answers, no accept buttom for comments ;-)