hi @ksahlin , your approach looks really promising. I already started to work with it and do some benchmarks in respect with downstream analyses. Would it be possible to allow StrobeAlign to directly print the sam alignments to standard output? that would be very useful to pipe downstream tasks, such as marking duplicates and sorting with samtools. Thanks in advance!

I have illumina paired reads 150bps that is targeted. The average coverage is pretty high. There is a known duplication of 77bps that is homozygous. It looks very good in IGV see snapshot.

With the default parameters the alignment looks like this. Is there way I can adjust the settings to get the alignment to properly anchor both ends of the reads.

hi @ksahlin , sorry for coming back with a new possible issue! I used the very latest patch you've made to align a full sample against the public genome we discussed later. The output was correctly piped to samtools for further duplicate marking and sorting, so thank you very much for allowing the redirection !

Total mapping sites tried: 260811459
Total calls to ssw: 83238353
Calls to ksw (rescue mode): 24136154
Did not fit strobe start site: 452500
Tried rescue: 8059431
Total time mapping: 4317.84 s.
Total time reading read-file(s): 178.746 s.
Total time creating strobemers: 188.722 s.
Total time finding NAMs (non-rescue mode): 618.305 s.
Total time finding NAMs (rescue mode): 341.146 s.
Total time sorting NAMs (candidate sites): 75.2316 s.
Total time reverse compl seq: 0 s.
Total time base level alignment (ssw): 2239.28 s.
Total time writing alignment to files: 525.32 s.

Everything goes well, but looking at the standard output, I obtained the following message several times: ALIGNMENT TO REF LONGER THAN 2000bp - REPORT TO DEVELOPER.

For example: ALIGNMENT TO REF LONGER THAN 2000bp - REPORT TO DEVELOPER. Happened for read: TGCCAATCCTGTGGGACGAATATGGGGCTCAAACCTCAGCCAAAACTCAATAGACACAGTGACGAATGTCTGGTAAAAAATTCAGACCAAAATACCAAAGGAGTAAGGCGTAGCAAGTCCCAGACCGAGAGTGAATAAAACCGGTTTTCCG ref len:53438756

Do you have any idea about such a behavior? I ran strobealign with default parameters.

Hi, using GATK as a variant caller raises an error about the sort order of the chromosomes and indeed they are not sorted although samtools sort was used. The header looks like this

@HD     VN:1.6  SO:coordinate
@SQ     SN:StSOLv1.1ch07_RagTag LN:51859799
@SQ     SN:StSOLv1.1ch01_RagTag LN:89994189
@SQ     SN:StSOLv1.1ch12_RagTag LN:108042198
@SQ     SN:Chr0_RagTag  LN:74722543
@SQ     SN:StSOLv1.1ch08_RagTag LN:122763581
@SQ     SN:StSOLv1.1ch02_RagTag LN:52958123
@SQ     SN:StSOLv1.1ch11_RagTag LN:82684344
@SQ     SN:StSOLv1.1ch09_RagTag LN:110327934
@SQ     SN:StSOLv1.1ch06_RagTag LN:77416214
@SQ     SN:StSOLv1.1ch04_RagTag LN:108266590
@SQ     SN:StSOLv1.1ch05_RagTag LN:91627361
@SQ     SN:StSOLv1.1ch10_RagTag LN:63742240
@SQ     SN:StSOLv1.1ch03_RagTag LN:45113244
@PG     ID:strobealign  PN:strobealign  VN:0.7  CL:strobealign
@PG     ID:samtools     PN:samtools     PP:strobealign  VN:1.13 CL:samtools fixmate [email protected] -u -m - -
@PG     ID:samtools.1   PN:samtools     PP:samtools     VN:1.13 CL:samtools view -bhS -
@PG     ID:samtools.2   PN:samtools     PP:samtools.1   VN:1.13 CL:samtools sort [email protected]
@PG     ID:samtools.3   PN:samtools     PP:samtools.2   VN:1.13 CL:samtools view -H Altus.ColombaNRGene.ragtag.strobealign.bam

any ideas what is going wrong ? Other aligners like bwa produce a proper sorted bam file

Dear StrobeAlign author,

Use default settings, the mapping sam file output cannot be transformed into sorted bam file by samtools:

(base) [[email protected] Competitive_mapping]$ time StrobeAlign -t 24 T4Aer_MAG.fasta T4AerOil_R1.fa T4AerOil_R2.fa Using n: 2 k: 22 s: 18 t: 24 R: 2 w_min: 6 w_max: 14 [w_min, w_max] under thinning w roughly corresponds to sampling from downstream read coordinates (under random minimizer sampling): [30, 70] Time reading references: 7.76452 s

ref vector approximate size: 9290379 Ref vector actual size: 8664874 Unique strobemers: 8651483 Total time generating flat vector: 2.85387 s

Flat vector size: 8664874 Total strobemers count: 8664874 Total strobemers occur once: 8639067 Total strobemers highly abundant > 100: 0 Total strobemers mid abundance (between 2-100): 12415 Total distinct strobemers stored: 8651483 Ratio distinct to non distinct: 696 Filtered cutoff index: 1730 Filtered cutoff count: 2

Total time generating hash table index: 0.72899 s

Total time indexing: 11.3475 s

Using rescue cutoff: 4 Running PE mode Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 104.07, stddev: 124.369) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 96.6191, stddev: 175.71) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 90.6641, stddev: 208.097) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 117.778, stddev: 245.558) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 106.226, stddev: 274.345) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 101.158, stddev: 301.601) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 97.5571, stddev: 323.28) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 102.831, stddev: 345.341) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 109.84, stddev: 371.438) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 114.511, stddev: 392.006) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 117.354, stddev: 407.773) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 138.058, stddev: 426.87) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 114.614, stddev: 440.569) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 108.915, stddev: 452.172) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 132.276, stddev: 467.806) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 129.727, stddev: 483.589) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 109.746, stddev: 500.041) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 107.748, stddev: 512.066) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 107.68, stddev: 525.201) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 117.557, stddev: 539.586) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 107.239, stddev: 553.863) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 116.547, stddev: 566.07) Mapping chunk of 1000000 query sequences... Estimated diff in start coordinates b/t mates, (mean: 110.742, stddev: 578.011) Mapping chunk of 148502 query sequences... Estimated diff in start coordinates b/t mates, (mean: 110.161, stddev: 589.893) Total mapping sites tried: 13055987 Total calls to ksw: 3218137 Calls to ksw (rescue mode): 1276608 Did not fit strobe start site: 1776340 Tried rescue: 33208162 Total time mapping: 103.476 s. Total time reading read-file(s): 25.1783 s. Total time creating strobemers: 18.6664 s. Total time finding NAMs (non-rescue mode): 3.60437 s. Total time finding NAMs (rescue mode): 0.409914 s. Total time sorting NAMs (candidate sites): 0.0890606 s. Total time reverse compl seq: 0 s. Total time extending alignment: 17.081 s. Total time writing alignment to files: 5.61499 s.

real 1m55.191s user 17m49.681s sys 0m14.051s (base) [[email protected] Competitive_mapping]$ ls -lhs mapped.sam 7.3G -rw-r--r-- 1 jzhao399 p-ktk3 7.3G Nov 21 20:23 mapped.sam (base) [[email protected] Competitive_mapping]$ samtools view -bS [email protected] 24 mapped.sam | samtools sort [email protected] 24 -O bam -o mapped_sorted.bam samtools view: error reading file "mapped.sam": Input/output error samtools view: error closing "mapped.sam": -5 [bam_sort_core] merging from 0 files and 24 in-memory blocks...

Any idea?

Thanks,

Jianshu

Hi @ksahlin , I am giving a try to strobealign v0.6. I was able to compile it correctly on my computing nodes. Genome indexing works fine but in the PE mode, it crashes from the beginning of the mapping part:

...
Using rescue cutoff: 208
@SQ     SN:Chr01        LN:51433939
@SQ     SN:Chr04        LN:48048378
@SQ     SN:Chr08        LN:63048260
@SQ     SN:Chr11        LN:53580169
@SQ     SN:Chr09        LN:38250102
@SQ     SN:Chr05        LN:40923498
@SQ     SN:Chr02        LN:49670989
@SQ     SN:Chr10        LN:44302882
@SQ     SN:Chr07        LN:40041001
@SQ     SN:Chr06        LN:31236378
@SQ     SN:Chr03        LN:53438756
@PG     ID:strobealign  PN:strobealign  VN:0.6  CL:strobealign
Running PE mode
Mapping chunk of 1000000 query sequences...
terminate called after throwing an instance of 'std::out_of_range'
  what():  basic_string::substr
Aborted (core dumped)

This is strange, because running it in SE mode with either the forward or the reverse reads work perfectly. Any idea about that possible issue?

Thanks!

Hello,

first of all thanks for this efficient software. I was wondering if there is currently any option to retain all the reads with asociated flags on .sam output file, especially the FLAG tag specifying if the read is mapped 0, 1 or >1 time. This behavior is the same like -a option of Bowtie2.

Multimapping reads seems to be discadreded at the time am I right ?

Thanks in advance

From issue #3:

"Small suggestion: remove the fasta annotation (only keep after > and before the first tab) from reference in sam output so that the sam output can be a little bit smaller. Or add an option to only keep them or an option to include annotation."

Hi there,

tried to compile v0.4 (on Ubuntu) from source via cmake, but - https://github.com/ksahlin/StrobeAlign/blob/main/CMakeLists.txt#L6 refers to source/edlib.cpp source/edlib.h which are not present in the repository/release tarball.

Compiling via g++ -std=c++14 main.cpp source/index.cpp source/xxhash.c source/ksw2_extz2_sse.c source/ssw_cpp.cpp source/ssw.c -lz -fopenmp -o strobealign -O3 -mavx2, as outlined in the README, works as intended, though.

Also, could you consider lowering cmake_minimum_required? 3.19 isn't widely deployed, Ubuntu 20.04 LTS currently ships with 3.16.3, for example.

Thanks a lot!

Hi @ksahlin,

for automated pipelines that employ strobealign, it's difficult to come up with a sensible setting for the -r parameter (even though it has a default value) since datasets to be processed will likely have different read lengths. Would it make sense to e.g. automatically infer a value here by looking at the first n reads to be mapped?

Hi @ksahlin,

I've now also added an install target to cmake; by default, make install would attempt installation to /usr/local, but this can be adapted with -DCMAKE_INSTALL_PREFIX:PATH=${PREFIX} when invoking cmake.

It would be great if StrobeAlign could (optionally) add an @RG header and RG tags to each read. Having RG tags is a requirement in some of our pipelines and for some programs, and I would consider it best practice. It’s possible, but awkward to add read group tags afterwards with samtools addreplacrg.

For compatibility with BWA(-MEM)/minimap2, I would suggest to use command-line parameter -R, but at the moment, that is already used to set the "Rescue level". And although there would certainly be nicer ways to provide read group information than the way in which BWA does it, it would be good to at least also support the syntax that it uses (-R '@RG\tID:foo\tSM:bar').

The program exits with a segfault when the reference FASTA file does not exist:

$ strobealign ref.fasta dummy.fasta
...
Unique strobemers: 1
Total time generating flat vector: 3.2896e-05 s

Flat vector size: 0
Segmentation fault (core dumped)

The same happens when a compressed FASTA (.gz) was provided (even if it exists).

Note to developer:

Using 32-bit hashes may reduce both vector and hash table size. For human, it is expected to reduce peak memory from 32Gb (64bit) down to around 20-24Gb

At the final hash value computation h1/2+h2/2 when the seed has been decided, simply hash this 64 bit down to 32 bit with wyhash or XXH32.

This will result in more collisions on large genomes though as there are only 4.2billion unique slots. For human, we store about 450M unique values, but for larger genomes this may be a problem (can this be solved with 32/64 bit template?).

I see quite a reduction in both space and alignment time if we change to 32bit hash value (on a small dataset!).

Note to developer:

The extension step (nucleotide level alignment) is the bottleneck in strobealign. There are different three ways to reduce this:

1. Direction 1 (change the alignment module):
   1. Change to base level alignment with WFA (WFA publ) as is done in Accelalign
2. Direction 2 Speedup the current module used (SSW):
   1. By using 8bit slots in alignment matrix?
   2. By not computing alignment twice - is ssw does this?
3. Direction 3 (partitioned SSW)
   1. Finish implementing partitioned SW (split alignment into several small hamming or SW alignments) if seeds are in middle.

Note to developer:

It is not clear if strobealign has optimal mapq scores for variant calling. For example, having an accurate MAPQ score is crucial for SNP calling. See this tweet thread and this subthread.

For info, a relevant paper about this is found here.

Note to developer:

1. Use Randstobes order n=3.
2. Collinear chaining: is needed, perhaps heuristic by sorting based on reference (see find_NAMs_alt() function)
3. How to implement split mapping?
4. Need partitioned local alignments of reads (regions extracted from seeds) - not SW-aligning the whole read at once