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!

… call a helper that returns a reference to a vetor of alignments and align_PE simply returns the best (by SW score as it did before) while align_SE_secondary_hits does it's usual thing with the vector.

I am relatively certain that this will do the same thing as before, however since we now have unit tests, adding some unit tests that check that alignments look as they should with test data could be really good.

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.

I did two things:

1. Collected the index into a struct.
2. Added the possibility to read/write the index from/to disk.

I have not created test cases (I don't know how in this project) but have tested the code, it produces identical output for phix compared to previous runs.

I tested to run on the human transcriptome: Read of reference: 4.3 s Generation of index: 31.1 s Writing index: 6.8 s - this is a bit slow, it is probably possible to speed it up but I left it for now. Reading index: 3.6 s Index size: 1.5 GB - the kallisto index for the same is 2.4 GB, somewhat the same order of magnitude, the fasta file is 350 MB.

I developed the code in Visual Studio, worked fine. I have not tested to compile with GCC, but it should work. Would be good if someone with that environment up could test it quickly.

There are three things we discovered that should be changed later at some point, but it is better to do these as separate actions:

• idx_to_acc - change this to a vector, no point in having a map with <i,data>, where i is just the index - access will be identical, i.e., index.acc_map[i]
• Remove the first uint64_t in the mers_vector - according to Kristoffer it should not be needed, which will reduce the memory taken up by the index (and file size).
• I think we should rename unsigned int etc. to uint32_t - the number of bits in types like int has traditionally been different between machine architectures to optimize speed, so the size is a bit unclear, although it is pretty safe to assume they will be 32 bits.

I can fix these things once this is approved.

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

This is to document which accuracy changes I found since v0.7.1.

Test data

Only paired-end reads.

• simulated Drosophila
  • Queries: First 1 million 50 bp reads from /proj/snic2022-6-31/nobackup/strobemap_eval/reads_PE/drosophila/
  • Reference: /proj/snic2022-6-31/nobackup/strobemap_eval/genomes/drosophila.fa
• real Drosophila
  • Queries: First 1 million reads of SRR6055476
  • Reference: GCF_000001215.4
• CHM13
  • Queries: First 100'000 reads of /proj/snic2022-6-31/nobackup/strobemap_eval/reads_PE/CHM13/200_?.fq
  • Reference: /proj/snic2022-6-31/nobackup/strobemap_eval/genomes/CHM13.fa

Command used: strobealign-${commithash} -t 1 ${reference} ${r1} ${r2} | samtools view --no-PG -o out.bam

Commits with accuracy or other output changes

|Commit|Accuracy in CHM13|Description|Fix/discussion|Trigger in real Drosophila| |-|-|-|-|-| |v0.7.1|94.3725||baseline|| |cd46f8f|94.363|Use std::{max,min,abs} instead of ternary operators|#136|SRR6055476.524| |14d512f|94.3650|Fix a probable logic error due to a typo||| |d987274|94.169|Do not convert reference to uppercase in seq_to_randstrobes2|#138|?| |13bdc5d|94.169|Factor out rescue_read|#139|SRR6055476.15| |6973009|92.875|Invert conditions to allow early break|#127|?| |d82ce31|93.997|Bug in only this one commit: filter_cutoff was set incorrectly||| |720887d|92.876|Fix incorrect refactor, part 1||| |8d512a9|94.169|Fix incorrect refactor, part 2||| |e758839|94.1765|Merge pull request #136 from ksahlin/fix-max|| |464a852|94.3745|Merge pull request https://github.com/ksahlin/StrobeAlign/pull/138 from ksahlin/uppercase||

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

strobealign StaphylococcusLT963435.fasta fastpSRR13329724_1.fastq.gz fastpSRR13329724_2.fastq.gz > outputStaphreads.sam Using k: 20 s: 16 w_min: 5 w_max: 11 Read length (r): 129 Maximum seed length: 79 Threads: 3 R: 2 Expected [w_min, w_max] in #syncmers: [5, 11] Expected [w_min, w_max] in #nucleotides: [25, 55] A: 2 B: 8 O: 12 E: 1 Time reading references: 0.01117 s

ref vector approximate size: 555450 Illegal instruction (strobealign) [email protected]:~/Turtle$

To solve #31 and #56, we should create an AlignmentRecord or so class (or something like that). The idea is that objects of this class represent the information from one row in a SAM/BAM file, and that there would be a method that writes out the data in SAM format. All the functions that want to write SAM output then would create one instance per row that they want to write and let the class do the formatting. That way, the formatting (serialization) is done in one place only and could, at least in theory, be much more easily switched for something else, for example, if we wanted to support writing BAM or CRAM (not that I suggest this, I think writing SAM is totally fine).

Possibly some kind of SamFile or so class could be involved in order to store the read group ID that we want to add as part of #31.

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!

Picard ValidateSamFile has many complaints like this:

ERROR::INVALID_ALIGNMENT_START:Record 23701, Read name SRR6055476.11851, Mate Alignment start (3496047) must be <= reference sequence length (44411) on reference NW_007931074.1

Commit 5b462d6 introduced a change in alignment results. This is the diff of the alignments (- is before, + is after) for the first of two affected read pairs in the simulated Drosophila dataset:

-simulated.262476/1 97  211000022278483 515 60  2S48=   3L  14468499    14468024    TGGAAATTAAATTTTTTGGCATTATTTTGCAAATTTTGATGACCCCCCTC  HIHHGGIHHDHHFFIGIHIDIIHIEFG>IGEIFDIEHIBIG;I7>[email protected]    NM:i:0  AS:i:96
+simulated.262476/1 65  211000022278483 515 60  2S48=   211000022278398 443 -122    TGGAAATTAAATTTTTTGGCATTATTTTGCAAATTTTGATGACCCCCCTC  HIHHGGIHHDHHFFIGIHIDIIHIEFG>IGEIFDIEHIBIG;I7>[email protected]    NM:i:0  AS:i:96

-simulated.262476/2 145 3L              14468499  60  9S10=1X30=  211000022278483 515 -14468024   ATATGGAATGTGTTGGAATNTACTATTCAACCTACAAAAATAACGTTAAA  HIII7CIIIFIIIEIFIIGIHIB8IIIIIIBIIFIIFIEGIGHIIIGHIH    NM:i:1  AS:i:72
+simulated.262476/2 129 211000022278398      443  60  30=1X9=10S  211000022278483 515       122   TTTAACGTTATTTTTGTAGGTTGAATAGTANATTCCAACACATTCCATATthe new version chooses an alignment that is worse than the best (given how the 'best' is defined) - as none of the alignment configurations are 'mapped in proper pair'. (of course the second one is better if inversion is more likely than translocation).

This could be an issue that comes back and haunts us.  HIHGIIIHGIGEIFIIFIIBIIIIII8BIHIGIIFIEIIIFIIIC7IIIH    NM:i:1  AS:i:70

I have no idea at the moment why such a change occurs. It appears that the robin_hood::unordered_map behaves differently when it is cleared and then re-filled compared to creating a new one from scratch.

Originally posted by @marcelm in https://github.com/ksahlin/StrobeAlign/issues/140#issuecomment-1302126477

And @ksahlin replied:

[...] the new version chooses an alignment that is worse than the best (given how the 'best' is defined) - as none of the alignment configurations are 'mapped in proper pair'. (of course the second one is better if inversion is more likely than translocation).

This could be an issue that comes back and haunts us.

Hi, I'm using strobealing to test alignment performance and I get this error from samtools view command with sam created by strobealign (see below command and output) The .fna reference file is Homo sapiens chromosome 18, GRCh38.p14 the sam creation was succesful, aligning with a fastq that is a pbsim simulated directly from the same NC000018.fna What am I doing wrong ? Thanks in advance .

The command samtools view was: samtools view -bS -T NC000018.fna strobealign_ps01_NC000018.sam -O BAM -o strobealign_ps01_NC000018.bam [email protected] 12

and the output was

samtools view: error reading file "NC000018.sam"
[E::read_ncigar] No CIGAR operations
samtools view: error closing "NC000018.sam": -5

This is the command and output of sam ceration with strobealign

 /packages/StrobeAlign/build/strobealign NC000018.fna ps01_NC000018.fastq -o strobealign_ps01_NC000018.sam -t 12

Using
k: 23
s: 17
w_min: 5
w_max: 15
Read length (r): 10032
Maximum seed length: 278
Threads: 12
R: 2
Expected [w_min, w_max] in #syncmers: [5, 15]
Expected [w_min, w_max] in #nucleotides: [35, 105]
A: 2
B: 8
O: 12
E: 1
Time reading references: 0.965166 s

ref vector approximate size: 11481897
Ref vector actual size: 10783788
Time generating seeds: 2.99098 s

Time sorting seeds: 1.18937 s

Unique strobemers: 9922754
Total time generating flat vector: 4.18048 s

Flat vector size: 10783788
Total strobemers count: 10783788
Total strobemers occur once: 9840601
Fraction Unique: 0.991721
Total strobemers highly abundant > 100: 783
Total strobemers mid abundance (between 2-100): 81369
Total distinct strobemers stored: 9922754
Ratio distinct to highly abundant: 12672
Ratio distinct to non distinct: 120
Filtered cutoff index: 1984
Filtered cutoff count: 41


Total time generating hash table index: 0.584856 s

Total time indexing: 5.73055 s

Using rescue cutoff: 82
Running SE mode
Done!
Total mapping sites tried: 502929
Total calls to ssw: 502929
Calls to ksw (rescue mode): 0
Did not fit strobe start site: 1585
Tried rescue: 16813
Total time mapping: 91.5354 s.
Total time reading read-file(s): 13.2334 s.
Total time creating strobemers: 2.68648 s.
Total time finding NAMs (non-rescue mode): 0.603002 s.
Total time finding NAMs (rescue mode): 0.0492406 s.
Total time sorting NAMs (candidate sites): 0.0038492 s.
Total time reverse compl seq: 5.79428e-311 s.
Total time base level alignment (ssw): 3.01523 s.
Total time writing alignment to files: 6.92938e-310 s.

SAM spec:

TLEN: [...]. It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).

A FASTQ record with header @readname/1 or @readname/2 should be named readname in the SAM output.

The SAM spec says:

QNAME: Query template NAME. Reads/segments having identical QNAME are regarded to come from the same template. That is, if we want paired-end reads to be considered coming from the same template, then they need to get the same QNAME.

BWA-MEM does this.