engreitzlab / variant-flowfish Goto Github PK

https://github.com/EngreitzLab/variant-flowfish/blob/main/workflow/rules/make_count_tables.smk#L27

We should have it pick the most frequent one (as intended?) or account for all of the unique reference sequences.

trim_count_table rule takes the head N lines of the count table file; don't think the make_count_tables function is currently sorting the output so we should make sure that is happening.

• outline:
  • for each sample:
    • match variants with .str.contains() and remove from sample
    • from the rest, choose the most frequent refAllele as the "true" reference and quantify mismatches with Aligned_Sequence
  • collect variants and references from each sample, aggregate and create count tables.
• todos:
  • add mismatch function from @engreitz
  • plot mismatches in reference

PhiX FASTA: https://www.ncbi.nlm.nih.gov/nuccore/NC_001422.1?report=fasta

Variant matching to MappingSequence currently occurs in normalize_allele_effects.py after MLE, should move this to the count table creation step.

In https://github.com/EngreitzLab/variant-flowfish/blob/main/workflow/rules/make_count_tables.smk:

• In this line -

  variant-flowfish/workflow/rules/make_count_tables.smk

  Line 27 in 1212739

  ref_seq = allele_tbl.loc[allele_tbl['Aligned_Sequence'] == allele_tbl['Reference_Sequence'], 'Reference_Sequence'].values[0]

  is it supposed to be exact matching Aligned_Sequence and Reference_Sequence? Seems like in tests, there aren't any matches which breaks the code.
  • From @lampburglar: If our amplicon is ~250 nucleotides long then the probability having an amplicon without any substitutions is .99^250 or ~8% which makes it seem like you should still be able to find an exact match amongst the 100K reads that there will be. There must be some unstable region in the amplicon that yields lots of background error.
  • I am rerunning CRISPResso and the rest of the pipeline. We were thinking that we can just check the region that there is supposed to be the edit, and if just this region is an exact match then the sequence can be considered an unedited reference sequence.
• separate out functions as another script in scripts/