the indel table

[thkuo]

Problems to solve:

• core_genes_50.txt
  implemented in makeGroupAln.py
• currently not able to expand by genes
  solved with "dynamic" function of snakemake
• vcf2indel.py
  keep Aaron's algorithm but make it communicate with snakemake more smoothly
• 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...