Link FASTQ files for a sample.
357362
def link(self, family: str, sample: str, analysis_type: str, files: List[str]):
        """"""
        root_dir = Path(self.families_dir) / family / analysis_type / sample / 'fastq'
        root_dir.mkdir(parents=True, exist_ok=True)
        for fastq_data in files:
            fastq_path = Path(fastq_data['path'])
            fastq_name = self.name_file(
                lane=fastq_data['lane'],
                flowcell=fastq_data['flowcell'],
                sample=sample,
                read=fastq_data['read'],
                undetermined=fastq_data['undetermined'],
            )
            dest_path = root_dir / fastq_name
            if not dest_path.exists():
                log.info(f"linking: {fastq_path} -> {dest_path}")
                dest_path.symlink_to(fastq_path)
            else:
                log.debug(f"destination path already exists: {dest_path}")
653019
def _convert_fastq(fastq_file, pair_file, rg_name, out_file, config):
    
    out_file = "{0}-fq.mkb".format(os.path.splitext(out_file)[0])
    if not file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            cl = [config_utils.get_program("mosaik", config,
                                           default="MosaikAligner").replace("Aligner", "Build")]
            cl += ["-q", fastq_file,
                   "-out", tx_out_file,
                   "-st", config["algorithm"].get("platform", "illumina").lower()]
            if pair_file:
                cl += ["-q2", pair_file]
            if rg_name:
                cl += ["-id", rg_name]
            env_set = "export MOSAIK_TMP={0}".format(os.path.dirname(tx_out_file))
            subprocess.check_call(env_set + " && " + " ".join(cl), shell=True)
    return out_file
352361
def splitFASTA(params):
    
    assert params['fastaFile'][-1] == 'a', ('You must specify a file in '
                                            'fasta-format that ends in '
                                            '.fasta')

    fileCount = count = seqCount = 0
    outfp = None
    with open(params['fastaFile']) as infp:
        for seq in SeqIO.parse(infp, 'fasta'):
            seqCount += 1
            if count == params['seqsPerJob']:
                outfp.close()
                count = 0
            if count == 0:
                outfp = open('%d.fasta' % fileCount, 'w')
                fileCount += 1
            count += 1
            outfp.write('>%s\n%s\n' % (seq.description, str(seq.seq)))
    outfp.close()
    return fileCount, seqCount
761331
def prepare_sample(data):
    
    data = utils.to_single_data(data)
    logger.debug("Preparing %s" % data["rgnames"]["sample"])
    data["files"] = get_fastq_files(data)
    # get_fastq_files swaps over quality scores to standard, unless trimming
    if not(dd.get_trim_reads(data)):
        data = dd.set_quality_format(data, "standard")
    return [[data]]
944603
def cat_files(infiles, outfile):
    ''''''
    f_out = pyfastaq.utils.open_file_write(outfile)

    for filename in infiles:
        if os.path.exists(filename):
            f_in = pyfastaq.utils.open_file_read(filename)
            for line in f_in:
                print(line, end='', file=f_out)
            pyfastaq.utils.close(f_in)

    pyfastaq.utils.close(f_out)
342686
def matchingfiles(self, projectpath):
        """"""
        results = []

        if projectpath[-1] == '/':
            inputpath = projectpath + 'input/'
        else:
            inputpath = projectpath + '/input/'

        for linkf,realf in clam.common.util.globsymlinks(inputpath + '/.*.INPUTTEMPLATE.' + self.id + '.*'):
            seqnr = int(linkf.split('.')[-1])
            results.append( (seqnr, realf[len(inputpath):], self) )
        results = sorted(results)
        if self.unique and len(results) != 1:
            return []
        else:
            return results
39877
def sequence_prep(self):
        
        # Create a sorted list of all the FASTA files in the sequence path
        strains = sorted(glob(os.path.join(self.fastapath, '*.fa*'.format(self.fastapath))))
        for sample in strains:
            # Create the object
            metadata = MetadataObject()
            # Set the sample name to be the file name of the sequence by removing the path and file extension
            sample_name = os.path.splitext(os.path.basename(sample))[0]
            if sample_name in self.strainset:
                # Extract the OLNID from the dictionary using the SEQID
                samplename = self.straindict[sample_name]
                # samplename = sample_name
                # Set and create the output directory
                outputdir = os.path.join(self.path, samplename)
                make_path(outputdir)
                # Set the name of the JSON file
                json_metadata = os.path.join(outputdir, '{name}.json'.format(name=samplename))
                if not os.path.isfile(json_metadata):
                    # Create the name and output directory attributes
                    metadata.name = samplename
                    metadata.seqid = sample_name
                    metadata.outputdir = outputdir
                    metadata.jsonfile = json_metadata
                    # Set the name of the FASTA file to use in the analyses
                    metadata.bestassemblyfile = os.path.join(metadata.outputdir,
                                                             '{name}.fasta'.format(name=metadata.name))
                    # Symlink the original file to the output directory
                    relative_symlink(sample, outputdir, '{sn}.fasta'.format(sn=metadata.name))
                    # Associate the corresponding FASTQ files with the assembly
                    metadata.fastqfiles = sorted(glob(os.path.join(self.fastqpath,
                                                                   '{name}*.gz'.format(name=metadata.name))))
                    metadata.forward_fastq, metadata.reverse_fastq = metadata.fastqfiles
                    # Write the object to file
                    self.write_json(metadata)
                else:
                    metadata = self.read_json(json_metadata)
                # Add the metadata object to the list of objects
                self.metadata.append(metadata)
814889
def get_fastq_files(data):
    
    assert "files" in data, "Did not find `files` in input; nothing to process"
    ready_files = []
    should_gzip = True

    # Bowtie does not accept gzipped fastq
    if 'bowtie' in data['reference'].keys():
        should_gzip = False
    for fname in data["files"]:
        if fname.endswith(".bam"):
            if _pipeline_needs_fastq(data["config"], data):
                ready_files = convert_bam_to_fastq(fname, data["dirs"]["work"],
                                                   data, data["dirs"], data["config"])
            else:
                ready_files = [fname]
        elif objectstore.is_remote(fname):
            ready_files.append(fname)
        # Trimming does quality conversion, so if not doing that, do an explicit conversion
        elif not(dd.get_trim_reads(data)) and dd.get_quality_format(data) != "standard":
            out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "fastq_convert"))
            ready_files.append(fastq.groom(fname, data, out_dir=out_dir))
        else:
            ready_files.append(fname)
    ready_files = [x for x in ready_files if x is not None]
    if should_gzip:
        out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "fastq"))
        ready_files = [_gzip_fastq(x, out_dir) for x in ready_files]
    for in_file in ready_files:
        if not objectstore.is_remote(in_file):
            assert os.path.exists(in_file), "%s does not exist." % in_file
    return ready_files
1289615
def fileprep(self):
        """"""
        # Create and start threads
        for i in range(self.cpus):
            # Send the threads to the appropriate destination function
            threads = Thread(target=self.prep, args=())
            # Set the daemon to true - something to do with thread management
            threads.setDaemon(True)
            # Start the threading
            threads.start()
        for sample in self.metadata:
            # Set the name of the decompressed, combined .fastq file
            sample.general.combined = os.path.join(sample.general.outputdirectory, '{sample_name}_combined.fastq'
                                                   .format(sample_name=sample.name))
            self.queue.put(sample)
        self.queue.join()
532517
def filer(args):
        
        # List to store all the metadata objects
        samples = list()
        # Find all the sequence files in the path
        fastas = sorted(glob(os.path.join(args.sequencepath, '*.fa*')))
        for fasta in fastas:
            # Create a metadata object for each sample
            metadata = MetadataObject()
            # Populate the metadata object with the required attributes
            metadata.name = os.path.splitext(os.path.basename(fasta))[0]
            metadata.general = GenObject()
            metadata.commands = GenObject()
            metadata.general.bestassemblyfile = fasta
            metadata.general.outputdirectory = os.path.join(args.sequencepath, metadata.name)
            samples.append(metadata)
        return samples
548696
def fastq_to_csv(in_file, fastq_format, work_dir):
    
    out_file = "%s.csv" % (os.path.splitext(os.path.basename(in_file))[0])
    out_file = os.path.join(work_dir, out_file)
    if not (os.path.exists(out_file) and os.path.getsize(out_file) > 0):
        with open(in_file) as in_handle:
            with open(out_file, "w") as out_handle:
                writer = csv.writer(out_handle)
                for rec in SeqIO.parse(in_handle, fastq_format):
                    writer.writerow([rec.id] + rec.letter_annotations["phred_quality"])
    return out_file
973393
def _run_cortex(fastq, indexes, params, out_base, dirs, config):
    
    print(out_base)
    fastaq_index = "{0}.fastaq_index".format(out_base)
    se_fastq_index = "{0}.se_fastq".format(out_base)
    pe_fastq_index = "{0}.pe_fastq".format(out_base)
    reffasta_index = "{0}.list_ref_fasta".format(out_base)
    with open(se_fastq_index, "w") as out_handle:
        out_handle.write(fastq + "\n")
    with open(pe_fastq_index, "w") as out_handle:
        out_handle.write("")
    with open(fastaq_index, "w") as out_handle:
        out_handle.write("{0}\t{1}\t{2}\t{2}\n".format(params["sample"], se_fastq_index,
                                                       pe_fastq_index))
    with open(reffasta_index, "w") as out_handle:
        for x in indexes["fasta"]:
            out_handle.write(x + "\n")
    os.environ["PERL5LIB"] = "{0}:{1}:{2}".format(
        os.path.join(dirs["cortex"], "scripts/calling"),
        os.path.join(dirs["cortex"], "scripts/analyse_variants/bioinf-perl/lib"),
        os.environ.get("PERL5LIB", ""))
    kmers = sorted(params["kmers"])
    kmer_info = ["--first_kmer", str(kmers[0])]
    if len(kmers) > 1:
        kmer_info += ["--last_kmer", str(kmers[-1]),
                      "--kmer_step", str(kmers[1] - kmers[0])]
    subprocess.check_call(["perl", os.path.join(dirs["cortex"], "scripts", "calling", "run_calls.pl"),
                           "--fastaq_index", fastaq_index,
                           "--auto_cleaning", "yes", "--bc", "yes", "--pd", "yes",
                           "--outdir", os.path.dirname(out_base), "--outvcf", os.path.basename(out_base),
                           "--ploidy", str(config["algorithm"].get("ploidy", 2)),
                           "--stampy_hash", indexes["stampy"],
                           "--stampy_bin", os.path.join(dirs["stampy"], "stampy.py"),
                           "--refbindir", os.path.dirname(indexes["cortex"][0]),
                           "--list_ref_fasta",  reffasta_index,
                           "--genome_size", str(params["genome_size"]),
                           "--max_read_len", "30000",
                           #"--max_var_len", "4000",
                           "--format", "FASTQ", "--qthresh", "5", "--do_union", "yes",
                           "--mem_height", "17", "--mem_width", "100",
                           "--ref", "CoordinatesAndInCalling", "--workflow", "independent",
                           "--vcftools_dir", dirs["vcftools"],
                           "--logfile", "{0}.logfile,f".format(out_base)]
                          + kmer_info)
    final = glob.glob(os.path.join(os.path.dirname(out_base), "vcfs",
                                   "{0}*FINALcombined_BC*decomp.vcf".format(os.path.basename(out_base))))
    # No calls, need to setup an empty file
    if len(final) != 1:
        print("Did not find output VCF file for {0}".format(out_base))
        return None
    else:
        return final[0]
1354632
def run(bam_file, sample, out_dir):
    """"""
    out = {}
    # if "rchipqc" in dd.get_tools_on(sample):
    #    out = chipqc(bam_file, sample, out_dir)

    peaks = sample.get("peaks_files", {}).get("main")
    if peaks:
        out.update(_reads_in_peaks(bam_file, peaks, sample))
    return out
720725
def demultiplex_samples(data):
    
    work_dir = os.path.join(dd.get_work_dir(data), "umis")
    sample_dir = os.path.join(work_dir, dd.get_sample_name(data))
    demulti_dir = os.path.join(sample_dir, "demultiplexed")

    files = data["files"]
    if len(files) == 2:
        logger.error("Sample demultiplexing doesn't handle paired-end reads, but "
            "we can add it. Open an issue here https://github.com/bcbio/bcbio-nextgen/issues if you need this and we'll add it.")
        sys.exit(1)
    else:
        fq1 = files[0]
    # check if samples need to be demultiplexed
    with open_fastq(fq1) as in_handle:
        read = next(in_handle)
        if "SAMPLE_" not in read:
            return [[data]]

    bcfile = get_sample_barcodes(dd.get_sample_barcodes(data), sample_dir)
    demultiplexed = glob.glob(os.path.join(demulti_dir, "*.fq*"))
    if demultiplexed:
        return [split_demultiplexed_sampledata(data, demultiplexed)]
    umis = _umis_cmd(data)
    cmd = ("{umis} demultiplex_samples --nedit 1 --barcodes {bcfile} "
           "--out_dir {tx_dir} {fq1}")
    msg = "Demultiplexing {fq1}."
    with file_transaction(data, demulti_dir) as tx_dir:
        do.run(cmd.format(**locals()), msg.format(**locals()))
    demultiplexed = glob.glob(os.path.join(demulti_dir, "*.fq*"))
    return [split_demultiplexed_sampledata(data, demultiplexed)]
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))
1346064
def load_fasta_file_as_dict_of_seqrecords(filename):
    

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

    return results
164656
def fastqfilter(self):
        """"""
        printtime('Creating filtered .fastqfiles', self.start)
        # Create and start threads
        for i in range(self.cpus):
            # Send the threads to the appropriate destination function
            threads = Thread(target=self.filterfastq, args=())
            # Set the daemon to true - something to do with thread management
            threads.setDaemon(True)
            # Start the threading
            threads.start()
        for sample in self.runmetadata.samples:
            self.filterqueue.put(sample)
        self.filterqueue.join()
        # Print the metadata to file
        metadataprinter.MetadataPrinter(self)