Run BLAST analyses of the subsampled FASTQ reads against the NCBI 16S reference database
356290
def blast(self):
        
        logging.info('BLASTing FASTA files against {} database'.format(self.analysistype))
        for _ in range(self.cpus):
            threads = Thread(target=self.blastthreads, args=())
            threads.setDaemon(True)
            threads.start()
        with progressbar(self.runmetadata.samples) as bar:
            for sample in bar:
                if sample.general.bestassemblyfile != 'NA':
                    # Set the name of the BLAST report
                    sample[self.analysistype].blastreport = os.path.join(
                        sample[self.analysistype].outputdir,
                        '{}_{}_blastresults.csv'.format(sample.name, self.analysistype))
                    # Use the NCBI BLASTn command line wrapper module from BioPython to set the parameters of the search
                    blastn = NcbiblastnCommandline(query=sample[self.analysistype].fasta,
                                                   db=os.path.splitext(sample[self.analysistype].baitfile)[0],
                                                   max_target_seqs=1,
                                                   num_threads=self.threads,
                                                   outfmt="'6 qseqid sseqid positive mismatch gaps evalue "
                                                          "bitscore slen length qstart qend qseq sstart send sseq'",
                                                   out=sample[self.analysistype].blastreport)
                    # Add a string of the command to the metadata object
                    sample[self.analysistype].blastcall = str(blastn)
                    # Add the object and the command to the BLAST queue
                    self.blastqueue.put((sample, blastn))
        self.blastqueue.join()
211318
def fastaAlignmentRead(fasta, mapFn=(lambda x : x), l=None):
    
    if l is None:
        l = _getMultiFastaOffsets(fasta)
    else:
        l = l[:]
    seqNo = len(l)
    for i in xrange(0, seqNo):
        j = open(fasta, 'r')
        j.seek(l[i])
        l[i] = j
    column = [sys.maxint]*seqNo
    if seqNo != 0:
        while True:
            for j in xrange(0, seqNo):
                i = l[j].read(1)
                while i == '\n':
                    i = l[j].read(1)
                column[j] = i
            if column[0] == '>' or column[0] == '':
                for j in xrange(1, seqNo):
                    assert column[j] == '>' or column[j] == ''
                break
            for j in xrange(1, seqNo):
                 assert column[j] != '>' and column[j] != ''
                 column[j] = mapFn(column[j])
            yield column[:]
    for i in l:
        i.close()
984494
def fasttagcount(sam, out, genemap, positional, minevidence, cb_histogram, 
                 cb_cutoff, subsample, parse_tags, gene_tags, umi_matrix):
    ''' 
    '''
    from pysam import AlignmentFile

    from io import StringIO
    import pandas as pd

    from utils import weigh_evidence

    if sam.endswith(".sam"):
        logger.error("To use the fasttagcount subcommand, the alignment file must be a "
                     "coordinate sorted, indexed BAM file.")
        sys.exit(1)

    logger.info('Reading optional files')

    gene_map = None
    if genemap:
        with open(genemap) as fh:
            try:
                gene_map = dict(p.strip().split() for p in fh)
            except ValueError:
                logger.error('Incorrectly formatted gene_map, need to be tsv.')
                sys.exit()

    if positional:
        tuple_template = '{0},{1},{2},{3}'
    else:
        tuple_template = '{0},{1},{3}'

    if not cb_cutoff:
        cb_cutoff = 0

    if cb_histogram and cb_cutoff == "auto":
        cb_cutoff = guess_depth_cutoff(cb_histogram)

    cb_cutoff = int(cb_cutoff)

    cb_hist = None
    filter_cb = False
    if cb_histogram:
        cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
        total_num_cbs = cb_hist.shape[0]
        cb_hist = cb_hist[cb_hist > cb_cutoff]
        logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
        filter_cb = True

    parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')

    if subsample:
        logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
        start_sampling  = time.time()

        reservoir = collections.defaultdict(list)
        cb_hist_sampled = 0 * cb_hist
        cb_obs = 0 * cb_hist

        track = stream_bamfile(sam)
        current_read = 'none_observed_yet'
        for i, aln in enumerate(track):
            if aln.qname == current_read:
                continue

            current_read = aln.qname

            if parse_tags:
                CB = aln.get_tag('CR')
            else:
                match = parser_re.match(aln.qname)
                CB = match.group('CB')

            if CB not in cb_hist.index:
                continue

            cb_obs[CB] += 1
            if len(reservoir[CB]) < subsample:
                reservoir[CB].append(i)
                cb_hist_sampled[CB] += 1
            else:
                s = pd.np.random.randint(0, cb_obs[CB])
                if s < subsample:
                    reservoir[CB][s] = i

        index_filter = set(itertools.chain.from_iterable(reservoir.values()))
        sam_file.close()
        sampling_time = time.time() - start_sampling
        logger.info('Sampling done - {:.3}s'.format(sampling_time))

    evidence = collections.defaultdict(lambda: collections.defaultdict(float))
    bare_evidence = collections.defaultdict(float)
    logger.info('Tallying evidence')
    start_tally = time.time()

    sam_mode = 'r' if sam.endswith(".sam") else 'rb'
    sam_file = AlignmentFile(sam, mode=sam_mode)
    transcript_map = collections.defaultdict(set)
    sam_transcripts = [x["SN"] for x in sam_file.header["SQ"]]
    if gene_map:
        for transcript, gene in gene_map.items():
            if transcript in sam_transcripts:
                transcript_map[gene].add(transcript)
    else:
        for transcript in sam_transcripts:
            transcript_map[transcript].add(transcript)
    missing_transcripts = set()
    alignments_processed = 0
    unmapped = 0
    kept = 0
    nomatchcb = 0
    current_read = 'none_observed_yet'
    current_transcript = None
    count_this_read = True
    transcripts_processed = 0
    genes_processed = 0
    cells = list(cb_hist.index)
    targets_seen = set()

    if umi_matrix:
        bare_evidence_handle = open(umi_matrix, "w")
        bare_evidence_handle.write(",".join(["gene"] + cells) + "\n")

    with open(out, "w") as out_handle:
        out_handle.write(",".join(["gene"] + cells) + "\n")
        for gene, transcripts in transcript_map.items():
            for transcript in transcripts:
                for aln in sam_file.fetch(transcript):
                    alignments_processed += 1

                    if aln.is_unmapped:
                        unmapped += 1
                        continue

                    if gene_tags and not aln.has_tag('GX'):
                        unmapped += 1
                        continue

                    if aln.qname != current_read:
                        current_read = aln.qname
                        if subsample and i not in index_filter:
                            count_this_read = False
                            continue
                        else:
                            count_this_read = True
                    else:
                        if not count_this_read:
                            continue

                    if parse_tags:
                        CB = aln.get_tag('CR')
                    else:
                        match = parser_re.match(aln.qname)
                        CB = match.group('CB')

                    if filter_cb:
                        if CB not in cb_hist.index:
                            nomatchcb += 1
                            continue

                    if parse_tags:
                        MB = aln.get_tag('UM')
                    else:
                        MB = match.group('MB')

                    if gene_tags:
                        target_name = aln.get_tag('GX').split(',')[0]
                    else:
                        txid = sam_file.getrname(aln.reference_id)
                        if gene_map:
                            if txid in gene_map:
                                target_name = gene_map[txid]
                            else:
                                missing_transcripts.add(txid)
                                continue
                        else:
                            target_name = txid
                    targets_seen.add(target_name)

                    # Scale evidence by number of hits
                    evidence[CB][MB] += weigh_evidence(aln.tags)
                    bare_evidence[CB] += weigh_evidence(aln.tags)
                    kept += 1
                transcripts_processed += 1
                if not transcripts_processed % 1000:
                    logger.info("%d genes processed." % genes_processed)
                    logger.info("%d transcripts processed." % transcripts_processed)
                    logger.info("%d alignments processed." % alignments_processed)

            earray = []
            for cell in cells:
                umis = [1 for _, v in evidence[cell].items() if v >= minevidence]
                earray.append(str(sum(umis)))
            out_handle.write(",".join([gene] + earray) + "\n")
            earray = []
            if umi_matrix:
                for cell in cells:
                    earray.append(str(int(bare_evidence[cell])))
                bare_evidence_handle.write(",".join([gene] + earray) + "\n")

            evidence = collections.defaultdict(lambda: collections.defaultdict(int))
            bare_evidence = collections.defaultdict(int)
            genes_processed += 1

    if umi_matrix:
        bare_evidence_handle.close()

    # fill dataframe with missing values, sort and output
    df = pd.read_csv(out, index_col=0, header=0)
    targets = pd.Series(index=set(transcript_map.keys()))
    targets = targets.sort_index()
    df = df.reindex(targets.index.values, fill_value=0)
    df = df.sort_index()
    df.to_csv(out)

    if umi_matrix:
        df = pd.read_csv(umi_matrix, index_col=0, header=0)
        df = df.reindex(targets.index.values, fill_value=0)
        df = df.sort_index()
        df.to_csv(umi_matrix)
848059
def main(self):
        
        # Find the target files
        self.targets()
        kmer = 15 if self.analysistype == 'GDCS' else 17
        # Use bbduk to bait the FASTQ reads matching the target sequences
        self.bait(maskmiddle='t', k=kmer)
        # If desired, use bbduk to bait the target sequences with the previously baited FASTQ files
        if self.revbait:
            self.reversebait(maskmiddle='t', k=kmer)
        # Run the bowtie2 read mapping module
        self.mapping()
        # Use samtools to index the sorted bam file
        self.indexing()
        # Parse the results
        # self.parsing()
        self.parsebam()
1190268
def blast(args):
    
    from jcvi.apps.align import run_megablast

    p = OptionParser(blast.__doc__)
    p.add_option("-n", type="int", default=2,
            help="Take best N hits [default: %default]")
    opts, args = p.parse_args(args)

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

    allfasta, clonename = args
    fastadir = "fasta"
    infile = op.join(fastadir, clonename + ".fasta")
    if not op.exists(infile):
        entrez([clonename, "--skipcheck", "--outdir=" + fastadir])

    outfile = "{0}.{1}.blast".format(clonename, allfasta.split(".")[0])
    run_megablast(infile=infile, outfile=outfile, db=allfasta, \
            pctid=GoodPct, hitlen=GoodOverlap)

    blasts = [BlastLine(x) for x in open(outfile)]
    besthits = []
    for b in blasts:
        if b.query.count("|") >= 3:
            b.query = b.query.split("|")[3]

        if b.subject.count("|") >= 3:
            b.subject = b.subject.split("|")[3]

        b.query = b.query.rsplit(".", 1)[0]
        b.subject = b.subject.rsplit(".", 1)[0]

        if b.query == b.subject:
            continue

        if b.subject not in besthits:
            besthits.append(b.subject)
        if len(besthits) == opts.n:
            break

    for b in besthits:
        overlap([clonename, b, "--dir=" + fastadir])
1009590
def _sample_template(sample, out_dir):
    """"""
    bam_fn = dd.get_work_bam(sample)
    genome = dd.get_genome_build(sample)
    if genome in supported:
        peaks = sample.get("peaks_files", []).get("main")
        if peaks:
            r_code = ("library(ChIPQC);\n"
                      "sample = ChIPQCsample(\"{bam_fn}\","
                      "\"{peaks}\", "
                      "annotation = \"{genome}\","
                      ");\n"
                      "ChIPQCreport(sample);\n")
            r_code_fn = os.path.join(out_dir, "chipqc.r")
            with open(r_code_fn, 'w') as inh:
                inh.write(r_code.format(**locals()))
            return r_code_fn
509737
def blastp(seqs, blast_db="nr", e_value="1e-20", max_hits=200,
           working_dir=tempfile.gettempdir(), blast_mat_root=None,
           extra_params={}):
    

    # set up params to use with blastp
    params = {
        # matrix
        "-M":"BLOSUM62",

        # max procs
        "-a":"1",

        # expectation
        "-e":e_value,

        # max seqs to show
        "-b":max_hits,

        # max one line descriptions
        "-v":max_hits,

        # program
        "-p":"blastp"
    }
    params.update(extra_params)

    # blast
    blast_res =  blast_seqs(seqs,
        Blastall,
        blast_mat_root=blast_mat_root,
        blast_db=blast_db,
        params=params,
        add_seq_names=False,
        WorkingDir=working_dir
        )

    # get prot id map
    if blast_res['StdOut']:
        lines = [x for x in blast_res['StdOut']]
        return BlastResult(lines)

    return None
62850
def cluster_seqs(seqs,
                 neighbor_join=False,
                 params={},
                 add_seq_names=True,
                 WorkingDir=tempfile.gettempdir(),
                 SuppressStderr=None,
                 SuppressStdout=None,
                 max_chars=1000000,
                 max_hours=1.0,
                 constructor=PhyloNode,
                 clean_up=True
                 ):
    
    num_seqs = len(seqs)
    if num_seqs < 2:
        raise ValueError, "Muscle requres 2 or more sequences to cluster."


    num_chars = sum(map(len, seqs))
    if num_chars > max_chars:
        params["-maxiters"] = 2
        params["-diags1"] = True
        params["-sv"] = True
        #params["-distance1"] = "kmer6_6"
        #params["-distance1"] = "kmer20_3"
        #params["-distance1"] = "kbit20_3"
        print "lots of chars, using fast align", num_chars


    params["-maxhours"] = max_hours
    #params["-maxiters"] = 10

    #cluster_type = "upgmb"
    #if neighbor_join:
    #    cluster_type = "neighborjoining"

    params["-clusteronly"] = True
    params["-tree1"] = get_tmp_filename(WorkingDir)

    muscle_res = muscle_seqs(seqs,
                 params=params,
                 add_seq_names=add_seq_names,
                 WorkingDir=WorkingDir,
                 SuppressStderr=SuppressStderr,
                 SuppressStdout=SuppressStdout)

    tree = DndParser(muscle_res["Tree1Out"], constructor=constructor)

    if clean_up:
        muscle_res.cleanUp()
    return tree
133334
def extract_pairs(fastqfile, p1fw, p2fw, fragsfw, p, suffix=False):
    
    fp = open(fastqfile)
    currentID = 0
    npairs = nfrags = 0
    for x, lib in izip(p.r1, p.libs):
        while currentID != x:
            fragsfw.writelines(islice(fp, 4))  # Exhaust the iterator
            currentID += 1
            nfrags += 1
        a = list(islice(fp, 4))
        b = list(islice(fp, 4))
        if suffix:
            name = a[0].rstrip()
            a[0] = name + "/1\n"
            b[0] = name + "/2\n"
        else:
            b[0] = a[0]  # Keep same read ID for pairs

        p1fw[lib].writelines(a)
        p2fw[lib].writelines(b)
        currentID += 2
        npairs += 2

    # Write the remaining single reads
    while True:
        contents = list(islice(fp, 4))
        if not contents:
            break
        fragsfw.writelines(contents)
        nfrags += 1

    logging.debug("A total of {0} paired reads written to `{1}`.".\
                  format(npairs, ",".join(x.name for x in p1fw + p2fw)))
    logging.debug("A total of {0} single reads written to `{1}`.".\
                  format(nfrags, fragsfw.name))

    # Validate the numbers
    expected_pairs = 2 * p.npairs
    expected_frags = p.nreads - 2 * p.npairs
    assert npairs == expected_pairs, "Expect {0} paired reads, got {1} instead".\
              format(expected_pairs, npairs)
    assert nfrags == expected_frags, "Expect {0} single reads, got {1} instead".\
              format(expected_frags, nfrags)
947816
def rnaseq_prep_samples(config, run_info_yaml, parallel, dirs, samples):
    
    pipeline = dd.get_in_samples(samples, dd.get_analysis)
    trim_reads_set = any([tz.get_in(["algorithm", "trim_reads"], d) for d in dd.sample_data_iterator(samples)])
    resources = ["picard"]
    needs_trimming = (_is_smallrnaseq(pipeline) or trim_reads_set)
    if needs_trimming:
        resources.append("atropos")
    with prun.start(_wres(parallel, resources),
                    samples, config, dirs, "trimming",
                    max_multicore=1 if not needs_trimming else None) as run_parallel:
        with profile.report("organize samples", dirs):
            samples = run_parallel("organize_samples", [[dirs, config, run_info_yaml,
                                                            [x[0]["description"] for x in samples]]])
            samples = run_parallel("prepare_sample", samples)
        if needs_trimming:
            with profile.report("adapter trimming", dirs):
                if _is_smallrnaseq(pipeline):
                    samples = run_parallel("trim_srna_sample", samples)
                else:
                    samples = run_parallel("trim_sample", samples)
    return samples
373782
def sspro8_results(self):
        
        return ssbio.protein.sequence.utils.fasta.load_fasta_file_as_dict_of_seqs(self.out_sspro8)
1367498
def _flatten_samples(samples, base_file, get_retriever):
    
    flat_data = []
    for data in samples:
        data["reference"] = _indexes_to_secondary_files(data["reference"], data["genome_build"])
        cur_flat = {}
        for key_path in [["analysis"], ["description"], ["rgnames"], ["config", "algorithm"],
                         ["metadata"], ["genome_build"], ["resources"],
                         ["files"], ["reference"], ["genome_resources"], ["vrn_file"]]:
            cur_key = "__".join(key_path)
            for flat_key, flat_val in _to_cwldata(cur_key, tz.get_in(key_path, data), get_retriever):
                cur_flat[flat_key] = flat_val
        flat_data.append(cur_flat)
    out = {}
    for key in sorted(list(set(reduce(operator.add, [list(d.keys()) for d in flat_data])))):
        # Periods in keys cause issues with WDL and some CWL implementations
        clean_key = key.replace(".", "_")
        out[clean_key] = []
        for cur_flat in flat_data:
            out[clean_key].append(cur_flat.get(key))
    # special case for back-compatibility with fasta specifications -- yuck
    if "reference__fasta__base" not in out and "reference__fasta" in out:
        out["reference__fasta__base"] = out["reference__fasta"]
        del out["reference__fasta"]
    return _samplejson_to_inputs(out), out
486473
def tiered_alignment(in_bam, tier_num, multi_mappers, extra_args,
                     genome_build, pair_stats,
                     work_dir, dirs, config):
    
    nomap_fq1, nomap_fq2 = select_unaligned_read_pairs(in_bam, "tier{}".format(tier_num),
                                                       work_dir, config)
    if nomap_fq1 is not None:
        base_name = "{}-tier{}out".format(os.path.splitext(os.path.basename(in_bam))[0],
                                          tier_num)
        config = copy.deepcopy(config)
        dirs = copy.deepcopy(dirs)
        config["algorithm"]["bam_sort"] = "queryname"
        config["algorithm"]["multiple_mappers"] = multi_mappers
        config["algorithm"]["extra_align_args"] = ["-i", int(pair_stats["mean"]),
                                               int(pair_stats["std"])] + extra_args
        out_bam, ref_file = align_to_sort_bam(nomap_fq1, nomap_fq2,
                                              lane.rg_names(base_name, base_name, config),
                                              genome_build, "novoalign",
                                              dirs, config,
                                              dir_ext=os.path.join("hydra", os.path.split(nomap_fq1)[0]))
        return out_bam
    else:
        return None
655228
def runblast(self, assembly, allele, sample):
        
        genome = os.path.split(assembly)[1].split('.')[0]
        # Run the BioPython BLASTn module with the genome as query, fasta(target gene) as db.
        # Do not re-perform the BLAST search each time
        make_path(sample[self.analysistype].reportdir)
        try:
            report = glob('{}{}*rawresults*'.format(sample[self.analysistype].reportdir, genome))[0]
            size = os.path.getsize(report)
            if size == 0:
                os.remove(report)
                report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                          time.strftime("%Y.%m.%d.%H.%M.%S"))
        except IndexError:
            report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                      time.strftime("%Y.%m.%d.%H.%M.%S"))
        db = allele.split('.')[0]
        # BLAST command line call. Note the mildly restrictive evalue, and the high number of alignments.
        # Due to the fact that all the targets are combined into one database, this is to ensure that all potential
        # alignments are reported. Also note the custom outfmt: the doubled quotes are necessary to get it work
        blastn = NcbiblastnCommandline(query=assembly, db=db, evalue='1E-20', num_alignments=1000000,
                                       num_threads=12,
                                       outfmt="'6 qseqid sseqid positive mismatch gaps "
                                              "evalue bitscore slen length qstart qend qseq sstart send'",
                                       out=report)
        # Save the blast command in the metadata
        sample[self.analysistype].blastcommand = str(blastn)
        sample[self.analysistype].blastreport = report
        if not os.path.isfile(report):
            # Run BLAST
            blastn()
        # Run the blast parsing module
        self.blastparser(report, sample)
812193
def write_16S_rRNA_fasta(self, org_list):

        '''
        
        '''

        fasta = []
        for org in org_list:

            handle = open(org.genome_path, "rU")


            seq_record = SeqIO.read(handle, "genbank")

            for feat in seq_record.features:
                if feat.type == 'rRNA':
                    if '16S ribosomal' in feat.qualifiers['product'][0]:
                        start = feat.location.start.position
                        end = feat.location.end.position

                        if ((end - start) > 1400) & ((end - start) < 1700) :

                            print 'rRNA sequence extracted for', org.accession

                            fasta.append('>' + org.accession +
                                         '\n' + 
                                         str(feat.extract(seq_record.seq)) +
                                          '\n')

                            org.rRNA16S_sequence = str(feat.extract(seq_record.seq))
                            
                            break

                    
        faastring = "".join(fasta)
        filename = '16S-rRNA.fasta'
        write_fasta = open(filename, 'w')
        write_fasta.write(faastring)
        write_fasta.close()
751401
def _run_genotype_gvcfs_genomicsdb(genomics_db, region, out_file, data):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            broad_runner = broad.runner_from_config(data["config"])
            params = ["-T", "GenotypeGVCFs",
                      "--variant", "gendb://%s" % genomics_db,
                      "-R", dd.get_ref_file(data),
                      "--output", tx_out_file,
                      "-L", bamprep.region_to_gatk(region)]
            params += ["-ploidy", str(ploidy.get_ploidy([data], region))]
            # Avoid slow genotyping runtimes with improved quality score calculation in GATK4
            # https://gatkforums.broadinstitute.org/gatk/discussion/11471/performance-troubleshooting-tips-for-genotypegvcfs/p1
            params += ["--use-new-qual-calculator"]
            resources = config_utils.get_resources("gatk", data["config"])
            params += [str(x) for x in resources.get("options", [])]
            cores = dd.get_cores(data)
            memscale = {"magnitude": 0.9 * cores, "direction": "increase"} if cores > 1 else None
            broad_runner.run_gatk(params, memscale=memscale)
    return vcfutils.bgzip_and_index(out_file, data["config"])
1285542
def run(n = 10)
      raise NotEnoughHitsError if hits.length < opt[:min_blast_hits]
      raise unless prediction.is_a?(Query) && !prediction.raw_sequence.nil? &&
                   hits[0].is_a?(Query)

      start = Time.new
      # get the first n hits
      n_hits = [n - 1, @hits.length].min
      less_hits = @hits[0..n_hits]
      # get raw sequences for less_hits
      less_hits.delete_if do |hit|
        if hit.raw_sequence.nil?
          hit.raw_sequence = FetchRawSequences.run(hit.identifier,
                                                   hit.accession_no)
        end
        hit.raw_sequence.nil? ? true : false
      end

      raise NoInternetError if less_hits.length.zero?

      averages = []

      less_hits.each do |hit|
        coverage = Array.new(hit.length_protein, 0)
        # each residue of the protein should be evluated once only
        ranges_prediction = []

        hit.hsp_list.each do |hsp|
          # align subsequences from the hit and prediction that match
          if !hsp.hit_alignment.nil? && !hsp.query_alignment.nil?
            hit_alignment   = hsp.hit_alignment
            query_alignment = hsp.query_alignment
          else
            align = find_local_alignment(hit, prediction, hsp)
            hit_alignment   = align[0]
            query_alignment = align[1]
          end

          coverage = check_multiple_coverage(hit_alignment, query_alignment,
                                             hsp, coverage, ranges_prediction)

          ranges_prediction << (hsp.match_query_from..hsp.match_query_to)
        end
        overlap = coverage.reject(&:zero?)
        if overlap != []
          averages.push((overlap.inject(:+) / (overlap.length + 0.0)).round(2))
        end
      end

      # if all hsps match only one time
      if averages.reject { |x| x == 1 } == []
        @validation_report = DuplicationValidationOutput.new(@short_header,
                                                             @header,
                                                             @description, 1,
                                                             averages)
        @validation_report.run_time = Time.now - start
        return @validation_report
      end

      pval = wilcox_test(averages)

      @validation_report = DuplicationValidationOutput.new(@short_header,
                                                           @header,
                                                           @description, pval,
                                                           averages)
      @run_time = Time.now - start
      @validation_report
    rescue NotEnoughHitsError
      @validation_report = ValidationReport.new('Not enough evidence', :warning,
                                                @short_header, @header,
                                                @description)
    rescue NoMafftInstallationError
      @validation_report = ValidationReport.new('Mafft error', :error,
                                                @short_header, @header,
                                                @description)
      @validation_report.errors.push NoMafftInstallationError
    rescue NoInternetError
      @validation_report = ValidationReport.new('Internet error', :error,
                                                @short_header, @header,
                                                @description)
      @validation_report.errors.push NoInternetError
    rescue StandardError
      @validation_report = ValidationReport.new('Unexpected error', :error,
                                                @short_header, @header,
                                                @description)
      @validation_report.errors.push 'Unexpected Error'
    end