Categorize variant by location in transcript (5'utr, exon, intron, 3'utr)

        :return "exon", "intron", "five_utr", "three_utr", "whole_gene"
        :rtype str
967563
def _get_variant_region(self):
        
        if self._var_c.posedit.pos.start.datum == Datum.CDS_END and self._var_c.posedit.pos.end.datum == Datum.CDS_END:
            result = self.T_UTR
        elif self._var_c.posedit.pos.start.base < 0 and self._var_c.posedit.pos.end.base < 0:
            result = self.F_UTR
        elif self._var_c.posedit.pos.start.base < 0 and self._var_c.posedit.pos.end.datum == Datum.CDS_END:
            result = self.WHOLE_GENE
        elif self._var_c.posedit.pos.start.offset != 0 or self._var_c.posedit.pos.end.offset != 0:
            # leave out anything intronic for now
            result = self.INTRON
        else:    # anything else that contains an exon
            result = self.EXON
        return result
103912
def _get_variant_effect(cls, variants, ref_sequence):
        ''''''
        assert len(variants) != 0

        var_types = [x.var_type for x in variants]
        if len(set(var_types)) != 1:
            return 'MULTIPLE', '.', '.'

        var_type = var_types[0]

        assert set([x.ref_name for x in variants]) == set([ref_sequence.id])
        codon_starts = [AssemblyVariants._get_codon_start(0, x.ref_start) for x in variants]
        assert len(set(codon_starts)) == 1
        codon_start = codon_starts[0]
        aa_start = codon_start // 3
        ref_codon = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start:codon_start+3])
        ref_aa = ref_codon.translate()

        if var_type == pymummer.variant.SNP:
            new_codon = list(ref_codon.seq)
            for v in variants:
                new_codon[v.ref_start - codon_start] = v.qry_base
            new_codon = pyfastaq.sequences.Fasta('new', ''.join(new_codon))
            qry_aa = new_codon.translate()

            if ref_aa.seq == qry_aa.seq:
                return ('SYN', '.', aa_start)
            elif qry_aa.seq == '*':
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            else:
                return ('NONSYN', ref_aa.seq + str(aa_start + 1) + qry_aa.seq, aa_start)
        elif var_type in [pymummer.variant.INS, pymummer.variant.DEL]:
            if len(variants) > 1:
                return 'INDELS', '.', aa_start

            var = variants[0]

            if var_type == pymummer.variant.INS:
                new_seq = pyfastaq.sequences.Fasta('seq', var.qry_base)
            else:
                new_seq = pyfastaq.sequences.Fasta('seq', var.ref_base)

            if len(new_seq) % 3 != 0:
                return ('FSHIFT', ref_aa.seq + str(aa_start + 1) + 'fs', aa_start)

            new_seq_aa = new_seq.translate()
            if '*' in new_seq_aa.seq:
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            elif var_type == pymummer.variant.INS:
                ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                aa_after_ins = ref_codon_after_ins.translate()
                return ('INS', ref_aa.seq + str(aa_start + 1) + '_' + aa_after_ins.seq + str(aa_start + 2) + 'ins' + new_seq_aa.seq , aa_start)
            else:
                if len(new_seq) == 3:
                    return ('DEL', ref_aa.seq + str(aa_start + 1) + 'del', aa_start)
                else:
                    assert len(new_seq) % 3 == 0
                    ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                    aa_after_ins = ref_codon_after_ins.translate()
                    return ('DEL', ref_aa.seq + str(aa_start + 1)+ '_' + aa_after_ins.seq + str(aa_start + 2) + 'del', aa_start)

        else:
            return ('UNKNOWN', '.', aa_start)
691958
def translateIndelGenotypeFromVCF(vcf_genotypes, ref):
    ''''''

    # indels
    def getPrefix(s1, s2):
        '''get common prefix of strings s1 and s2.'''
        n = min(len(s1), len(s2))
        for x in range(n):
            if s1[x] != s2[x]:
                return s1[:x]
        return s1[:n]

    def getSuffix(s1, s2):
        '''get common sufix of strings s1 and s2.'''
        n = min(len(s1), len(s2))
        if s1[-1] != s2[-1]:
            return ""
        for x in range(-2, -n - 1, -1):
            if s1[x] != s2[x]:
                return s1[x + 1:]
        return s1[-n:]

    def getGenotype(variant, ref):

        if variant == ref:
            return "*", 0

        if len(ref) > len(variant):
            # is a deletion
            if ref.startswith(variant):
                return "-%s" % ref[len(variant):], len(variant) - 1
            elif ref.endswith(variant):
                return "-%s" % ref[:-len(variant)], -1
            else:
                prefix = getPrefix(ref, variant)
                suffix = getSuffix(ref, variant)
                shared = len(prefix) + len(suffix) - len(variant)
                # print "-", prefix, suffix, ref, variant, shared, len(prefix), len(suffix), len(ref)
                if shared < 0:
                    raise ValueError()
                return "-%s" % ref[len(prefix):-(len(suffix) - shared)], len(prefix) - 1

        elif len(ref) < len(variant):
            # is an insertion
            if variant.startswith(ref):
                return "+%s" % variant[len(ref):], len(ref) - 1
            elif variant.endswith(ref):
                return "+%s" % variant[:len(ref)], 0
            else:
                prefix = getPrefix(ref, variant)
                suffix = getSuffix(ref, variant)
                shared = len(prefix) + len(suffix) - len(ref)
                if shared < 0:
                    raise ValueError()

                return "+%s" % variant[len(prefix):-(len(suffix) - shared)], len(prefix)
        else:
            assert 0, "snp?"

    # in pileup, the position refers to the base
    # after the coordinate, hence subtract 1
    # pos -= 1

    genotypes, offsets = [], []
    is_error = True

    for variant in vcf_genotypes:
        try:
            g, offset = getGenotype(variant, ref)
        except ValueError:
            break

        genotypes.append(g)
        if g != "*":
            offsets.append(offset)

    else:
        is_error = False

    if is_error:
        raise ValueError()

    assert len(set(offsets)) == 1, "multiple offsets for indel"
    offset = offsets[0]

    genotypes = "/".join(genotypes)
    return genotypes, offset
494180
def get_transcript_source_set(gtf):
    
    gene_to_source = get_gene_source_set(gtf)
    transcript_to_source = {}
    db = get_gtf_db(gtf)
    for feature in complete_features(db):
        gene_id = feature['gene_id'][0]
        transcript_to_source[feature['transcript_id'][0]] = gene_to_source[gene_id]
    return transcript_to_source
632753
def get_start_stops(transcript_sequence, start_codons=None, stop_codons=None):
        
        transcript_sequence = transcript_sequence.upper()  # for comparison with codons below
        if not start_codons:
            start_codons = ['ATG']
        if not stop_codons:
            stop_codons = ['TAA', 'TAG', 'TGA']

        seq_frames = {1: {'starts': [], 'stops': []},
                      2: {'starts': [], 'stops': []},
                      3: {'starts': [], 'stops': []}}

        for codons, positions in ((start_codons, 'starts'),
                                  (stop_codons, 'stops')):
            if len(codons) > 1:
                pat = re.compile('|'.join(codons))
            else:
                pat = re.compile(codons[0])

            for m in re.finditer(pat, transcript_sequence):
                # Increment position by 1, Frame 1 starts at position 1 not 0
                start = m.start() + 1
                rem = start % 3
                if rem == 1:  # frame 1
                    seq_frames[1][positions].append(start)
                elif rem == 2:  # frame 2
                    seq_frames[2][positions].append(start)
                elif rem == 0:  # frame 3
                    seq_frames[3][positions].append(start)
        return seq_frames
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
1321978
def variant_filtration(call_file, ref_file, vrn_files, data, items):
    
    caller = data["config"]["algorithm"].get("variantcaller")
    if "gvcf" not in dd.get_tools_on(data):
        call_file = ploidy.filter_vcf_by_sex(call_file, items)
    if caller in ["freebayes"]:
        return vfilter.freebayes(call_file, ref_file, vrn_files, data)
    elif caller in ["platypus"]:
        return vfilter.platypus(call_file, data)
    elif caller in ["samtools"]:
        return vfilter.samtools(call_file, data)
    elif caller in ["gatk", "gatk-haplotype", "haplotyper"]:
        if dd.get_analysis(data).lower().find("rna-seq") >= 0:
            from bcbio.rnaseq import variation as rnaseq_variation
            return rnaseq_variation.gatk_filter_rnaseq(call_file, data)
        else:
            return gatkfilter.run(call_file, ref_file, vrn_files, data)
    # no additional filtration for callers that filter as part of call process
    else:
        return call_file
34558
def transcript_to_gene(gtf):
    
    gene_lookup = {}
    for feature in complete_features(get_gtf_db(gtf)):
        gene_id = feature.attributes.get('gene_id', [None])[0]
        transcript_id = feature.attributes.get('transcript_id', [None])[0]
        gene_lookup[transcript_id] = gene_id
    return gene_lookup
705163
def _retranslate(seq):
        '''
        
        '''
        if len(seq['vdj_nt']) % 3 != 0:
            trunc = len(seq['vdj_nt']) % 3
            seq['vdj_nt'] = seq['vdj_nt'][:-trunc]
        seq['vdj_aa'] = Seq(seq['vdj_nt'], generic_dna).translate()
257105
def get_sequence(self,ref):
    """"""
    strand = '+'
    if not self._options.direction:
      sys.stderr.write("WARNING: no strand information for the transcript\n")
    if self._options.direction: strand = self._options.direction
    seq = ''
    for e in [x.range for x in self.exons]:
      seq += str(ref[e.chr][e.start-1:e.end])
    if strand == '-':  seq = rc(seq)
    return Sequence(seq.upper())
976365
def get_transcript_ids_for_ensembl_gene_ids(self, gene_ids, hgnc_symbols):
        
        
        chroms = {"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", \
             "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", \
              "X", "Y"}
        
        headers = {"content-type": "application/json"}
        
        transcript_ids = []
        for gene_id in gene_ids:
            self.attempt = 0
            ext = "/overlap/id/{}?feature=transcript".format(gene_id)
            r = self.ensembl_request(ext, headers)
            
            for item in json.loads(r):
                # ignore non-coding transcripts
                if item["biotype"] not in ["protein_coding", "polymorphic_pseudogene"]:
                    continue
                
                # ignore transcripts not on the standard chromosomes
                # (non-default chroms fail to map the known de novo variants
                # to the gene location
                if item["Parent"] != gene_id or item["seq_region_name"] not in \
                        chroms or \
                        all([symbol not in item["external_name"] for symbol in hgnc_symbols]):
                    continue
                transcript_ids.append(item["id"])
        
        return transcript_ids
958671
def _init_exons(self):
        
        exon_starts = [self.chrom_start + int(s)
                       for s in self.bed_tuple.blockStarts.strip(',').split(',')]
        exon_sizes = list(map(int, self.bed_tuple.blockSizes.strip(',').split(',')))

        # get chromosome intervals
        exons = [(exon_starts[i], exon_starts[i] + exon_sizes[i])
                 for i in range(len(exon_starts))]
        no_utr_exons = self._filter_utr(exons)
        self.exons = no_utr_exons
        self.exon_lens = [e[1] - e[0] for e in self.exons]
        self.num_exons = len(self.exons)
        self.cds_len = sum(self.exon_lens)
        self.five_ss_len = 2*(self.num_exons-1)
        self.three_ss_len = 2*(self.num_exons-1)
        self._init_splice_site_pos()
776694
def parse_reports(self):
    

    self.qualimap_rnaseq_genome_results = dict()
    regexes = {
        'reads_aligned': r"read(?:s| pairs) aligned\s*=\s*([\d,]+)",
        'total_alignments': r"total alignments\s*=\s*([\d,]+)",
        'non_unique_alignments': r"non-unique alignments\s*=\s*([\d,]+)",
        'reads_aligned_genes': r"aligned to genes\s*=\s*([\d,]+)",
        'ambiguous_alignments': r"ambiguous alignments\s*=\s*([\d,]+)",
        'not_aligned': r"not aligned\s*=\s*([\d,]+)",
        '5_3_bias': r"5'-3' bias\s*=\s*([\d,\.]+)$",
        'reads_aligned_exonic': r"exonic\s*=\s*([\d,]+)",
        'reads_aligned_intronic': r"intronic\s*=\s*([\d,]+)",
        'reads_aligned_intergenic': r"intergenic\s*=\s*([\d,]+)",
        'reads_aligned_overlapping_exon': r"overlapping exon\s*=\s*([\d,]+)",
    }
    for f in self.find_log_files('qualimap/rnaseq/rnaseq_results'):
        d = dict()

        # Get the sample name
        s_name_regex = re.search(r"bam file\s*=\s*(.+)", f['f'], re.MULTILINE)
        if s_name_regex:
            d['bam_file'] = s_name_regex.group(1)
            s_name = self.clean_s_name(d['bam_file'], f['root'])
        else:
            log.warn("Couldn't find an input filename in genome_results file {}/{}".format(f['root'], f['fn']))
            return None

        # Check for and 'fix' European style decimal places / thousand separators
        comma_regex = re.search(r"exonic\s*=\s*[\d\.]+ \(\d{1,3},\d+%\)", f['f'], re.MULTILINE)
        if comma_regex:
            log.debug("Trying to fix European comma style syntax in Qualimap report {}/{}".format(f['root'], f['fn']))
            f['f'] = f['f'].replace('.','')
            f['f'] = f['f'].replace(',','.')

        # Go through all numeric regexes
        for k, r in regexes.items():
            r_search = re.search(r, f['f'], re.MULTILINE)
            if r_search:
                try:
                    d[k] = float(r_search.group(1).replace(',',''))
                except UnicodeEncodeError:
                    # Qualimap reports infinity (\u221e) when 3' bias denominator is zero
                    pass
                except ValueError:
                    d[k] = r_search.group(1)

        # Add to general stats table
        for k in ['5_3_bias', 'reads_aligned']:
            try:
                self.general_stats_data[s_name][k] = d[k]
            except KeyError:
                pass

        # Save results
        if s_name in self.qualimap_rnaseq_genome_results:
            log.debug("Duplicate genome results sample name found! Overwriting: {}".format(s_name))
        self.qualimap_rnaseq_genome_results[s_name] = d
        self.add_data_source(f, s_name=s_name, section='rna_genome_results')


    #### Coverage profile
    self.qualimap_rnaseq_cov_hist = dict()
    for f in self.find_log_files('qualimap/rnaseq/coverage', filehandles=True):
        s_name = self.get_s_name(f)
        d = dict()
        for l in f['f']:
            if l.startswith('#'):
                continue
            coverage, count = l.split(None, 1)
            coverage = int(round(float(coverage)))
            count = float(count)
            d[coverage] = count

        if len(d) == 0:
            log.debug("Couldn't parse contents of coverage histogram file {}".format(f['fn']))
            return None

        # Save results
        if s_name in self.qualimap_rnaseq_cov_hist:
            log.debug("Duplicate coverage histogram sample name found! Overwriting: {}".format(s_name))
        self.qualimap_rnaseq_cov_hist[s_name] = d
        self.add_data_source(f, s_name=s_name, section='rna_coverage_histogram')

    # Filter to strip out ignored sample names
    self.qualimap_rnaseq_genome_results = self.ignore_samples(self.qualimap_rnaseq_genome_results)
    self.qualimap_rnaseq_cov_hist = self.ignore_samples(self.qualimap_rnaseq_cov_hist)

    #### Plots

    # Genomic Origin Bar Graph
    # NB: Ignore 'Overlapping Exon' in report - these make the numbers add up to > 100%
    if len(self.qualimap_rnaseq_genome_results) > 0:
        gorigin_cats = OrderedDict()
        gorigin_cats['reads_aligned_exonic'] = {'name': 'Exonic'}
        gorigin_cats['reads_aligned_intronic'] = {'name': 'Intronic'}
        gorigin_cats['reads_aligned_intergenic'] = {'name': 'Intergenic'}
        gorigin_pconfig = {
            'id': 'qualimap_genomic_origin',
            'title': 'Qualimap RNAseq: Genomic Origin',
            'ylab': 'Number of reads',
            'cpswitch_c_active': False
        }
        genomic_origin_helptext = '''
        There are currently three main approaches to map reads to transcripts in an
        RNA-seq experiment: mapping reads to a reference genome to identify expressed
        transcripts that are annotated (and discover those that are unknown), mapping
        reads to a reference transcriptome, and <i>de novo</i> assembly of transcript
        sequences (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).

        For RNA-seq QC analysis, QualiMap can be used to assess alignments produced by
        the first of these approaches. For input, it requires a GTF annotation file
        along with a reference genome, which can be used to reconstruct the exon
        structure of known transcripts. This allows mapped reads to be grouped by
        whether they originate in an exonic region (for QualiMap, this may include
        5&#8242; and 3&#8242; UTR regions as well as protein-coding exons), an intron,
        or an intergenic region (see the <a href="http://qualimap.bioinfo.cipf.es/doc_html/index.html"
        target="_blank">Qualimap 2 documentation</a>).

        The inferred genomic origins of RNA-seq reads are presented here as a bar graph
        showing either the number or percentage of mapped reads in each read dataset
        that have been assigned to each type of genomic region. This graph can be used
        to assess the proportion of useful reads in an RNA-seq experiment. That
        proportion can be reduced by the presence of intron sequences, especially if
        depletion of ribosomal RNA was used during sample preparation (<a href="https://doi.org/10.1038/nrg3642"
        target="_blank">Sims et al. 2014</a>). It can also be reduced by off-target
        transcripts, which are detected in greater numbers at the sequencing depths
        needed to detect poorly-expressed transcripts (<a href="https://doi.org/10.1101/gr.124321.111"
        target="_blank">Tarazona et al. 2011</a>).'''
        self.add_section (
            name = 'Genomic origin of reads',
            anchor = 'qualimap-reads-genomic-origin',
            description = 'Classification of mapped reads as originating in exonic, intronic or intergenic regions. These can be displayed as either the number or percentage of mapped reads.',
            helptext = genomic_origin_helptext,
            plot = bargraph.plot(self.qualimap_rnaseq_genome_results, gorigin_cats, gorigin_pconfig)
        )

    if len(self.qualimap_rnaseq_cov_hist) > 0:
        coverage_profile_helptext = '''
        There are currently three main approaches to map reads to transcripts in an
        RNA-seq experiment: mapping reads to a reference genome to identify expressed
        transcripts that are annotated (and discover those that are unknown), mapping
        reads to a reference transcriptome, and <i>de novo</i> assembly of transcript
        sequences (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).

        For RNA-seq QC analysis, QualiMap can be used to assess alignments produced by
        the first of these approaches. For input, it requires a GTF annotation file
        along with a reference genome, which can be used to reconstruct the exon
        structure of known transcripts. QualiMap uses this information to calculate the
        depth of coverage along the length of each annotated transcript. For a set of
        reads mapped to a transcript, the depth of coverage at a given base position is
        the number of high-quality reads that map to the transcript at that position
        (<a href="https://doi.org/10.1038/nrg3642" target="_blank">Sims et al. 2014</a>).

        QualiMap calculates coverage depth at every base position of each annotated
        transcript. To enable meaningful comparison between transcripts, base positions
        are rescaled to relative positions expressed as percentage distance along each
        transcript (*0%, 1%, &#8230;, 99%*). For the set of transcripts with at least
        one mapped read, QualiMap plots the cumulative mapped-read depth (y-axis) at
        each relative transcript position (x-axis). This plot shows the gene coverage
        profile across all mapped transcripts for each read dataset. It provides a
        visual way to assess positional biases, such as an accumulation of mapped reads
        at the 3&#8242; end of transcripts, which may indicate poor RNA quality in the
        original sample (<a href="https://doi.org/10.1186/s13059-016-0881-8"
        target="_blank">Conesa et al. 2016</a>).'''
        self.add_section (
            name = 'Gene Coverage Profile',
            anchor = 'qualimap-genome-fraction-coverage',
            description = 'Mean distribution of coverage depth across the length of all mapped transcripts.',
            helptext = coverage_profile_helptext,
            plot = linegraph.plot(self.qualimap_rnaseq_cov_hist, {
                'id': 'qualimap_gene_coverage_profile',
                'title': 'Qualimap RNAseq: Coverage Profile Along Genes (total)',
                'ylab': 'Coverage',
                'xlab': 'Transcript Position (%)',
                'ymin': 0,
                'xmin': 0,
                'xmax': 100,
                'tt_label': '<b>{point.x} bp</b>: {point.y:.0f}%',
            })
        )


    #### General Stats
    self.general_stats_headers['5_3_bias'] = {
        'title': "5'-3' bias",
        'format': '{:,.2f}',
    }
    self.general_stats_headers['reads_aligned'] = {
        'title': '{} Aligned'.format(config.read_count_prefix),
        'description': 'Reads Aligned ({})'.format(config.read_count_desc),
        'min': 0,
        'scale': 'RdBu',
        'shared_key': 'read_count',
        'modify': lambda x: x * config.read_count_multiplier
    }

    # Return the number of reports we found
    return len(self.qualimap_rnaseq_genome_results.keys())
79767
def get_transcript(self,exon_bounds='max'):
    """"""
    out = Transcript()
    out.junctions = [x.get_junction() for x in self.junction_groups]
    # check for single exon transcript
    if len(out.junctions) == 0:
      leftcoord = min([x.exons[0].range.start for x in self.transcripts])
      rightcoord = max([x.exons[-1].range.end for x in self.transcripts])
      e = Exon(GenomicRange(x.exons[0].range.chr,leftcoord,rightcoord))
      e.set_is_leftmost()
      e.set_is_rightmost()
      out.exons.append(e)
      return out
    # get internal exons
    self.exons = []
    for i in range(0,len(self.junction_groups)-1):
      j1 = self.junction_groups[i].get_junction()
      j2 = self.junction_groups[i+1].get_junction()
      e = Exon(GenomicRange(j1.right.chr,j1.right.end,j2.left.start))
      e.set_left_junc(j1)
      e.set_right_junc(j2)
      #print str(i)+" to "+str(i+1)
      out.exons.append(e)
    # get left exon
    left_exons = [y for y in [self.transcripts[e[0]].junctions[e[1]].get_left_exon() for e in self.junction_groups[0].evidence] if y]
    if len(left_exons) == 0:
      sys.stderr.write("ERROR no left exon\n")
      sys.exit()
    e_left = Exon(GenomicRange(out.junctions[0].left.chr,\
                               min([x.range.start for x in left_exons]),
                               out.junctions[0].left.start))
    e_left.set_right_junc(out.junctions[0])
    out.exons.insert(0,e_left)
    # get right exon
    right_exons = [y for y in [self.transcripts[e[0]].junctions[e[1]].get_right_exon() for e in self.junction_groups[-1].evidence] if y]
    if len(right_exons) == 0:
      sys.stderr.write("ERROR no right exon\n")
      sys.exit()
    e_right = Exon(GenomicRange(out.junctions[-1].right.chr,\
                               out.junctions[-1].right.end,\
                               max([x.range.end for x in right_exons])))
    e_right.set_left_junc(out.junctions[-1])
    out.exons.append(e_right)
    return out
68603
def collapse_substrings(variant_sequences):
    
    if len(variant_sequences) <= 1:
        # if we don't have at least two VariantSequences then just
        # return your input
        return variant_sequences

    # dictionary mapping VariantSequence objects to lists of reads
    # they absorb from substring VariantSequences
    extra_reads_from_substrings = defaultdict(set)
    result_list = []
    # sort by longest to shortest total length
    for short_variant_sequence in sorted(
            variant_sequences,
            key=lambda seq: -len(seq)):
        found_superstring = False
        for long_variant_sequence in result_list:
            found_superstring = long_variant_sequence.contains(short_variant_sequence)
            if found_superstring:
                extra_reads_from_substrings[long_variant_sequence].update(
                    short_variant_sequence.reads)
        if not found_superstring:
            result_list.append(short_variant_sequence)
    # add to each VariantSequence the reads it absorbed from dropped substrings
    # and then return
    return [
        variant_sequence.add_reads(
            extra_reads_from_substrings[variant_sequence])
        for variant_sequence in result_list
    ]
805578
def from23andme(args):
    
    p = OptionParser(from23andme.__doc__)
    p.set_ref()
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    txtfile, seqid = args
    ref_dir = opts.ref
    fastafile = op.join(ref_dir, "hs37d5.fa")
    fasta = Fasta(fastafile)

    pf = txtfile.rsplit(".", 1)[0]
    px = CM[seqid]
    chrvcf = pf + ".{0}.vcf".format(px)
    legend = op.join(ref_dir, "1000GP_Phase3/{0}.rsids".format(px))
    register = read_rsid(seqid, legend)

    fw = open(chrvcf, "w")
    print(get_vcfstanza(fastafile, fasta, txtfile), file=fw)

    fp = open(txtfile)
    seen = set()
    duplicates = skipped = missing = 0
    for row in fp:
        if row[0] == '#':
            continue
        rsid, chr, pos, genotype = row.split()
        if chr != seqid:
            continue
        pos = int(pos)
        if (chr, pos) in seen:
            duplicates += 1
            continue
        seen.add((chr, pos))
        genotype = list(genotype)
        if "-" in genotype:  # missing daa
            missing += 1
            continue

        # Y or MT
        if not register:
            assert len(genotype) == 1
            ref = fasta[chr][pos - 1].seq.upper()
            if "D" in genotype or "I" in genotype:
                skipped += 1
                continue
            genotype = genotype[0]
            code = "0/0" if ref == genotype else "1/1"
            alt = "." if ref == genotype else genotype
            print("\t".join(str(x) for x in \
                    (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)
            continue

        # If rsid is seen in the db, use that
        if rsid in register:
            pos, ref, alt = register[rsid]
        elif pos in register:
            pos, ref, alt = register[pos]
        else:
            skipped += 1  # Not in reference panel
            continue

        assert fasta[chr][pos - 1:pos + len(ref) - 1].seq.upper() == ref
        # Keep it bi-allelic
        not_seen = [x for x in alt if x not in genotype]
        while len(alt) > 1 and not_seen:
            alt.remove(not_seen.pop())
        if len(alt) > 1:
            alt = [alt[0]]
        alleles = [ref] + alt

        if len(genotype) == 1:
            genotype = [genotype[0]] * 2

        alt = ",".join(alt) or "."
        if "D" in genotype or "I" in genotype:
            max_allele = max((len(x), x) for x in alleles)[1]
            alleles = [("I" if x == max_allele else "D") for x in alleles]
            assert "I" in alleles and "D" in alleles
        a, b = genotype
        try:
            ia, ib = alleles.index(a), alleles.index(b)
        except ValueError:  # alleles not seen
            logging.error("{0}: alleles={1}, genotype={2}".\
                            format(rsid, alleles, genotype))
            skipped += 1
            continue
        code = "/".join(str(x) for x in sorted((ia, ib)))

        print("\t".join(str(x) for x in \
                (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)

    logging.debug("duplicates={0} skipped={1} missing={2}".\
                    format(duplicates, skipped, missing))
1291224
def fixpartials(args):
    
    p = OptionParser(fixpartials.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    gffile, gfasta, partials, = args

    gff = make_index(gffile)
    genome = Fasta(gfasta, index=True)
    partials = LineFile(partials, load=True).lines

    # all_transcripts = [f.id for f in gff.features_of_type("mRNA", \
    #     order_by=("seqid", "start"))]
    seen = set()
    fw = must_open(opts.outfile, "w")
    for gene in gff.features_of_type('gene',  order_by=('seqid', 'start')):
        children = AutoVivification()
        cflag = False
        transcripts = list(gff.children(gene, level=1, order_by=('start')))
        for transcript in transcripts:
            trid, seqid, strand = transcript.id, transcript.seqid, transcript.strand

            for child in gff.children(transcript, order_by=('start')):
                ftype = child.featuretype
                if ftype not in children[trid]: children[trid][ftype] = []
                children[trid][ftype].append(child)

            five_prime, three_prime = True, True
            nstart, nstop = (None, None), (None, None)
            cds_span = [children[trid]['CDS'][0].start, \
                children[trid]['CDS'][-1].stop]
            new_cds_span = [x for x in cds_span]

            start_codon = (cds_span[0], cds_span[0] + 2)
            stop_codon = (cds_span[1] - 2, cds_span[1])
            if strand == '-': start_codon, stop_codon = stop_codon, start_codon

            if trid in partials:
                seen.add(trid)
                start_codon_fasta = _fasta_slice(genome, seqid, \
                    start_codon[0], start_codon[1], strand)
                stop_codon_fasta = _fasta_slice(genome, seqid, \
                    stop_codon[0], stop_codon[1], strand)

                if not is_valid_codon(start_codon_fasta, type='start'):
                    five_prime = False
                    nstart = scan_for_valid_codon(start_codon, strand, \
                        seqid, genome, type='start')

                if not is_valid_codon(stop_codon_fasta, type='stop'):
                    three_prime = False
                    nstop = scan_for_valid_codon(stop_codon, strand, \
                        seqid, genome, type='stop')

                logging.debug("feature={0} ({1})".format(trid, strand) + \
                ", 5'={0}, 3'={1}".format(five_prime, three_prime) + \
                ", {0} <== {1} ==> {2}".format(nstart if strand == '+' else nstop, \
                    cds_span, nstop if strand == '+' else nstart))

            if not five_prime or not three_prime:
                if nstart != (None, None) and (start_codon != nstart):
                    i = 0 if strand == '+' else 1
                    new_cds_span[i] = nstart[i]
                if nstop != (None, None) and (stop_codon != nstop):
                    i = 1 if strand == '+' else 0
                    new_cds_span[i] = nstop[i]
                new_cds_span.sort()

            if set(cds_span) != set(new_cds_span):
                cflag = True
                # if CDS has been extended, appropriately adjust all relevent
                # child feature (CDS, exon, UTR) coordinates
                for ftype in children[trid]:
                    for idx in range(len(children[trid][ftype])):
                        child_span = (children[trid][ftype][idx].start, \
                            children[trid][ftype][idx].stop)
                        if ftype in ('exon', 'CDS'):
                            # if exons/CDSs, adjust start and stop according to
                            # new CDS start and stop, respectively
                            if child_span[0] == cds_span[0]:
                                children[trid][ftype][idx].start = new_cds_span[0]
                            if child_span[1] == cds_span[1]:
                                children[trid][ftype][idx].stop = new_cds_span[1]
                        elif ftype.endswith('UTR'):
                            # if *_prime_UTR, adjust stop according to new CDS start and
                            #                 adjust start according to new CDS stop
                            if child_span[1] == cds_span[0]:
                                children[trid][ftype][idx].stop = new_cds_span[0]
                            if child_span[0] == cds_span[1]:
                                children[trid][ftype][idx].start = new_cds_span[1]

                transcript.start, transcript.stop = \
                    children[trid]['exon'][0].start, children[trid]['exon'][-1].stop

        if cflag:
            _gene_span = range_minmax([(tr.start, tr.stop) for tr in transcripts])
            gene.start, gene.stop = _gene_span[0], _gene_span[1]

        # print gff file
        print(gene, file=fw)
        for transcript in transcripts:
            trid = transcript.id
            print(transcript, file=fw)
            for cftype in children[trid]:
                for child in children[trid][cftype]:
                    print(child, file=fw)

    fw.close()