This is a Toil pipeline for the UNC best practice RNA-Seq analysis.
    RNA-seq fastqs are combined, aligned, sorted, filtered, and quantified.

    Please read the README.md located in the same directory.
321262
def main():
    
    # Define Parser object and add to toil
    parser = build_parser()
    Job.Runner.addToilOptions(parser)
    args = parser.parse_args()
    # Store inputs from argparse
    inputs = {'config': args.config,
              'config_fastq': args.config_fastq,
              'input': args.input,
              'unc.bed': args.unc,
              'hg19.transcripts.fa': args.fasta,
              'composite_exons.bed': args.composite_exons,
              'normalize.pl': args.normalize,
              'output_dir': args.output_dir,
              'rsem_ref.zip': args.rsem_ref,
              'chromosomes.zip': args.chromosomes,
              'ebwt.zip': args.ebwt,
              'ssec': args.ssec,
              's3_dir': args.s3_dir,
              'sudo': args.sudo,
              'single_end_reads': args.single_end_reads,
              'upload_bam_to_s3': args.upload_bam_to_s3,
              'uuid': None,
              'sample.tar': None,
              'cpu_count': None}

    # Launch jobs
    Job.Runner.startToil(Job.wrapJobFn(download_shared_files, inputs), args)
1074404
def align_transcriptome(fastq_file, pair_file, ref_file, data):
    
    work_bam = dd.get_work_bam(data)
    base, ext = os.path.splitext(work_bam)
    out_file = base + ".transcriptome" + ext
    if utils.file_exists(out_file):
        data = dd.set_transcriptome_bam(data, out_file)
        return data
    # bwa mem needs phred+33 quality, so convert if it is Illumina
    if dd.get_quality_format(data).lower() == "illumina":
        logger.info("bwa mem does not support the phred+64 quality format, "
                    "converting %s and %s to phred+33.")
        fastq_file = fastq.groom(fastq_file, data, in_qual="fastq-illumina")
        if pair_file:
            pair_file = fastq.groom(pair_file, data, in_qual="fastq-illumina")
    bwa = config_utils.get_program("bwa", data["config"])
    gtf_file = dd.get_gtf_file(data)
    gtf_fasta = index_transcriptome(gtf_file, ref_file, data)
    args = " ".join(_bwa_args_from_config(data["config"]))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    samtools = config_utils.get_program("samtools", data["config"])
    cmd = ("{bwa} mem {args} -a -t {num_cores} {gtf_fasta} {fastq_file} "
           "{pair_file} ")
    with file_transaction(data, out_file) as tx_out_file:
        message = "Aligning %s and %s to the transcriptome." % (fastq_file, pair_file)
        cmd += "| " + postalign.sam_to_sortbam_cl(data, tx_out_file, name_sort=True)
        do.run(cmd.format(**locals()), message)
    data = dd.set_transcriptome_bam(data, out_file)
    return data
12299
def map_alignment(self, prot_seq, nucl_seq):
        
        if prot_seq.id != nucl_seq.id:
            logging.warning(
                'ID mismatch: %s != %s. Are the sequences in the same order?',
                prot_seq.id, nucl_seq.id)

        # Ungap nucleotides
        codons = batch(str(nucl_seq.seq.ungap('-')), 3)
        codons = [''.join(i) for i in codons]
        codon_iter = iter(codons)

        ungapped_prot = str(prot_seq.seq).replace('-', '')

        if len(ungapped_prot) != len(codons):
            table = self.translation_table.forward_table
            prot_str = ' '.join(' ' + p + ' ' for p in ungapped_prot)
            codon_str = ' '.join(codons)
            trans_str = ' '.join(' ' + table.get(codon, 'X') + ' '
                                 for codon in codons)
            raise ValueError("""Length of codon sequence ({0}) does not match \
length of protein sequence ({1}) for {2}
Protein:       {3}
Codons:        {4}
Trans. Codons: {5}""".format(len(codons), len(ungapped_prot), nucl_seq.id, prot_str,
                            codon_str, trans_str))

        try:
            nucl_align = ['---' if p == '-' else next(codon_iter)
                          for p in str(prot_seq.seq)]
        except StopIteration:
            assert False  # Should be checked above

        result = SeqRecord(Seq(''.join(nucl_align)), id=nucl_seq.id,
                           description=nucl_seq.description)

        # Validate
        self._validate_translation(prot_seq.seq.upper(), result.seq.upper())

        return result
332779
def read_sequencing(directory):
    '''

    '''
    dirfiles = os.listdir(directory)
    seq_exts = ['.seq', '.abi', '.ab1']
    # Exclude files that aren't sequencing results
    seq_paths = [x for x in dirfiles if os.path.splitext(x)[1] in seq_exts]
    paths = [os.path.join(directory, x) for x in seq_paths]
    sequences = [read_dna(x) for x in paths]

    return sequences
1262295
def transcript_sequence(self, transcript_id):
        
        if self.transcript_sequences is None:
            raise ValueError(
                "No transcript FASTA supplied to this Genome: %s" % self)
        return self.transcript_sequences.get(transcript_id)
79606
def find_orfs(fa, seqs):
    
    faa = '%s.prodigal.faa' % (fa)
    fna = '%s.prodigal.fna' % (fa)
    gbk = '%s.prodigal.gbk' % (fa)
    if os.path.exists(faa) is False:
        p = subprocess.Popen('prodigal -q -i %s -a %s -d %s -c -f gbk -m -n -o %s -p meta' \
                % (fa, faa, fna, gbk), shell = True)
        p.communicate()
    for orf in parse_fasta(faa):
        if orf[0] == []:
            continue
        id = orf[0].split('>')[1].split('_', 1)[0]
        pos = sorted([int(i) for i in orf[0].split()[2:5] if i != '#'])
        if id not in seqs:
            continue
        for i, ins in enumerate(seqs[id][2]):
            if check_overlap(pos, ins, 0.90) is True:
                seqs[id][2][i][4].append(orf)
    return seqs, faa
576175
def gatk_rnaseq_calling(data):
    
    from bcbio.bam import callable
    data = utils.deepish_copy(data)
    tools_on = dd.get_tools_on(data)
    if not tools_on:
        tools_on = []
    tools_on.append("gvcf")
    data = dd.set_tools_on(data, tools_on)
    data = dd.set_jointcaller(data, ["%s-joint" % v for v in dd.get_variantcaller(data)])
    out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data),
                                              "variation", "rnaseq", "gatk-haplotype"))
    data = _setup_variant_regions(data, out_dir)
    out_file = os.path.join(out_dir, "%s-gatk-haplotype.vcf.gz" % dd.get_sample_name(data))
    if not utils.file_exists(out_file):
        region_files = []
        regions = []
        for cur_region in callable.get_split_regions(dd.get_variant_regions(data), data):
            str_region = "_".join([str(x) for x in cur_region])
            region_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data),
                                                                    "variation", "rnaseq", "gatk-haplotype",
                                                                    "regions")),
                                    "%s-%s-gatk-haplotype.vcf.gz" % (dd.get_sample_name(data), str_region))
            region_file = gatk.haplotype_caller([dd.get_split_bam(data)], [data], dd.get_ref_file(data), {},
                                                region=cur_region, out_file=region_file)
            region_files.append(region_file)
            regions.append(cur_region)
        out_file = vcfutils.concat_variant_files(region_files, out_file, regions,
                                                 dd.get_ref_file(data), data["config"])
    return dd.set_vrn_file(data, out_file)
1100333
def process_uclust_pw_alignment_results(fasta_pairs_lines, uc_lines):
    
    alignments = get_next_two_fasta_records(fasta_pairs_lines)
    for hit in get_next_record_type(uc_lines, 'H'):
        matching_strand = hit[4]
        if matching_strand == '-':
            strand_id = '-'
            target_rev_match = True
        elif matching_strand == '+':
            strand_id = '+'
            target_rev_match = False
        elif matching_strand == '.':
            # protein sequence, so no strand information
            strand_id = ''
            target_rev_match = False
        else:
            raise UclustParseError("Unknown strand type: %s" % matching_strand)
        uc_query_id = hit[8]
        uc_target_id = hit[9]
        percent_id = float(hit[3])

        fasta_pair = alignments.next()

        fasta_query_id = fasta_pair[0][0]
        aligned_query = fasta_pair[0][1]

        if fasta_query_id != uc_query_id:
            raise UclustParseError("Order of fasta and uc files do not match." +
                                   " Got query %s but expected %s." %
                                   (fasta_query_id, uc_query_id))

        fasta_target_id = fasta_pair[1][0]
        aligned_target = fasta_pair[1][1]

        if fasta_target_id != uc_target_id + strand_id:
            raise UclustParseError("Order of fasta and uc files do not match." +
                                   " Got target %s but expected %s." %
                                   (fasta_target_id, uc_target_id + strand_id))

        if target_rev_match:
            query_id = uc_query_id + ' RC'
            aligned_query = DNA.rc(aligned_query)
            target_id = uc_target_id
            aligned_target = DNA.rc(aligned_target)
        else:
            query_id = uc_query_id
            aligned_query = aligned_query
            target_id = uc_target_id
            aligned_target = aligned_target

        yield (query_id, target_id, aligned_query, aligned_target, percent_id)
769588
def cb_filter(fastq, bc1, bc2, bc3, cores, nedit):
    ''' 
    '''

    with open_gzipsafe(bc1) as bc1_fh:
        bc1 = set(cb.strip() for cb in bc1_fh)
    if bc2:
        with open_gzipsafe(bc2) as bc2_fh:
            bc2 = set(cb.strip() for cb in bc2_fh)
    if bc3:
        with open_gzipsafe(bc3) as bc3_fh:
            bc3 = set(cb.strip() for cb in bc3_fh)

    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)

    if nedit == 0:
        filter_cb = partial(exact_barcode_filter, bc1=bc1, bc2=bc2, bc3=bc3,
                            re_string=re_string)
    else:
        bc1hash = MutationHash(bc1, nedit)
        bc2hash = None
        bc3hash = None
        if bc2:
            bc2hash = MutationHash(bc2, nedit)
        if bc3:
            bc3hash = MutationHash(bc3, nedit)
        filter_cb = partial(correcting_barcode_filter, bc1hash=bc1hash,
                            bc2hash=bc2hash, bc3hash=bc3hash, re_string=re_string)
    p = multiprocessing.Pool(cores)

    chunks = tz.partition_all(10000, read_fastq(fastq))
    bigchunks = tz.partition_all(cores, chunks)
    for bigchunk in bigchunks:
        for chunk in p.map(filter_cb, list(bigchunk)):
            for read in chunk:
                sys.stdout.write(read)
1187703
def main():
    
    parser = argparse.ArgumentParser(
        description='monoseq: pretty-printing DNA and protein sequences',
        epilog='If INPUT is in FASTA format, each record is pretty-printed '
        'after printing its name and ANNOTATION (if supplied) is used by '
        'matching chromosome/record name. If INPUT contains a raw sequence, '
        'only the first chromosome in ANNOTATION is used.')
    parser.add_argument(
        'sequence_file', metavar='INPUT', nargs='?', default=sys.stdin,
        type=argparse.FileType('r'), help='file to read sequence(s) from, '
        'can be in FASTA format (default: standard input)')
    parser.add_argument(
        '-b', '--block-length', metavar='LENGTH', dest='block_length',
        type=int, default=10, help='block length in letters (default: 10)')
    parser.add_argument(
        '-l', '--blocks-per-line', metavar='BLOCKS', dest='blocks_per_line',
        type=int, default=6, help='blocks per line (default: 6)')
    parser.add_argument(
        '-a', '--annotation', metavar='POS', dest='annotation', nargs=2,
        action='append', type=int, help='first and last positions of '
        'subsequence to annotate (allowed more than once)')
    parser.add_argument(
        '-e', '--bed', metavar='ANNOTATION', dest='annotation_file',
        type=argparse.FileType('r'), help='file to read annotation from in '
        'BED format')

    args = parser.parse_args()
    pprint(_until_eof(args.sequence_file), annotation=args.annotation,
           annotation_file=args.annotation_file,
           block_length=args.block_length,
           blocks_per_line=args.blocks_per_line)
293558
def parse_rna(rna, s2bins, min_rna):
    
    rna_cov = {}
    for seq in parse_fasta(rna):
        # check that length passes threshold
        length = len(seq[1])
        if length < min_rna:
            continue
        # check if sequence is binnned
        s = seq[0].split('>')[1].split()[0]
        if s not in s2bins:
            continue
        if s not in rna_cov:
            rna_cov[s] = [0, 0, []]
        position = [int(i) for i in seq[0].rsplit('pos=', 1)[1].split()[0].split('-')]
        rna_cov[s][2].append(position)
        rna_cov[s][1] += length
    return rna_cov
94703
def find_introns(fa, seqs, sequences, threads):
    
    db = '%s/rfam/Rfam.cm.1_1.introns' % (os.environ['databases'])
    out = '%s.rfam-introns.cmscan' % (fa)
    tblout = '%s.rfam-introns.cmscan-tblout' % (fa)
    if os.path.exists(tblout) is False:
        p = subprocess.Popen('cmscan --cpu %s --tblout %s %s %s > %s'\
                % (threads, tblout, db, fa, out), shell = True)
        p.communicate()
    for line in open(tblout):
        if line.startswith('#'):
            continue
        line = line.strip().split()
        if line[16] == '?': # does not pass inclusion threshold
            continue
        id = line[2]
        type, start, stop, strand = line[0], int(line[7]), int(line[8]), line[9]
        if 'intron' not in type.lower():
            continue
        pos = sorted([start, stop])
        if id not in seqs:
            continue
        for i, ins in enumerate(seqs[id][2]):
            if check_overlap(pos, ins, 0.25) is True:
                seqs[id][2][i][5].append(['>%s_%s %s %s %s-%s' % (id, (i + 1), type, strand, start, stop), sequences[id][1][pos[0]-1:pos[1]]])
    return seqs
163510
def extract_ngram_filter(pos_seq, regex=SimpleNP, minlen=1, maxlen=8):
	""""""
	ss = coarse_tag_str(pos_seq)

	def gen():
		for s in xrange(len(ss)):
			for n in xrange(minlen, 1 + min(maxlen, len(ss) - s)):
				e = s + n
				substr = ss[s:e]
				if re.match(regex + "$", substr):
					yield (s, e)

	return list(gen())
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)
369046
def combinealleles(self, alleles):
        
        printtime('Creating combined core allele file', self.start)
        if not os.path.isfile(os.path.join(self.coregenelocation, 'core_combined.tfa')):
            with open(os.path.join(self.coregenelocation, 'core_combined.tfa'), 'w') as combinedfile:
                # Open each allele file
                for allele in sorted(alleles):
                    for record in SeqIO.parse(open(allele, "rU"), "fasta"):
                        # Extract the sequence record from each entry in the multifasta
                        # Remove and dashes or 'N's from the sequence data - makeblastdb can't handle sequences
                        # with gaps
                        # noinspection PyProtectedMember
                        record.seq._data = record.seq._data.replace('-', '').replace('N', '')
                        # Clear the name and description attributes of the record
                        record.name = ''
                        record.description = ''
                        # Write each record to the combined file
                        SeqIO.write(record, combinedfile, 'fasta')
237241
def blat(self, db=None, sequence=None, seq_type="DNA"):
        
        from . blat_blast import blat, blat_all
        assert sequence in (None, "cds", "mrna")
        seq = self.sequence() if sequence is None else ("".join(self.cds_sequence if sequence == "cds" else self.mrna_sequence))
        if isinstance(db, (tuple, list)):
            return blat_all(seq, self.gene_name, db, seq_type)
        else:
            return blat(seq, self.gene_name, db or self.db, seq_type)
1062372
def parse_fasta_annotations(fastas, annot_tables, trans_table):
    
    if annot_tables is not False:
        annots = {}
        for table in annot_tables:
            for cds in open(table):
                ID, start, end, strand = cds.strip().split()
                annots[ID] = [start, end, int(strand)]
    for fasta in fastas:
        for seq in parse_fasta(fasta):
            if ('# ;gc_cont' not in seq[0] and '# ID=' not in seq[0]) and annot_tables is False:
                print('# specify fasta from Prodigal or annotations table (-t)', file=sys.stderr)
                exit()
            if 'ID=' in seq[0]:
                ID = seq[0].rsplit('ID=', 1)[1].split(';', 1)[0]
                contig = seq[0].split()[0].split('>')[1].rsplit('_%s' % (ID), 1)[0]
            else:
                contig = seq[0].split()[0].split('>')[1].rsplit('_', 1)[0]
            locus = seq[0].split()[0].split('>')[1]
            # annotation info from Prodigal
            if ('# ;gc_cont' in seq[0] or '# ID=' in seq[0]):
                info = seq[0].split(' # ')
                start, end, strand = int(info[1]), int(info[2]), info[3]
                if strand == '1':
                    strand = 1
                else:
                    strand = -1
                product = [''.join(info[4].split()[1:])]
            # annotation info from table
            else:
                start, end, strand = annots[locus]
                product = seq[0].split(' ', 1)[1]
            info = {'transl_table':[trans_table], \
                    'translation':[seq[1]], \
                    'product':product}
            yield contig, [locus, [start, end, strand], info]