Convert BAM file to fastq.
499756
def picard_bam_to_fastq(picard, in_bam, fastq_one, fastq_two=None):
    
    if not file_exists(fastq_one):
        with tx_tmpdir(picard._config) as tmp_dir:
            with file_transaction(picard._config, fastq_one) as tx_out1:
                opts = [("INPUT", in_bam),
                        ("FASTQ", tx_out1),
                        ("TMP_DIR", tmp_dir)]
                if fastq_two is not None:
                    opts += [("SECOND_END_FASTQ", fastq_two)]
                picard.run("SamToFastq", opts)
    return (fastq_one, fastq_two)
998383
def trim_fastq(credentials, instance_config, instance_name,
               script_dir, input_file, output_file,
               trim_crop, trim_headcrop=0, self_destruct=True,
               **kwargs):
    """"""
    template = _TEMPLATE_ENV.get_template('trim_fastq.sh')
    startup_script = template.render(
        script_dir=script_dir,
        input_file=input_file,
        output_file=output_file,
        trim_crop=trim_crop,
        trim_headcrop=trim_headcrop,
        self_destruct=self_destruct)

    if len(startup_script) > 32768:
        raise ValueError('Startup script larger than 32,768 bytes!')

    #print(startup_script)

    instance_config.create_instance(
        credentials, instance_name, startup_script=startup_script, **kwargs)
365232
def fastqmover(self):
        """"""
        # Create the sequence path if necessary
        make_path(self.sequencepath)
        # Iterate through all the sample names
        for sample in self.metadata.samples:
            # Make directory variables
            outputdir = os.path.join(self.sequencepath, sample.name)
            # Demultiplexed files will be present in the project name subfolder within the fastq destination folder
            if self.demultiplex:
                glob_dir = self.projectpath
            # Undemultiplexed reads (Undetermined_S0_R1_001.fastq.gz) are present in the fastq destination folder
            else:
                glob_dir = self.fastqdestination
            # Sometimes the files are put in self.fastqdestination rather than self.projectpath? Matt has this problem
            # and I don't understands what's going on
            if not os.path.isdir(glob_dir):
                glob_dir = os.path.dirname(glob_dir)
            # Glob all the .gz files in the subfolders - projectpath/Sample_:sample.name/*.gz
            for fastq in sorted(glob(os.path.join(glob_dir, '*.gz'))):
                fastqname = os.path.basename(fastq)
                # Set the name of the destination file
                outputfile = os.path.join(self.sequencepath, fastqname)
                # Link the file if it doesn't already exist
                if not os.path.isfile(outputfile):
                    relative_symlink(src_file=fastq,
                                     output_dir=self.sequencepath)
            # Repopulate .strainfastqfiles with the freshly-linked/copied files
            if self.demultiplex:
                fastqfiles = glob(os.path.join(self.sequencepath, '{}*.fastq*'.format(sample.name)))
                fastqfiles = sorted([fastq for fastq in fastqfiles if 'trimmed' not in os.path.basename(fastq)])
            # Undemultiplexed files will not have the sample name in the file name
            else:
                fastqfiles = sorted(glob(os.path.join(glob_dir, '*.gz')))
            # Populate the metadata object with the name/path of the fastq files
            sample.general.fastqfiles = fastqfiles
            # Save the outputdir to the metadata object
            sample.run.outputdirectory = outputdir
            sample.general.outputdirectory = outputdir
            sample.general.bestassemblyfile = True
            sample.general.trimmedcorrectedfastqfiles = sorted(sample.general.fastqfiles)
            sample.general.logout = os.path.join(sample.general.outputdirectory, 'logout')
            sample.general.logerr = os.path.join(sample.general.outputdirectory, 'logerr')
            sample.commands = GenObject()
251426
def _cram_to_fastq_region(cram_file, work_dir, base_name, region, data):
    
    ref_file = tz.get_in(["reference", "fasta", "base"], data)
    resources = config_utils.get_resources("bamtofastq", data["config"])
    cores = tz.get_in(["config", "algorithm", "num_cores"], data, 1)
    max_mem = config_utils.convert_to_bytes(resources.get("memory", "1G")) * cores
    rext = "-%s" % region.replace(":", "_").replace("-", "_") if region else "full"
    out_s, out_p1, out_p2, out_o1, out_o2 = [os.path.join(work_dir, "%s%s-%s.fq.gz" %
                                                          (base_name, rext, fext))
                                             for fext in ["s1", "p1", "p2", "o1", "o2"]]
    if not utils.file_exists(out_p1):
        with file_transaction(data, out_s, out_p1, out_p2, out_o1, out_o2) as \
             (tx_out_s, tx_out_p1, tx_out_p2, tx_out_o1, tx_out_o2):
            cram_file = objectstore.cl_input(cram_file)
            sortprefix = "%s-sort" % utils.splitext_plus(tx_out_s)[0]
            cmd = ("bamtofastq filename={cram_file} inputformat=cram T={sortprefix} "
                   "gz=1 collate=1 colsbs={max_mem} exclude=SECONDARY,SUPPLEMENTARY "
                   "F={tx_out_p1} F2={tx_out_p2} S={tx_out_s} O={tx_out_o1} O2={tx_out_o2} "
                   "reference={ref_file}")
            if region:
                cmd += " ranges='{region}'"
            do.run(cmd.format(**locals()), "CRAM to fastq %s" % region if region else "")
    return [[out_p1, out_p2, out_s]]
1172788
def _prepare_bam_file(bam_file, tmp_dir, config):
    
    sort_mode = _get_sort_order(bam_file, config)
    if sort_mode != "queryname":
        bam_file = sort(bam_file, config, "queryname")
    return bam_file
340860
def _pipeline_needs_fastq(config, data):
    
    aligner = config["algorithm"].get("aligner")
    support_bam = aligner in alignment.metadata.get("support_bam", [])
    return aligner and not support_bam
242949
def to_fastq_str(self):
    
    return "@" + self.name + "\n" + self.sequenceData +\
           "\n" + "+" + self.name + "\n" + self.seq_qual
160448
def to_bam(in_file, out_file, data):
    
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = ["samtools", "view", "-O", "BAM", "-o", tx_out_file, in_file]
            do.run(cmd, "Convert CRAM to BAM")
    bam.index(out_file, data["config"])
    return out_file
574342
def _convert_fastq(srafn, outdir, single=False):
    ""
    cmd = "fastq-dump --split-files --gzip {srafn}"
    cmd = "%s %s" % (utils.local_path_export(), cmd)
    sraid = os.path.basename(utils.splitext_plus(srafn)[0])
    if not srafn:
        return None
    if not single:
        out_file = [os.path.join(outdir, "%s_1.fastq.gz" % sraid),
                    os.path.join(outdir, "%s_2.fastq.gz" % sraid)]
        if not utils.file_exists(out_file[0]):
            with utils.chdir(outdir):
                do.run(cmd.format(**locals()), "Covert to fastq %s" % sraid)
        if not utils.file_exists(out_file[0]):
            raise IOError("SRA %s didn't convert, something happened." % srafn)
        return [out for out in out_file if utils.file_exists(out)]
    else:
        raise ValueError("Not supported single-end sra samples for now.")
314623
def prep_fastq_inputs(in_files, data):
    
    if len(in_files) == 1 and _is_bam_input(in_files):
        out = _bgzip_from_bam(in_files[0], data["dirs"], data)
    elif len(in_files) == 1 and _is_cram_input(in_files):
        out = _bgzip_from_cram(in_files[0], data["dirs"], data)
    elif len(in_files) in [1, 2] and _ready_gzip_fastq(in_files, data):
        out = _symlink_in_files(in_files, data)
    else:
        if len(in_files) > 2:
            fpairs = fastq.combine_pairs(in_files)
            pair_types = set([len(xs) for xs in fpairs])
            assert len(pair_types) == 1
            fpairs.sort(key=lambda x: os.path.basename(x[0]))
            organized = [[xs[0] for xs in fpairs]]
            if len(fpairs[0]) > 1:
                organized.append([xs[1] for xs in fpairs])
            in_files = organized
        parallel = {"type": "local", "num_jobs": len(in_files),
                    "cores_per_job": max(1, data["config"]["algorithm"]["num_cores"] // len(in_files))}
        inputs = [{"in_file": x, "read_num": i, "dirs": data["dirs"], "config": data["config"],
                   "is_cwl": "cwl_keys" in data,
                   "rgnames": data["rgnames"]}
                  for i, x in enumerate(in_files) if x]
        out = run_multicore(_bgzip_from_fastq_parallel, [[d] for d in inputs], data["config"], parallel)
    return out
319485
def tofasta(args):
    
    p = OptionParser(tofasta.__doc__)
    p.add_option("--prefix", default="gbfasta",
            help="prefix of output files [default: %default]")
    filenames, accessions, idfile, opts, args = preparegb(p, args)
    prefix = opts.prefix

    GenBank(filenames=filenames, accessions=accessions, idfile=idfile).\
        write_fasta(output=prefix, individual=opts.individual)

    if opts.individual:
        logging.debug("Output written dir {0}".format(prefix))
    else:
        logging.debug("Output written to {0}.fasta".format(prefix))
405846
def _read_bam(bam_fn, precursors):
    
    mode = "r" if bam_fn.endswith("sam") else "rb"
    handle = pysam.Samfile(bam_fn, mode)
    reads = defaultdict(realign)
    for line in handle:
        chrom = handle.getrname(line.reference_id)
        # print("%s %s %s %s" % (line.query_name, line.reference_start, line.query_sequence, chrom))
        query_name = line.query_name
        if query_name not in reads:
            reads[query_name].sequence = line.query_sequence
        iso = isomir()
        iso.align = line
        iso.start = line.reference_start
        iso.subs, iso.add = _realign(reads[query_name].sequence, precursors[chrom], line.reference_start)
        reads[query_name].set_precursor(chrom, iso)

    reads = _clean_hits(reads)
    return reads
726467
def convert_cufflinks_mapq(in_bam, out_bam=None):
    
    CUFFLINKSMAPQ = 255
    if not out_bam:
        out_bam = os.path.splitext(in_bam)[0] + "-cufflinks.bam"
    if utils.file_exists(out_bam):
        return out_bam
    maxmapq = estimate_max_mapq(in_bam)
    if maxmapq == CUFFLINKSMAPQ:
        return in_bam
    logger.info("Converting MAPQ scores in %s to be Cufflinks compatible." % in_bam)
    with pysam.Samfile(in_bam, "rb") as in_bam_fh:
        with pysam.Samfile(out_bam, "wb", template=in_bam_fh) as out_bam_fh:
            for read in in_bam_fh:
                if read.mapq == maxmapq and not read.is_unmapped:
                    read.mapq = CUFFLINKSMAPQ
                out_bam_fh.write(read)
    return out_bam
178254
def seqs_from_fastacmd(acc_list, blast_db,is_protein=True):
    """"""
    fasta_cmd_res = fasta_cmd_get_seqs(acc_list, blast_db=blast_db, \
        is_protein=is_protein)
    recs = FastaCmdFinder(fasta_cmd_res['StdOut'])
    result = {}
    for rec in recs:
        try:
            result[rec[0][1:].strip()] = ''.join(map(strip, rec[1:]))
        except IndexError:  #maybe we didn't get a sequence?
            pass
    fasta_cmd_res.cleanUp()
    return result
109282
def trimquality(self):
        """"""
        logging.info("Trimming fastq files")
        # Iterate through strains with fastq files
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # As the metadata can be populated with 'NA' (string) if there are no fastq files, only process if
                # :fastqfiles is a list
                if type(sample.general.fastqfiles) is list:
                    # Check to see if the fastq files exist
                    fastqfiles = sorted(sample.general.fastqfiles)
                    # Define the output directory
                    outputdir = sample.general.outputdirectory
                    # Define the name of the trimmed fastq files
                    cleanforward = os.path.join(outputdir, '{}_R1_trimmed.fastq.gz'.format(sample.name))
                    cleanreverse = os.path.join(outputdir, '{}_R2_trimmed.fastq.gz'.format(sample.name))
                    # Incorporate read length into the minlength parameter - set it to 50 unless one or more of the
                    # reads has a lower calculated length than 50
                    try:
                        lesser_length = min(int(sample.run.forwardlength), int(sample.run.reverselength))
                    except ValueError:
                        lesser_length = int(sample.run.forwardlength)
                    min_len = 50 if lesser_length >= 50 else lesser_length
                    # Initialise a variable to store the number of bases to automatically trim from the beginning of
                    # each read, as these bases tend to have lower quality scores. If trimming the reads will cause
                    trim_left = 0
                    # If, for some reason, only the reverse reads are present, use the appropriate output file name
                    try:
                        if 'R2' in fastqfiles[0]:
                            if not os.path.isfile(cleanreverse):
                                out, \
                                    err, \
                                    bbdukcall = bbtools.bbduk_trim(forward_in=fastqfiles[0],
                                                                   reverse_in=None,
                                                                   forward_out=cleanreverse,
                                                                   trimq=10,
                                                                   minlength=min_len,
                                                                   forcetrimleft=trim_left,
                                                                   returncmd=True)
                            else:
                                bbdukcall = str()
                                out = str()
                                err = str()
                        else:
                            if not os.path.isfile(cleanforward):
                                    out, \
                                        err, \
                                        bbdukcall = bbtools.bbduk_trim(forward_in=fastqfiles[0],
                                                                       forward_out=cleanforward,
                                                                       trimq=10,
                                                                       minlength=min_len,
                                                                       forcetrimleft=trim_left,
                                                                       returncmd=True)
                            else:
                                bbdukcall = str()
                                out = str()
                                err = str()
                    except (IndexError, CalledProcessError):
                        bbdukcall = str()
                        out = str()
                        err = str()
                    # Write the command, stdout, and stderr to the logfile
                    write_to_logfile(bbdukcall, bbdukcall, self.logfile, sample.general.logout, sample.general.logerr,
                                     None, None)
                    write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                    # Add the trimmed fastq files to a list
                    trimmedfastqfiles = sorted(glob(os.path.join(sample.general.outputdirectory, '*trimmed.fastq.gz')))
                    # Populate the metadata if the files exist
                    sample.general.trimmedfastqfiles = trimmedfastqfiles if trimmedfastqfiles else list()
        # Add all the trimmed files to the metadata
        logging.info('Fastq files trimmed')
1337233
def _get_fastq_in_region(region, align_bam, out_base):
    
    out_file = "{0}.fastq".format(out_base)
    if not file_exists(out_file):
        with pysam.Samfile(align_bam, "rb") as in_pysam:
            with file_transaction(out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    contig, start, end = region
                    for read in in_pysam.fetch(contig, int(start), int(end)):
                        seq = Seq.Seq(read.seq)
                        qual = list(read.qual)
                        if read.is_reverse:
                            seq = seq.reverse_complement()
                            qual.reverse()
                        out_handle.write("@{name}\n{seq}\n+\n{qual}\n".format(
                                name=read.qname, seq=str(seq), qual="".join(qual)))
    return out_file
808248
def fastqRead(fileHandleOrFile):
    
    fileHandle = _getFileHandle(fileHandleOrFile)
    line = fileHandle.readline()
    while line != '':
        if line[0] == '@':
            name = line[1:-1]
            seq = fileHandle.readline()[:-1]
            plus = fileHandle.readline()
            if plus[0] != '+':
                raise RuntimeError("Got unexpected line: %s" % plus)
            qualValues = [ ord(i) for i in fileHandle.readline()[:-1] ]
            if len(seq) != len(qualValues):
                logger.critical("Got a mismatch between the number of sequence characters (%s) and number of qual values (%s) for sequence: %s, ignoring returning None" % (len(seq), len(qualValues), name))
                qualValues = None
            else:
                for i in qualValues:
                    if i < 33 or i > 126:
                        raise RuntimeError("Got a qual value out of range %s (range is 33 to 126)" % i)
            for i in seq:
                #For safety and sanity I only allows roman alphabet characters in fasta sequences.
                if not ((i >= 'A' and i <= 'Z') or (i >= 'a' and i <= 'z') or i == '-'):
                    raise RuntimeError("Invalid FASTQ character, ASCII code = \'%d\', found in input sequence %s" % (ord(i), name))
            yield name, seq, qualValues
        line = fileHandle.readline()
    if isinstance(fileHandleOrFile, "".__class__):
        fileHandle.close()