I have a very big reads data from pacbio seq2, and I want to get overlap for all reads. So I split the huge read fasta file, and use minimap2 ava-pb to map each pair of sub files. But I found when I use A as query and B as target will get nothing , but use B as query and A as target will get a good results. Why ??
my code:
minimap2 -t 56 -x ava-pb 031.fa 000.fa > 031_vs_000.paf
this cmd get a good results file
[M::mm_idx_gen::73.558*1.94] collected minimizers
[M::mm_idx_gen::81.744*4.10] sorted minimizers
[M::main::81.745*4.10] loaded/built the index for 250000 target sequence(s)
[M::mm_mapopt_update::90.531*3.80] mid_occ = 357
[M::mm_idx_stat] kmer size: 19; skip: 5; is_hpc: 1; #seq: 250000
[M::mm_idx_stat::95.469*3.65] distinct minimizers: 242745064 (44.26% are singletons); average occurrences: 3.701; average spacing: 4.351
[M::worker_pipeline::153.218*15.64] mapped 32255 sequences
[M::worker_pipeline::187.490*18.34] mapped 31983 sequences
[M::worker_pipeline::220.958*20.28] mapped 33009 sequences
[M::worker_pipeline::252.800*21.56] mapped 33889 sequences
[M::worker_pipeline::286.707*22.49] mapped 33790 sequences
[M::worker_pipeline::319.235*23.42] mapped 34288 sequences
[M::worker_pipeline::345.419*22.96] mapped 35404 sequences
[M::worker_pipeline::355.369*22.34] mapped 15382 sequences
[M::main] Version: 2.17-r941
[M::main] CMD: minimap2 -t 56 -x ava-pb 031.fa 000.fa
[M::main] Real time: 357.013 sec; CPU: 7941.179 sec; Peak RSS: 21.388 GB
minimap2 -t 56 -x ava-pb 000.fa 031.fa > 000_vs_031.paf
this cmd get a empty output paf
[M::mm_idx_gen::70.207*1.95] collected minimizers
[M::mm_idx_gen::78.683*4.33] sorted minimizers
[M::main::78.683*4.33] loaded/built the index for 250000 target sequence(s)
[M::mm_mapopt_update::88.403*3.96] mid_occ = 372
[M::mm_idx_stat] kmer size: 19; skip: 5; is_hpc: 1; #seq: 250000
[M::mm_idx_stat::93.258*3.81] distinct minimizers: 233769422 (45.47% are singletons); average occurrences: 3.633; average spacing: 4.374
[M::worker_pipeline::101.071*5.86] mapped 31360 sequences
[M::worker_pipeline::105.144*7.67] mapped 31463 sequences
[M::worker_pipeline::109.309*9.48] mapped 32001 sequences
[M::worker_pipeline::113.724*10.76] mapped 32064 sequences
[M::worker_pipeline::118.359*12.07] mapped 31908 sequences
[M::worker_pipeline::122.208*13.44] mapped 32472 sequences
[M::worker_pipeline::126.029*14.71] mapped 32316 sequences
[M::worker_pipeline::128.879*15.59] mapped 26416 sequences
[M::main] Version: 2.17-r941
[M::main] CMD: minimap2 -t 56 -x ava-pb 000.fa 031.fa
[M::main] Real time: 130.109 sec; CPU: 2010.615 sec; Peak RSS: 17.669 GB
my subfasta file is very normal
000.fa
Nx Size Number
N90 8491 152518
N80 13211 117916
N70 17485 93628
N60 20827 74277
N50 23564 57542
N40 26399 42639
N30 29694 29362
N20 34052 17646
N10 40813 7596
Total Length 3714217773 250000
Maximum Length 110245
Minimum Length 50
031.fa
Nx Size Number
N90 9065 153364
N80 14104 119225
N70 18485 95239
N60 21673 75819
N50 24399 58842
N40 27290 43686
N30 30747 30171
N20 34974 18222
N10 41607 7899
Total Length 3908445356 250000
Maximum Length 132603
Minimum Length 50
You should set --dual=yes, if query fasta file is different from target fasta file.
In default parameters of -x ava-pb, --dual is set to no to reduce half of alignment work.
@wzboy1984 Thank you very much, I will try it : )
As @wzboy1984 said (Thank you!).
Most helpful comment
You should set --dual=yes, if query fasta file is different from target fasta file.
In default parameters of -x ava-pb, --dual is set to no to reduce half of alignment work.