diff --git a/conda.yml b/conda.yml index eb3341a..0939745 100755 --- a/conda.yml +++ b/conda.yml @@ -4,10 +4,8 @@ channels: - conda-forge - bioconda dependencies: - - bedtools=2.29.2 - - bowtie2=2.4.1 - - minimap2=2.17 - - samtools=1.9 - - bcftools=1.9 - - bedops=2.4.39 - - tabix=0.2.6 + - bedtools + - minimap2 + - samtools + - bcftools + - tabix diff --git a/variant_remapping_tools/reads_to_remapped_variants.py b/variant_remapping_tools/reads_to_remapped_variants.py index 5cfef47..87b7b82 100755 --- a/variant_remapping_tools/reads_to_remapped_variants.py +++ b/variant_remapping_tools/reads_to_remapped_variants.py @@ -8,6 +8,7 @@ from Bio.Alphabet import generic_dna import pysam +nucleotide_alphabet = {'A', 'T', 'C', 'G'} def reverse_complement(sequence): return str(Seq(sequence, generic_dna).reverse_complement()) @@ -27,6 +28,9 @@ def calculate_new_variant_definition(left_read, right_read, ref_fasta, original_ new_ref = fetch_bases(ref_fasta, left_read.reference_name, left_read.reference_end + 1, right_read.reference_start - left_read.reference_end).upper() + if len(set(new_ref).difference(nucleotide_alphabet)) != 0 : + failure_reason = 'Reference Allele not in ACGT' + new_pos = left_read.reference_end + 1 # 1. Handle reference strand change