what is a chimeric block

[sjanssen2]

I find in file /vol/jlab/tlin/all_project/nf_results/refbasedassemblereval/oak/identifychimericblocks/oak_reference_chimeric_blocks.tsv the line

chimblock_50670 GCF_900416725.1 ['NVRI1_RS03545', 'NVRI1_RS03545', 'NVRI1_RS03545', 'NVRI1_RS03545'] [0, 0, 1, 1] NZ_LS450958.2 791346.0 791760.0

which I don't understand. Does it state that chimblock_50670 in genome GCF_900416725.1 is detected, because there is an overlap of the gene NVRI1_RS03545 - but no others? If so, why is it a chimera?! Does not make sense to me?!

[tensulin]

tricky case, NVRI1_RS03545 can apparently be recovered as two distinct blocks, however on the genome the gene is duplicated and is located next to its copy, I also see that the entry for this gene for the same genomic entry is duplicated, I will check whether it is may be due to a missing deduplication step

[sjanssen2]

another edge? case:
gene IM45_RS00835 in genome NZ_CP008985.1/GCF_000754265.1:

     gene            187990..189088
                     /gene="prfB"
                     /locus_tag="IM45_RS00835"
                     /old_locus_tag="IM45_349"
     CDS             join(187990..188064,188066..189088)
                     /gene="prfB"
                     /locus_tag="IM45_RS00835"
                     /old_locus_tag="IM45_349"
                     /inference="COORDINATES: similar to AA
                     sequence:RefSeq:WP_015696003.1"
                     /ribosomal_slippage
                     /note="programmed frameshift; Derived by automated
                     computational analysis using gene prediction method:
                     Protein Homology."
                     /codon_start=1
                     /transl_table=11
                     /product="peptide chain release factor 2"
                     /protein_id="WP_144241937.1"
                     /translation="MMEKHYIQNQIRNILERIKVLRGYLNYDIKKKRIEEIQAKLEIT
                     ETWNKPREAQILKQELSAIKLITNTIDNISEETIDLELLLELAFTEKDEQTFLETKSE
                     LDNIVKKLEQLEFCRMFAGQYDSANCYIDIQAGSGGTEAQDWANMLLRMYLRWSESHG
                     FKTVIIAESLGEIAGIKSVTINIIGHYAYGWLRTETGIHRLVRKSPFDSGGRRHTSFS
                     SAFVYPEINNNINININQKDLRIDVFRASGAGGQHVNRTESAVRITHLPTNIVTQCQN
                     DRSQHKNKDQAMKQLKAKLYEFELHKKNIEKITQENTKSDINWGNQIRSYVLDDSRIK
                     DLRTGIEIRNIQAVLDGNLDKFIKASLKAGL"

is detected as chimeric:
chimblock_817 GCF_000754265.1 ['IM45_RS00835', 'IM45_RS00835'] [0, 0] NZ_CP008985.1 187990.0 189160.0

[sjanssen2]

I've manually checked these additional chimeric blocks. All point to ribosomal frameshifts

[sjanssen2]

I am still thinking about chimera and wonder if we don't need to have a min overlap size of approx. k-mer length?!

assume we have to CDS A and B and read size is 4. Reads 1 to 4 cover A, reads 8 and 9 cover B.

AAAAAAA     
       BBBBB
1111        
 2222       
  3333      
   4444     
       8888 
        9999

As A and B are next to each other, but not overlapping, there is also no read that overlaps with both CDS and I'd assume no assembler will ever join these CDS.

Even if we have small overlap of say 1bp like here:

AAAAAAA     
      BBBBB
1111       
 2222      
  3333     
   4444    
      8888 
       9999

although the reads 4 and 8 both overlap A and B, there is probably no k-mer long enough to actually span this overlap.

There is the cluster option from bedtools: https://bedtools.readthedocs.io/en/latest/content/tools/cluster.html where you can also specify a "Maximum distance between features" (which can be negative) to cluster CDS into "chimera". What is the typical k-mer size of the algorithms you tested?
How about setting this neg. max distance to 16?! This should result in even less chimera!

[sjanssen2]

Back to CDS duplication. I've checked all 615 gff files of the reference dataset. Whenever a CDS ID occurs multiple times, it is always due to a "ribosomal frameshift".
You can validate this by looking at the key-value pair exception=ribosomal slippage in these instances.

[tensulin]

Kmer lengths are varied for assemblers

• velvet/ Oases: 31
• ibda-mt (=ibd_tran kmer)
• idba_tran (=30)
• idba_ud (iterative k default: starting from 20 going to 100 each step is 10; inner k refinement with k=10 and k=15)
• megahit (iterative: [21,29,39,59,79,99,119,141])
• rnaspades (auto detect, chose 25 for Marbel datasets)
• dbg: 32
• SOAPdenovo-Trans (default=23)
• transabyss (default=32)
• Trinity (default=25)

[tensulin]

resolved in commit: e22b8a8