the indel table

[thkuo]

Problems to solve:

• [x] core_genes_50.txt implemented in makeGroupAln.py
• [x] currently not able to expand by genes solved with "dynamic" function of snakemake
• [x] vcf2indel.py keep Aaron's algorithm but make it communicate with snakemake more smoothly
• [x] generate_indel_features.py and roary_PA14_abricate.txt, Script rewritten and the intermediate file is no longer required, because the users need to provide one more file...

Strategy:

• Understand the roles of the intermediate files
• Rewrite the external scripts
• Try to reduce the intermediate files

[thkuo]

[vcf2indel.py] The current binary labels are ambiguous and can result in incorrect ML outcomes. The 1s may mean insertions, deletions, missing values, or mixture of indels. For a certain indel feature, furthermore, is 1 1 1 0 0 equal to 0 0 0 1 1? Are two features that have same values redundant? Key questions to answer:

• Should the definition of indels be reference-based?
• How to deal with the missing values (eg. the dots in ALT and GT fields)?

[thkuo]

In the four filters of vcf2indel.py,

case 1

not existing

if ref == ".": ...

case 2

majority-based assignment

elif alt == ".":...

case 3

mixture of majority-based and reference-based methods

elif len(sample_names.loc[samples == "."]) > sample_names.shape[0]/2:...

case 4:

reference-based

else:...

[thkuo]

Only update the connection with snakemake main script (for now)

[thkuo]

Even with the original methods, the indels were incorrectly detected in this case:

>CH2500
>F1659
>CH2502
ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTG
GTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTC
TCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGC
AGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTG
GTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGCGCCGCT
ACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCG
GCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAG
GACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCG
GTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGC
CACCTGCGCCTGCAGGACTGA
>CH2522
ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTG
GTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTC
TCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGC
AGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTG
GTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGTGCCGCT
ACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCG
GCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAG
GACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCG
GTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGC
CACCTGCGCCTGCAGGACTGA
>ESP088
ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTG
GTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTC
TCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGC
AGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTG
GTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGCGCCGCT
ACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCG
GCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAG
GACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCG
GTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGC
CACCTGCGCCTGCAGGACTGA
>MHH15083
ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTG
GTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTC
TCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGC
AGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTG
GTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGCGCCGCT
ACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCG
GCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAG
GACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCG
GTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGC
CACCTGCGCCTGCAGGACTGA

and the vcf outcome was:

##fileformat=VCFv4.2
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth.">
##contig=<ID=chrUn,length=561>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  CH2500  CH2502  CH2522  ESP088  F1659   MHH15083
chrUn   0   .   ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTGGTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTCTCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGCAGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTGGTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGCGCCGCTACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCGGCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAGGACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCGGTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGCCACCTGCGCCTGCAGGACTGA   ATGAGCCGCTTTGAAATCGCCTTTTCCGGCCAGTTGGTCGCCGGCGCCCGTCCCGAGGTGGTCAAGGCCAACCTGGCCAAGCTGTTCCAGGCCGACGCGCAGCGTATCGAACTGCTGTTCTCCGGCCGCCGGGTGGTGATCAAGAACAACCTCGATGCCGCCTCCGCGGAAAAATACCGCAGCGTGCTGGAGCGAGCGGGAGCGATCGCCGTGGTCGCCGAGATGGAGGTCGAGGAGGTGGTCATGGCGCCGCCGCCTGCGCAGACGACTCCCGTGGAGGCCCCGCAGACCCGTGCCGCTACTGGTACCAGCGCGCCCGCCGGACGCTTGCAGGTAGCGCCGCGGGACGGCTACATGGCGGCGTTCGCCGAGGTCGATGCGCCGGATTTCGGCCTGGCTCCGGTAGGCGCCGACCTACAGGACGCCAAGGCCGAAGCCGAGGCGCCGAAACTCGACCTGAGCCGCTTCAGCGTCGCCCCGGTCGGTAGCGACATGGGCCAGGCACGCTCCGAGCCAGCGGCTCCGGCTCCGGACACCAGCCACCTGCGCCTGCAGGACTGA   .   .   DP=4    GT:DP   ./. 0/0:1   1/1:1   0/0:1   ./. 0/0:1

The gene was lost in two strains, but it wasn't detected by the old method. It could be listed in the gpa table though...