Allow VCF to read a 1000G file

This should allow reading of 1000G VCFs, printing warnings rather than raising errors when there are problems. I guess that's more useful for some people?

Author: hyanwong

docs/tutorial.rst

@@ -380,33 +380,32 @@ variants from chromosome 24 of ten Norwegian and French house sparrows,
         """
         # You may want to change the following line, e.g. here we allow * (a spanning
         # deletion) to be a valid allele state
-        allowed_allele_chars = set("ATGCatgc*")
+        allele_chars = set("ATGCatgc*")
         pos = 0
         for variant in vcf:  # Loop over variants, each assumed at a unique site
             if pos == variant.POS:
-                raise ValueError("Duplicate positions for variant at position", pos)
+                print(f"Duplicate entries at position {pos}, ignoring all but the first")
+                continue
             else:
                 pos = variant.POS
             if any([not phased for _, _, phased in variant.genotypes]):
                 raise ValueError("Unphased genotypes for variant at position", pos)
-            alleles = [variant.REF] + variant.ALT
+            alleles = [variant.REF.upper()] + [v.upper() for v in variant.ALT]
             ancestral = variant.INFO.get("AA", ".")  # "." means unknown
             # some VCFs (e.g. from 1000G) have many values in the AA field: take the 1st
-            ancestral = ancestral.split("|")[0]
-            if ancestral == ".":
+            ancestral = ancestral.split("|")[0].upper()
+            if ancestral == "." or ancestral == "":
                 # use the reference as ancestral, if unknown (NB: you may not want this)
-                ancestral = variant.REF
+                ancestral = variant.REF.upper()
             # Ancestral state must be first in the allele list.
             ordered_alleles = [ancestral] + list(set(alleles) - {ancestral})
             # Check we have ATCG alleles
-            for allele in ordered_alleles:
-                if len(set(allele) - allowed_allele_chars) > 0:
-                    raise ValueError(
-                        "Site at pos {pos}: allele {allele} not in {allowed_allele_chars}"
-                    )
+            for a in ordered_alleles:
+                if len(set(a) - allele_chars) > 0:
+                    print(f"Ignoring site at pos {pos}: allele {a} not in {allele_chars}")
+                    continue
             allele_index = {
-                old_index: ordered_alleles.index(allele)
-                for old_index, allele in enumerate(alleles)
+                old_index: ordered_alleles.index(a) for old_index, a in enumerate(alleles)
             }
             # Map original allele indexes to their indexes in the new alleles list.
             genotypes = [