Extract metrics from a fastq file.

    Return average quality and read length
344149
def extract_from_fastq(fq):
    
    for rec in SeqIO.parse(fq, "fastq"):
        yield nanomath.ave_qual(rec.letter_annotations["phred_quality"]), len(rec)
411608
def iter(self):
        
        for _file in self._files:
            with asHandle(_file) as fp:
                # Use FastqGeneralIterator because it provides access to
                # the unconverted quality string (i.e., it doesn't try to
                # figure out the numeric quality values, which we don't
                # care about at this point).
                for sequenceId, sequence, quality in FastqGeneralIterator(fp):
                    yield self.readClass(sequenceId, sequence, quality)
370111
def filer(self):
        
        # Get a list of all the .fastq files in the sequence path
        fastq = sorted(glob(os.path.join(self.sequencepath, '*fastq*')))
        # Initialise a list to store .fasta formatted files
        fasta = list()
        # Get all the files in the sequence path
        for sequencefile in glob(os.path.join(self.sequencepath, '*')):
            # Split the extension from the file
            extension = os.path.splitext(sequencefile)[1]
            # If the extension is one of the acceptable FASTA extensions, add it to the list of FASTA files
            if extension in self.extensions:
                fasta.append(sequencefile)
        # Use filer from accessoryFunctions to find the base file name of each .fastq file/pair
        fastqnames = filer(fastq)
        # Get all the .fastq files into metadata objects
        for fastqname in sorted(fastqnames):
            # Create the metadata object
            sample = MetadataObject()
            sample.name = os.path.basename(fastqname)
            # Create an attribute for the current analysis
            setattr(sample, self.analysistype, GenObject())
            # Set the destination folder
            sample[self.analysistype].outputdir = os.path.join(self.path, fastqname)
            # Make the destination folder
            make_path(sample[self.analysistype].outputdir)
            # Get the fastq files specific to the fastqname
            specificfastq = glob(os.path.join(self.sequencepath, '{fq}*.fastq*'.format(fq=fastqname)))
            # Set the file type for the downstream analysis
            sample[self.analysistype].filetype = 'fastq'
            # Link the files to the output folder
            for fastq in specificfastq:
                relative_symlink(src_file=fastq,
                                 output_dir=sample[self.analysistype].outputdir)
            # Initialise the general and run categories
            sample.general = GenObject()
            # Populate the .fastqfiles category of :self.metadata
            sample.general.fastqfiles = [fastq for fastq in
                                         glob(os.path.join(sample[self.analysistype].outputdir, '{fq}*.fastq*'
                                                           .format(fq=fastqname))) if 'trimmed' not in fastq]
            # Populate certain attributes in order to be compatible with other local software
            sample.general.bestassemblyfile = sample.general.fastqfiles
            sample.general.trimmedcorrectedfastqfiles = sample.general.fastqfiles
            sample.general.outputdirectory = sample[self.analysistype].outputdir
            sample[self.analysistype].baitedfastq = os.path.join(
                sample[self.analysistype].outputdir,
                '{at}_targetMatches.fastq.gz'.format(at=self.analysistype))
            # Append the metadata to the list of samples
            self.metadata.append(sample)
        # Add any .fasta formatted files to the metadata object
        for fastafile in sorted(fasta):
            sample = MetadataObject()
            sample.name = os.path.splitext(os.path.split(fastafile)[-1])[0]
            setattr(sample, self.analysistype, GenObject())
            # Set the destination folder
            sample[self.analysistype].outputdir = os.path.join(self.path, 'detailed_reports', sample.name)
            # Make the destination folder
            make_path(sample[self.analysistype].outputdir)
            # Set the file type for the downstream analysis
            sample[self.analysistype].filetype = 'fasta'
            # Link the .fasta files to :self.path/:filename
            relative_symlink(src_file=fastafile,
                             output_dir=sample[self.analysistype].outputdir)
            # Initialise the general and run categories
            sample.general = GenObject()
            # Populate the .fastqfiles category of :self.metadata
            sample.general.fastqfiles = [fastq for fastq in glob(os.path.join(sample[self.analysistype].outputdir,
                                                                              '{sn}*{ext}*'
                                                                              .format(
                                                                                  sn=sample.name,
                                                                                  ext=os.path.splitext(fastafile)[1])))]
            sample.general.bestassemblyfile = sample.general.fastqfiles
            sample.general.trimmedcorrectedfastqfiles = sample.general.fastqfiles
            sample.general.outputdirectory = sample[self.analysistype].outputdir
            sample[self.analysistype].baitedfastq = os.path.join(
                sample[self.analysistype].outputdir,
                '{at}_targetMatches.fastq.gz'.format(at=self.analysistype))
            sample[self.analysistype].assemblyfile = fastafile
            # Append the metadata to the list of samples
            self.metadata.append(sample)
1058553
def process_fastq_minimal(fastq, **kwargs):
    """"""
    infastq = handle_compressed_input(fastq)
    try:
        df = pd.DataFrame(
            data=[rec for rec in fq_minimal(infastq) if rec],
            columns=["timestamp", "lengths"]
        )
    except IndexError:
        logging.error("Fatal: Incorrect file structure for fastq_minimal")
        sys.exit("Error: file does not match expected structure for fastq_minimal")
    return ut.reduce_memory_usage(df)
1105773
def read_file(self, fasta_path):
        
        fasta_dictionary = {}
        for (identifier, sequence) in self.iterate_over_file(fasta_path):
            fasta_dictionary[identifier] = sequence
        return fasta_dictionary
1294184
def process_fasta(fasta, **kwargs):
    """"""
    logging.info("Nanoget: Starting to collect statistics from a fasta file.")
    inputfasta = handle_compressed_input(fasta, file_type="fasta")
    return ut.reduce_memory_usage(pd.DataFrame(
        data=[len(rec) for rec in SeqIO.parse(inputfasta, "fasta")],
        columns=["lengths"]
    ).dropna())
423875
def fastq(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(fastq.__doc__)
    p.add_option("--outtie", dest="outtie", default=False, action="store_true",
            help="Are these outie reads? [default: %default]")
    p.set_phred()
    p.set_size()

    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(p.print_help())

    fastqfiles = [get_abs_path(x) for x in args]
    size = opts.size
    outtie = opts.outtie
    if size > 1000 and (not outtie):
        logging.debug("[warn] long insert size {0} but not outtie".format(size))

    mated = (size != 0)
    libname = op.basename(args[0]).split(".")[0]
    libname = libname.replace("_1_sequence", "")

    frgfile = libname + ".frg"
    mean, sv = get_mean_sv(opts.size)

    cmd = "fastqToCA"
    cmd += " -libraryname {0} ".format(libname)
    fastqs = " ".join("-reads {0}".format(x) for x in fastqfiles)
    if mated:
        assert len(args) in (1, 2), "you need one or two fastq files for mated library"
        fastqs = "-mates {0}".format(",".join(fastqfiles))
        cmd += "-insertsize {0} {1} ".format(mean, sv)
    cmd += fastqs

    offset = int(opts.phred) if opts.phred else guessoffset([fastqfiles[0]])
    illumina = (offset == 64)
    if illumina:
        cmd += " -type illumina"
    if outtie:
        cmd += " -outtie"

    sh(cmd, outfile=frgfile)
546452
def get_fasta_records(fasta, transcripts):
    
    records = {}
    f = pysam.FastaFile(fasta)
    for transcript in transcripts:
        try:
            sequence, length = f.fetch(transcript), f.get_reference_length(transcript)
        except KeyError:
            msg = 'Transcript "{}" does not exist in transcriptome FASTA file'.format(transcript)
            log.error(msg)
            raise ArgumentError(msg)
        records[transcript] = {'sequence': sequence, 'length': length}
    f.close()
    return records
538106
def process_file_metrics(context, file_processors):
    """"""
    file_metrics = OrderedDict()

    # TODO make available the includes and excludes feature
    gitignore = []
    if os.path.isfile('.gitignore'):
        with open('.gitignore', 'r') as ifile:
            gitignore = ifile.read().splitlines()

    in_files = glob_files(context['root_dir'], context['in_file_names'], gitignore=gitignore)
    # main loop
    for in_file, key in in_files:
        # print 'file %i: %s' % (i, in_file)
        try:
            with open(in_file, 'rb') as ifile:
                code = ifile.read()
            # lookup lexicographical scanner to use for this run
            try:
                lex = guess_lexer_for_filename(in_file, code, encoding='guess')
                # encoding is 'guess', chardet', 'utf-8'
            except:
                pass
            else:
                token_list = lex.get_tokens(code)  # parse code

                file_metrics[key] = OrderedDict()
                file_metrics[key].update(compute_file_metrics(file_processors, lex.name, key, token_list))
                file_metrics[key]['language'] = lex.name

        except IOError as e:
            sys.stderr.writelines(str(e) + " -- Skipping input file.\n\n")

    return file_metrics
900989
def _gzip_fastq(in_file, out_dir=None):
    
    if fastq.is_fastq(in_file) and not objectstore.is_remote(in_file):
        if utils.is_bzipped(in_file):
            return _bzip_gzip(in_file, out_dir)
        elif not utils.is_gzipped(in_file):
            if out_dir:
                gzipped_file = os.path.join(out_dir, os.path.basename(in_file) + ".gz")
            else:
                gzipped_file = in_file + ".gz"
            if file_exists(gzipped_file):
                return gzipped_file
            message = "gzipping {in_file} to {gzipped_file}.".format(
                in_file=in_file, gzipped_file=gzipped_file)
            with file_transaction(gzipped_file) as tx_gzipped_file:
                do.run("gzip -c {in_file} > {tx_gzipped_file}".format(**locals()),
                       message)
            return gzipped_file
    return in_file
993266
def write_fastq(filename):
    
    if filename:
        if filename.endswith('gz'):
            filename_fh = gzip.open(filename, mode='wb')
        else:
            filename_fh = open(filename, mode='w')
    else:
        filename_fh = None
    return filename_fh
827390
def to_fastq_apipe_cl(sdf_file, start=None, end=None):
    
    cmd = ["rtg", "sdf2fastq", "--no-gzip", "-o", "-"]
    if start is not None:
        cmd += ["--start-id=%s" % start]
    if end is not None:
        cmd += ["--end-id=%s" % end]
    if is_paired(sdf_file):
        out = []
        for ext in ["left", "right"]:
            out.append("<(%s)" % _rtg_cmd(cmd + ["-i", os.path.join(sdf_file, ext)]))
        return out
    else:
        cmd += ["-i", sdf_file]
        return ["<(%s)" % _rtg_cmd(cmd), None]
85645
def load_fasta_file_as_dict_of_seqs(filename):
    

    results = {}
    records = load_fasta_file(filename)
    for r in records:
        results[r.id] = str(r.seq)

    return results
596155
def _get_fastq_files(ldetail, read, fastq_dir):
    
    return glob.glob(os.path.join(fastq_dir, "Project_%s" % ldetail["project_name"],
                                  "Sample_%s" % ldetail["name"],
                                  "%s_*_%s_*.fastq.gz" % (ldetail["name"], read)))
634899
def _merge_fastqc(samples):
    
    fastqc_list = collections.defaultdict(list)
    seen = set()
    for data in samples:
        name = dd.get_sample_name(data)
        if name in seen:
            continue
        seen.add(name)
        fns = glob.glob(os.path.join(dd.get_work_dir(data), "qc", dd.get_sample_name(data), "fastqc") + "/*")
        for fn in fns:
            if fn.endswith("tsv"):
                metric = os.path.basename(fn)
                fastqc_list[metric].append([name, fn])

    for metric in fastqc_list:
        dt_by_sample = []
        for fn in fastqc_list[metric]:
            dt = pd.read_csv(fn[1], sep="\t")
            dt['sample'] = fn[0]
            dt_by_sample.append(dt)
        dt = utils.rbind(dt_by_sample)
        dt.to_csv(metric, sep="\t", index=False, mode ='w')
    return samples
307631
def validate_fastq(self):
        
        logging.info('Validating FASTQ files')
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # Tiny files can pass the validation tests - ensure that they don't
                size = os.path.getsize(sample.general.fastqfiles[0])
                if size >= 1000000:
                    # Try to run reformat.sh on the reads - on any errors try to run repair.sh
                    try:
                        out, err, cmd = bbtools.validate_reads(forward_in=sample.general.fastqfiles[0],
                                                               returncmd=True)
                        write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr,
                                         None, None)
                    except CalledProcessError:
                        # Set the file names for the reformatted and repaired files
                        outputfile1 = os.path.join(sample.general.outputdirectory, '{}_reformatted_R1.fastq.gz'
                                                   .format(sample.name))
                        repair_file1 = os.path.join(sample.general.outputdirectory, '{}_repaired_R1.fastq.gz'
                                                    .format(sample.name))
                        if len(sample.general.fastqfiles) == 2:
                            outputfile2 = os.path.join(sample.general.outputdirectory, '{}_reformatted_R2.fastq.gz'
                                                       .format(sample.name))
                            repair_file2 = os.path.join(sample.general.outputdirectory, '{}_repaired_R2.fastq.gz'
                                                        .format(sample.name))
                        else:
                            outputfile2 = str()
                            repair_file2 = str()
                        # Use reformat.sh to repair the reads. If this fails, discard the sample from the analyses
                        try:
                            logging.warning('Errors detected in FASTQ files for sample {sample}. '
                                            'Please check the following files for details {log} {logout} {logerr}. '
                                            'The pipeline will use reformat.sh to attempt to repair issues'
                                            .format(sample=sample.name,
                                                    log=self.logfile,
                                                    logout=sample.general.logout,
                                                    logerr=sample.general.logerr))
                            if not os.path.isfile(outputfile1):
                                # Run reformat.sh
                                out, err, cmd = bbtools.reformat_reads(forward_in=sample.general.fastqfiles[0],
                                                                       forward_out=outputfile1,
                                                                       returncmd=True)
                                write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr,
                                                 None, None)
                                # Run repair.sh (if necessary)
                                if outputfile2:
                                    out, err, cmd = bbtools.repair_reads(forward_in=outputfile1,
                                                                         forward_out=repair_file1,
                                                                         returncmd=True)
                                    write_to_logfile(out, err, self.logfile, sample.general.logout,
                                                     sample.general.logerr, None, None)
                            # Ensure that the output file(s) exist before declaring this a success
                            if os.path.isfile(outputfile1):
                                # Update the fastqfiles attribute to point to the repaired files
                                sample.general.fastqfiles = [repair_file1, repair_file2] if repair_file2 \
                                    else [outputfile1]
                        except CalledProcessError:
                            # The file(s) can be created even if there is STDERR from reformat.sh
                            if os.path.isfile(outputfile1) and outputfile2:
                                try:
                                    out, err, cmd = bbtools.repair_reads(forward_in=outputfile1,
                                                                         forward_out=repair_file1,
                                                                         returncmd=True)
                                    write_to_logfile(out, err, self.logfile, sample.general.logout,
                                                     sample.general.logerr, None, None)
                                    # Update the fastqfiles attribute to point to the repaired files
                                    sample.general.fastqfiles = [repair_file1, repair_file2] if repair_file2 else \
                                        [repair_file1]
                                except CalledProcessError:
                                    # Write in the logs that there was an error detected in the FASTQ files
                                    write_to_logfile('An error was detected in the FASTQ files for sample {}. '
                                                     'These files will not be processed further'.format(sample.name),
                                                     'An error was detected in the FASTQ files for sample {}. '
                                                     'These files will not be processed further'.format(sample.name),
                                                     self.logfile,
                                                     sample.general.logout,
                                                     sample.general.logerr,
                                                     None,
                                                     None)
                                    # Update metadata objects with error
                                    self.error(sample, 'fastq_error')
                            else:
                                # Write in the logs that there was an error detected in the FASTQ files
                                write_to_logfile('An error was detected in the FASTQ files for sample {}. '
                                                 'These files will not be processed further'.format(sample.name),
                                                 'An error was detected in the FASTQ files for sample {}. '
                                                 'These files will not be processed further'.format(sample.name),
                                                 self.logfile,
                                                 sample.general.logout,
                                                 sample.general.logerr,
                                                 None,
                                                 None)

                                # Update metadata objects with error
                                self.error(sample, 'fastq_error')
                else:
                    # Update metadata objects with error
                    self.error(sample, 'files_too_small')
        # Print the metadata to file
        metadataprinter.MetadataPrinter(self)
270120
def extract_length(self, lower_bound=None, upper_bound=None, new_path=None):
        """"""
        # Temporary path #
        if new_path is None: fraction = self.__class__(new_temp_path())
        elif isinstance(new_path, FASTA): fraction = new_path
        else:                fraction = self.__class__(new_path)
        # Generator #
        if lower_bound is None: lower_bound = 0
        if upper_bound is None: upper_bound = sys.maxint
        def fraction_iterator():
            for read in self:
                if lower_bound <= len(read) <= upper_bound:
                    yield read
        # Do it #
        fraction.write(fraction_iterator())
        fraction.close()
        return fraction