Make sure that reference nucleotides we expect to see on the reference
    transcript from a variant are the same ones we encounter.
975465
def variant_matches_reference_sequence(variant, ref_seq_on_transcript, strand):
    
    if strand == "-":
        ref_seq_on_transcript = reverse_complement_dna(ref_seq_on_transcript)
    return ref_seq_on_transcript == variant.ref
785335
def check_structure_chain_quality(self, seqprop, structprop, chain_id,
                                      seq_ident_cutoff=0.5, allow_missing_on_termini=0.2,
                                      allow_mutants=True, allow_deletions=False,
                                      allow_insertions=False, allow_unresolved=True):
        """"""
        alignment = self._get_seqprop_to_structprop_alignment(seqprop=seqprop, structprop=structprop, chain_id=chain_id)

        # Compare sequence to structure's sequence using the alignment
        chain_passes_quality_check = ssbio.protein.structure.properties.quality.sequence_checker(reference_seq_aln=alignment[0],
                                                                   structure_seq_aln=alignment[1],
                                                                   seq_ident_cutoff=seq_ident_cutoff,
                                                                   allow_missing_on_termini=allow_missing_on_termini,
                                                                   allow_mutants=allow_mutants,
                                                                   allow_deletions=allow_deletions,
                                                                   allow_insertions=allow_insertions,
                                                                   allow_unresolved=allow_unresolved)
        return chain_passes_quality_check
1093236
def initial_variant_sequences_from_reads(
        variant_reads,
        max_nucleotides_before_variant=None,
        max_nucleotides_after_variant=None):
    
    unique_sequence_groups = group_unique_sequences(
        variant_reads,
        max_prefix_size=max_nucleotides_before_variant,
        max_suffix_size=max_nucleotides_after_variant)

    return [
        VariantSequence(
            prefix=prefix,
            alt=alt,
            suffix=suffix,
            reads=reads)
        for ((prefix, alt, suffix), reads)
        in unique_sequence_groups.items()
    ]
620417
def add_flanking_seqs(self, ref_seq, new_start, new_end):
        ''''''
        if new_start > self.POS or new_end < self.ref_end_pos():
            raise Error('new start and end positions must not try to shrink VCF record. new_start=' + str(new_start) + ', new_end=' + str(new_end) + '. VCF=' + str(self))

        new_start_nucleotides = ref_seq[new_start:self.POS]
        new_end_nucleotodes = ref_seq[self.ref_end_pos() + 1:new_end + 1]
        self.POS = new_start
        self.REF = new_start_nucleotides + self.REF + new_end_nucleotodes
        self.ALT = [new_start_nucleotides + x + new_end_nucleotodes for x in self.ALT]
1371184
def _guess_sequence_type_from_string(self, seq):
        ''''''
        # Define expected residues for each sequence type
        aa_chars = ['P','V','L','I','M','F','Y','W','H','K','R','Q','N','E','D','S','X','*']
        aas = set(itertools.chain(aa_chars, [lower(a) for a in aa_chars]))
        na_chars = ['A', 'T', 'G', 'C', 'N', 'U']
        nas = set(itertools.chain(na_chars, [lower(a) for a in na_chars]))

        num_nucleotide = 0
        num_protein = 0
        count = 0
        for residue in seq:
            if residue in nas:
                num_nucleotide += 1
            elif residue in aas:
                num_protein += 1
            else:
                raise Exception('Encountered unexpected character when attempting to guess sequence type: %s' % (residue))
            count += 1
            if count >300: break
        if float(num_protein) / (num_protein+num_nucleotide) > 0.1:
            return self.PROTEIN_SEQUENCE_TYPE
        else:
            return self.NUCLEOTIDE_SEQUENCE_TYPE
91845
def get_variants(data, include_germline=False):
    
    data = utils.deepish_copy(data)
    supported = ["precalled", "vardict", "vardict-java", "vardict-perl",
                 "freebayes", "octopus", "strelka2"]
    # Right now mutect2 and mutect do not provide heterozygous germline calls
    # to be useful https://github.com/bcbio/bcbio-nextgen/issues/2464
    # supported += ["mutect2", "mutect"]
    if include_germline:
        supported.insert(1, "gatk-haplotype")
    out = []
    # CWL based input
    if isinstance(data.get("variants"), dict) and "samples" in data["variants"]:
        cur_vs = []
        # Unpack single sample list of files
        if (isinstance(data["variants"]["samples"], (list, tuple)) and
              len(data["variants"]["samples"]) == 1 and isinstance(data["variants"]["samples"][0], (list, tuple))):
            data["variants"]["samples"] = data["variants"]["samples"][0]
        for fname in data["variants"]["samples"]:
            variantcaller = utils.splitext_plus(os.path.basename(fname))[0]
            variantcaller = variantcaller.replace(dd.get_sample_name(data) + "-", "")
            for batch in dd.get_batches(data):
                variantcaller = variantcaller.replace(batch + "-", "")
            cur_vs.append({"vrn_file": fname, "variantcaller": variantcaller})
        data["variants"] = cur_vs
    for v in data.get("variants", []):
        if v["variantcaller"] in supported and v.get("vrn_file"):
            out.append((supported.index(v["variantcaller"]), v))
    out.sort()
    return [xs[1] for xs in out]
647307
def ensure_annotations(resources, data):
    
    transcript_gff = tz.get_in(["rnaseq", "transcripts"], resources)
    if transcript_gff and utils.file_exists(transcript_gff):
        out_dir = os.path.join(tz.get_in(["dirs", "work"], data),
                               "inputs", "data", "annotations")
        resources["rnaseq"]["gene_bed"] = gtf.gtf_to_bed(transcript_gff, out_dir)
    return resources
1095459
def _refine_v(seq, species):
        '''
        
        '''
        vgerm = germlines.get_germline(seq['v_gene']['full'], species)
        aln = global_alignment(seq['vdj_nt'], vgerm)
        prepend = ''
        for s, g in zip(aln.aligned_query, aln.aligned_target):
            if s != '-':
                break
            else:
                prepend += g
        seq['vdj_nt'] = prepend + seq['vdj_nt']
481607
def _get_biallelic_variant(self, variant, info, _check_alleles=True):
        """"""
        info = info.iloc[0, :]
        assert not info.multiallelic

        # Seeking and parsing the file
        self._impute2_file.seek(info.seek)
        genotypes = self._parse_impute2_line(self._impute2_file.readline())

        variant_alleles = variant._encode_alleles([
            genotypes.reference, genotypes.coded,
        ])
        if (_check_alleles and variant_alleles != variant.alleles):
            # Variant with requested alleles is unavailable.
            logging.variant_not_found(variant)
            return []

        return [genotypes]
1789
def run(self, args, pipeline_command):
        

        # output that must be moved but not renamed
        consistentNaming = ['alignments/normal_dna_fix_pg_sorted.bam',
                           'alignments/normal_dna_fix_pg_sorted.bam.bai',
                           'alignments/rna_genome_sorted.bam',
                           'alignments/rna_genome_sorted.bam.bai',
                           'alignments/rna_transcriptome.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam.bai',
                           'mutations/merged/all_merged.vcf',
                           'rankboost/mhcii_rankboost_concise_results.tsv',
                           'rankboost/mhci_rankboost_concise_results.tsv',
                           ]

        # output that must be renamed as well as moved
        # map of the original name to the final name
        renamingNeeded = {'binding_predictions': 'binding_predictions.tar',
                         'expression': 'expression.tar',
                         'haplotyping': 'haplotyping.tar',
                         'peptides': 'peptides.tar',
                         'rankboost': 'rankboost.tar',
                         'reports': 'reports.tar',
                         'mutations/snpeffed/mutations.vcf': 'all_snpeffed.vcf',
                         'mutations/transgened/mutations.vcf': 'all_transgened.vcf',
                         'mutations/merged': 'merged_perchrom.tar',
                         'mutations/muse': 'muse_perchrom.tar',
                         'mutations/mutect': 'mutect_perchrom.tar',
                         'mutations/radia': 'radia_perchrom.tar',
                         'mutations/somaticsniper': 'somaticsniper_perchrom.tar',
                         'mutations/strelka/snv': 'strelka_snv_perchrom.tar',
                         'mutations/strelka/indel': 'strelka_indel_perchrom.tar'}

        def make_output(output_dir, source_dir):
            """
            :param output_dir: dir to write the output to
            :param source_dir: dir containing the directory structure to be parsed
            :return:
            """
            def make_tar(dir, tar):
                with tarfile.open(tar, "w:gz") as tar:
                    tar.add(dir)
            # the output dir is where the real output directories are written
            protect_outputs = os.listdir(source_dir)
            for protectOut in protect_outputs:
                def getName(fileName):
                    return os.path.join(os.path.join(source_dir, protectOut), fileName)
                # move individual files out
                for fileName in consistentNaming:
                    shutil.copyfile(getName(fileName), os.path.join(output_dir, os.path.basename(fileName)))
                for src, dst in renamingNeeded.iteritems():
                    if dst.endswith('.tar'):
                        make_tar(getName(src), os.path.join(output_dir, dst))
                    else:
                        shutil.copyfile(getName(src), os.path.join(output_dir, dst))
            shutil.rmtree(source_dir)


        # prepare workdir
        mount = self._prepare_mount(args)
        self._workdir = os.path.join(mount, 'Toil-' + self._name)
        # insure the pairs are in the same directory, as protect expects
        # This is made more complicated by the fact CWLTool mounts inputs into random, read-only dirs
        # to get around this we copy all inputs into their own directories that we own
        tumor_dna_dir = os.path.expanduser('~/tumorDNA')
        tumor_rna_dir = os.path.expanduser('~/tumorRNA')
        normal_dna_dir = os.path.expanduser('~/normalDNA')
        os.mkdir(tumor_dna_dir)
        os.mkdir(tumor_rna_dir)
        os.mkdir(normal_dna_dir)
        shutil.copy(args.tumor_dna, tumor_dna_dir)
        shutil.copy(args.tumor_rna, tumor_rna_dir)
        shutil.copy(args.normal_dna, normal_dna_dir)
        shutil.copy(args.tumor_dna2, tumor_dna_dir)
        shutil.copy(args.tumor_rna2, tumor_rna_dir)
        shutil.copy(args.normal_dna2, normal_dna_dir)
        args.tumor_dna = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna))
        args.tumor_dna2 = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna2))
        args.tumor_rna = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna))
        args.tumor_rna2 = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna2))
        args.normal_dna = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna))
        args.normal_dna2 = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna2))

        # prepare config
        args_dict = vars(args)
        args_dict['output_dir'] = mount
        self._config = textwrap.dedent(self._config.format(**args_dict))
        self._sample_name = args_dict["sample_name"]
        config_path = os.path.join(self._workdir, 'config')
        command = self._make_prefix(os.path.join(self._workdir, 'jobStore'),
                                          config_path,
                                          self._workdir) + pipeline_command
        if self._resume and args.resume:
            command.append('--restart')
        self._create_workdir(args)
        with open(config_path, 'w') as f:
            f.write(self._config)

        try:
            subprocess.check_call(command)
        except subprocess.CalledProcessError as e:
            print(e, file=sys.stderr)
        finally:
            log.info('Pipeline terminated, changing ownership of output files from root to user.')
            stat = os.stat(self._mount)
            subprocess.check_call(['chown', '-R', '{}:{}'.format(stat.st_uid, stat.st_gid),
                                   self._mount])
            make_output(self._mount, os.path.join(self._mount, 'output'))
            if self._no_clean and args.no_clean:
                log.info('Flag "--no-clean" was used, therefore %s was not deleted.', self._workdir)
            else:
                log.info('Cleaning up temporary directory: %s', self._workdir)
                shutil.rmtree(self._workdir)
1208257
def from_variant_and_transcript_and_sequence_key(
            cls, variant, transcript, sequence_key):
        
        # get the interbase range of offsets which capture all reference
        # bases modified by the variant
        variant_start_offset, variant_end_offset = \
            interbase_range_affected_by_variant_on_transcript(
                variant=variant,
                transcript=transcript)

        start_codon_idx = min(transcript.start_codon_spliced_offsets)

        # skip any variants which occur in the 5' UTR or overlap the start codon
        # since
        #   (1) UTR variants have unclear coding effects and
        #   (2) start-loss variants may result in a new start codon / reading
        #       frame but we're not sure which one!
        if variant_start_offset < start_codon_idx + 3:
            logger.info(
                "Skipping transcript %s for variant %s, must be after start codon",
                transcript.name,
                variant)
            return None

        stop_codon_idx = min(transcript.stop_codon_spliced_offsets)

        # skip variants which affect the 3' UTR of the transcript since
        # they don't have obvious coding effects on the protein sequence
        if variant_start_offset >= stop_codon_idx + 3:
            logger.info(
                "Skipping transcript %s for variant %s, occurs in 3' UTR",
                transcript,
                variant)
            return None

        n_prefix = len(sequence_key.sequence_before_variant_locus)
        prefix_start_idx = variant_start_offset - n_prefix
        n_bases_between_start_and_variant = variant_start_offset - start_codon_idx
        n_full_codons_before_variant = n_bases_between_start_and_variant // 3

        # if the sequence before the variant contains more bases than the
        # distance to the start codon, then by definition it must contain
        # some untranslated bases
        contains_five_prime_utr = (n_prefix > n_bases_between_start_and_variant)
        # allows for the possibility that the first base in the sequence might
        # be the first nucleotide of the start codon
        contains_start_codon = (n_prefix >= n_bases_between_start_and_variant)
        # the sequence context might only include the 2nd or 3rd bases of
        # the start codon
        overlaps_start_codon = (n_prefix > n_bases_between_start_and_variant - 3)

        if contains_start_codon:
            offset_to_first_complete_codon = start_codon_idx - prefix_start_idx
            amino_acids_before_variant = transcript.protein_sequence[:n_full_codons_before_variant]
        else:
            reading_frame = (prefix_start_idx - start_codon_idx) % 3
            offset_to_first_complete_codon = reading_frame_to_offset(reading_frame)
            n_codons_in_prefix = (n_prefix - offset_to_first_complete_codon) // 3
            amino_acids_before_variant = transcript.protein_sequence[
                n_full_codons_before_variant - n_codons_in_prefix:
                n_full_codons_before_variant]

        return cls(
            strand=sequence_key.strand,
            sequence_before_variant_locus=sequence_key.sequence_before_variant_locus,
            sequence_at_variant_locus=sequence_key.sequence_at_variant_locus,
            sequence_after_variant_locus=sequence_key.sequence_after_variant_locus,
            offset_to_first_complete_codon=offset_to_first_complete_codon,
            contains_start_codon=contains_start_codon,
            overlaps_start_codon=overlaps_start_codon,
            contains_five_prime_utr=contains_five_prime_utr,
            amino_acids_before_variant=amino_acids_before_variant)
138082
def get_alt(self, individual=0, nucleotides_only=True):
        """"""
        #not i.startswith(',') is put in to handle cases like <DEL:ME:ALU> where we have no alternate allele
        #but some reference
        if isinstance(individual, str):
            individual = self.individuals[individual]
        if nucleotides_only:
            return [self.alt[i-1].replace('.', '') for i in self.genotype[individual] if i > 0 and not self.alt[i-1].startswith('<')]
        else:
            return [self.alt[i-1].replace('.', '') for i in self.genotype[individual] if i > 0]
52131
def _variant_filtration(in_file, ref_file, vrn_files, data, filter_type,
                        hard_filter_fn):
    
    # Algorithms multiplied by number of input files to check for large enough sample sizes
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if (not config_utils.use_vqsr(algs, in_file) or
          _already_cutoff_filtered(in_file, filter_type)):
        logger.info("Skipping VQSR, using cutoff-based filers: we don't have whole genome input data")
        return hard_filter_fn(in_file, data)
    elif not _have_training_data(vrn_files):
        logger.info("Skipping VQSR, using cutoff-based filers: genome build does not have sufficient training data")
        return hard_filter_fn(in_file, data)
    else:
        sensitivities = {"INDEL": "98.0", "SNP": "99.97"}
        recal_file, tranches_file = _run_vqsr(in_file, ref_file, vrn_files,
                                              sensitivities[filter_type], filter_type, data)
        if recal_file is None:  # VQSR failed
            logger.info("VQSR failed due to lack of training data. Using cutoff-based soft filtering.")
            return hard_filter_fn(in_file, data)
        else:
            return _apply_vqsr(in_file, ref_file, recal_file, tranches_file,
                               sensitivities[filter_type], filter_type, data)
405982
def check_genes(self):
        
        gene_dna_set = set([g.dna for g in self.genes])
        assert gene_dna_set == self.dna_choices_set
1099984
def _check_dedup(data):
    
    if dd.get_analysis(data).lower() in ["rna-seq", "smallrna-seq"] or not dd.get_aligner(data):
        dup_param = utils.get_in(data, ("config", "algorithm", "mark_duplicates"), False)
    else:
        dup_param = utils.get_in(data, ("config", "algorithm", "mark_duplicates"), True)
    if dup_param and isinstance(dup_param, six.string_types):
        logger.info("Warning: bcbio no longer support explicit setting of mark_duplicate algorithm. "
                    "Using best-practice choice based on input data.")
        dup_param = True
    return dup_param
1255989
def scan_for_valid_codon(codon_span, strand, seqid, genome, type='start'):
    
    s, e = codon_span[0], codon_span[1]
    while True:
        if (type == 'start' and strand == '+') or \
            (type == 'stop' and strand == '-'):
            s, e = s - 3, e - 3
        else:
            s, e = s + 3, e + 3

        codon = _fasta_slice(genome, seqid, s, e, strand)
        is_valid = is_valid_codon(codon, type=type)
        if not is_valid:
            if type == 'start':
                ## if we are scanning upstream for a valid start codon,
                ## stop scanning when we encounter a stop
                if is_valid_codon(codon, type='stop'):
                    return (None, None)
            elif type == 'stop':
                ## if we are scanning downstream for a valid stop codon,
                ## stop scanning when we encounter a start
                if is_valid_codon(codon, type='start'):
                    return (None, None)
            continue
        break

    return (s, e)
1211037
def _check_if_ins_is_dup(self, start, insertion):
        
        is_dup = False    # assume no
        variant_start = None

        dup_candidate_start = start - len(insertion) - 1
        dup_candidate = self._ref_seq[dup_candidate_start:dup_candidate_start + len(insertion)]
        if insertion == dup_candidate:
            is_dup = True
            variant_start = dup_candidate_start + 1

        return is_dup, variant_start