Test if reads can get aligned to the lambda genome,
    if not: write to stdout
483946
def align(aligner, reads):
    '''
    
    '''
    i = 0
    for record in SeqIO.parse(reads, "fastq"):
        try:
            next(aligner.map(str(record.seq)))
            i += 1
        except StopIteration:
            print(record.format("fastq"), end='')
    sys.stderr.write("NanoLyse: removed {} reads.\n".format(i))
298556
def CDSReader( fh, format='gff' ):
    

    known_formats = ( 'gff', 'gtf', 'bed')
    if format not in known_formats: 
        print('%s format not in %s' % (format, ",".join( known_formats )), file=sys.stderr)
        raise Exception('?')

    if format == 'bed':
        for line in fh:    
            f = line.strip().split()
            chrom = f[0]
            chrom_start = int(f[1])
            name = f[4]
            strand = f[5]
            cdsStart = int(f[6])
            cdsEnd = int(f[7])
            blockCount = int(f[9])
            blockSizes = [ int(i) for i in f[10].strip(',').split(',') ]
            blockStarts = [ chrom_start + int(i) for i in f[11].strip(',').split(',') ]

            # grab cdsStart - cdsEnd
            cds_exons = []
            cds_seq = ''
            genome_seq_index = []
            for base,offset in zip( blockStarts, blockSizes ):
                if (base + offset) < cdsStart: continue
                if base > cdsEnd: continue
                exon_start = max( base, cdsStart )
                exon_end = min( base+offset, cdsEnd ) 
                cds_exons.append( (exon_start, exon_end) )
            

    genelist = {}
    grouplist = []
    if format == 'gff' or format == 'gtf':
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len( fields ) < 9: continue
            if fields[2] not in ('CDS', 'stop_codon', 'start_codon'): continue

            # fields

            chrom = fields[0]
            ex_st = int( fields[3] ) - 1 # make zero-centered
            ex_end = int( fields[4] ) #+ 1 # make exclusive
            strand = fields[6]

            if format == 'gtf':
                group = fields[8].split(';')[0]
            else:
                group = fields[8]

            if group not in grouplist: grouplist.append( group )
            if group not in genelist:
                genelist[group] = (chrom, strand, [])
            
            genelist[group][2].append( ( ex_st, ex_end ) )

        sp = lambda a,b: cmp( a[0], b[0] )

        #for gene in genelist.values():
        for gene in grouplist:
            chrom, strand, cds_exons = genelist[ gene ]
            seqlen = sum([ a[1]-a[0] for a in cds_exons ])
            overhang = seqlen % 3
            if overhang > 0:
                #print >>sys.stderr, "adjusting ", gene  
                if strand == '+': 
                    cds_exons[-1] = ( cds_exons[-1][0], cds_exons[-1][1] - overhang )
                else:
                    cds_exons[0] = ( cds_exons[0][0] + overhang, cds_exons[0][1] )
            cds_exons = bitset_union( cds_exons )
            yield chrom, strand, cds_exons, gene
20549
def _run(self, stdout=subprocess.PIPE):
        '''
        
        '''
        assert self._pop is None
        if self._input is not None:
            self._input._run()

        cwd = (self._cwd if self._cwd else os.getcwd()) + '/'

        def glob_or(expr):
            res = [expr]
            if self._expand and ('*' in expr or '?' in expr):
                if expr and expr[0] != '/':
                    expr = cwd + expr
                res = glob.glob(expr)
                if not res:
                    res = [expr]
            return res

        args = [glob_or(arg) for arg in self._args]
        stderr = subprocess.STDOUT if self._err_to_out else None
        assert not self._stdin or self._input is None, "*** ERROR: Stdin must be none outside of head of chain."
        stdin = self._stdin if self._input is None else self._input._pop.stdout
        if stdin is not None and getattr(stdin, 'fileno', None) is None:
            stdin = subprocess.PIPE
        cmd = [self._cmd] + list(itertools.chain.from_iterable(args))
        if self.debug:
            print "*** DEBUG: Run", cmd
        try:
            self._pop = subprocess.Popen(cmd, stdin=stdin, stdout=stdout, stderr=stderr,
                                         close_fds=True, cwd=self._cwd, env=self._env)
            if stdin == subprocess.PIPE:
                try:
                    # TODO: Write a non-blocking reader-writer to ship data here
                    self._pop.stdin.write(self._stdin.read())
                except AttributeError:
                    self._pop.stdin.write(self._stdin)
        except OSError:
            print "*** ERROR: %r failed at %r" % (self, self._original_cmd)
            raise
177282
def parallel(args):
    
    from jcvi.formats.base import split

    p = OptionParser(parallel.__doc__)
    p.set_home("maker")
    p.set_tmpdir(tmpdir="tmp")
    p.set_grid_opts(array=True)
    opts, args = p.parse_args(args)

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

    genome, NN = args
    threaded = opts.threaded or 1
    tmpdir = opts.tmpdir

    mkdir(tmpdir)
    tmpdir = get_abs_path(tmpdir)

    N = int(NN)
    assert 1 <= N < 1000, "Required: 1 < N < 1000!"

    outdir = "outdir"
    fs = split([genome, outdir, NN])

    c = CTLFile("maker_opts.ctl")
    c.update_abs_path()
    if threaded > 1:
        c.update_tag("cpus", threaded)

    cwd = os.getcwd()
    dirs = []
    for name in fs.names:
        fn = get_abs_path(name)
        bn = op.basename(name)
        dirs.append(bn)
        c.update_tag("genome", fn)
        mkdir(bn)
        sh("cp *.ctl {0}".format(bn))

        os.chdir(bn)
        c.write_file("maker_opts.ctl")
        os.chdir(cwd)

    jobs = "jobs"
    fw = open(jobs, "w")
    print("\n".join(dirs), file=fw)
    fw.close()

    # Submit to grid
    ncmds = len(dirs)
    runfile = "array.sh"
    cmd = op.join(opts.maker_home, "bin/maker")
    if tmpdir:
        cmd += " -TMP {0}".format(tmpdir)

    engine = get_grid_engine()
    contents = arraysh.format(jobs, cmd) if engine == "SGE" \
                else arraysh_ua.format(N, threaded, jobs, cmd)
    write_file(runfile, contents)

    if engine == "PBS":
        return

    # qsub script
    outfile = "maker.\$TASK_ID.out"
    p = GridProcess(runfile, outfile=outfile, errfile=outfile,
                    arr=ncmds, grid_opts=opts)
    qsubfile = "qsub.sh"
    qsub = p.build()
    write_file(qsubfile, qsub)
1151992
def _prep_genome(out_dir, data):
    
    genome_name = utils.splitext_plus(os.path.basename(dd.get_ref_file(data)))[0]
    out_dir = utils.safe_makedir(os.path.join(out_dir, genome_name))
    ref_file = dd.get_ref_file(data)
    utils.symlink_plus(ref_file, os.path.join(out_dir, os.path.basename(ref_file)))
    with open(os.path.join(out_dir, "GenomeSize.xml"), "w") as out_handle:
        out_handle.write('<sequenceSizes genomeName="%s">' % genome_name)
        for c in pysam.AlignmentFile("%s.dict" % utils.splitext_plus(ref_file)[0]).header["SQ"]:
            cur_ploidy = ploidy.get_ploidy([data], region=[c["SN"]])
            out_handle.write('<chromosome fileName="%s" contigName="%s" totalBases="%s" knownBases="%s" '
                             'isCircular="false" ploidy="%s" md5="%s"/>' %
                             (os.path.basename(ref_file), c["SN"], c["LN"], c["LN"], cur_ploidy, c["M5"]))
        out_handle.write('</sequenceSizes>')
    return out_dir
1114717
def alignment(job, ids, input_args, sample):
    
    uuid, urls = sample
    # ids['bam'] = job.fileStore.getEmptyFileStoreID()
    work_dir = job.fileStore.getLocalTempDir()
    output_dir = input_args['output_dir']
    key_path = input_args['ssec']
    cores = multiprocessing.cpu_count()

    # I/O
    return_input_paths(job, work_dir, ids, 'ref.fa', 'ref.fa.amb', 'ref.fa.ann',
                                                     'ref.fa.bwt', 'ref.fa.pac', 'ref.fa.sa', 'ref.fa.fai')
    # Get fastqs associated with this sample
    for url in urls:
        download_encrypted_file(work_dir, url, key_path, os.path.basename(url))

    # Parameters for BWA and Bamsort
    docker_cmd = ['docker', 'run', '--rm', '-v', '{}:/data'.format(work_dir)]

    bwa_command = ["jvivian/bwa",
                   "mem",
                   "-R", "@RG\tID:{0}\tPL:Illumina\tSM:{0}\tLB:KapaHyper".format(uuid),
                   "-T", str(0),
                   "-t", str(cores),
                   "/data/ref.fa"] + [os.path.join('/data/',  os.path.basename(x)) for x in urls]

    bamsort_command = ["jeltje/biobambam",
                       "/usr/local/bin/bamsort",
                       "inputformat=sam",
                       "level=1",
                       "inputthreads={}".format(cores),
                       "outputthreads={}".format(cores),
                       "calmdnm=1",
                       "calmdnmrecompindetonly=1",
                       "calmdnmreference=/data/ref.fa",
                       "I=/data/{}".format(uuid + '.sam')]
    # Piping the output to a file handle
    with open(os.path.join(work_dir, uuid + '.sam'), 'w') as f_out:
        subprocess.check_call(docker_cmd + bwa_command, stdout=f_out)

    with open(os.path.join(work_dir, uuid + '.bam'), 'w') as f_out:
        subprocess.check_call(docker_cmd + bamsort_command, stdout=f_out)

    # Save in JobStore
    # job.fileStore.updateGlobalFile(ids['bam'], os.path.join(work_dir, uuid + '.bam'))
    ids['bam'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, uuid + '.bam'))
    # Copy file to S3
    if input_args['s3_dir']:
        job.addChildJobFn(upload_bam_to_s3, ids, input_args, sample, cores=32, memory='20 G', disk='30 G')
    # Move file in output_dir
    if input_args['output_dir']:
        move_to_output_dir(work_dir, output_dir, uuid=None, files=[uuid + '.bam'])
451196
def print_seqs(seqs, id2name, name):
    
    orfs = open('%s.orfs.faa' % (name.rsplit('.', 1)[0]), 'w')
    introns = open('%s.introns.fa' % (name.rsplit('.', 1)[0]), 'w')
    insertions = open('%s.insertions.fa' % (name.rsplit('.', 1)[0]), 'w')
    for seq in seqs:
        for i, ins in enumerate(seqs[seq][2], 1):
            model_pos = ins[1]
            if ins[6] is True: # orf(s) in ins[4]
                for orf in ins[4]:
                    orf_info = orf[0].split('>')[1].split()
                    id = orf_info[0].split('_', 1)[0]
                    name = id2name[id]
                    annotation = orf_info[1]
                    strand = orf_info[7]
                    if strand == '1':
                        strand = '+'
                    else:
                        strand = '-'
                    start, stop = sorted([int(orf_info[3]), int(orf_info[5])])
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s-%s model-pos=%s'\
                            % (name, i, annotation, strand, start, stop, model_pos)
                    print('\n'.join([header, orf[1]]), file=orfs)
            if ins[7] is True: # intron(s) in ins[5]
                for intron in ins[5]:
                    id, type, strand, pos = intron[0].split('>', 1)[1].split()
                    name = id2name[id.split('_')[0]]
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s model-pos=%s'\
                            % (name, i, type, strand, pos, model_pos)
                    print('\n'.join([header, intron[1]]), file=introns)
            insertion = ins[3]
            id, info = insertion[0].split('>')[1].split(' ', 1)
            name = id2name[id.split('_')[0]]
            header = '>%s %s' % (name, info)
            print('\n'.join([header, insertion[1]]), file=insertions)
    insertions.close()        
    orfs.close()
    introns.close()
125231
def alignUnalignedSeqsFromFile(filename,WorkingDir=None,SuppressStderr=None,\
    SuppressStdout=None):
    
    app = Clustalw(WorkingDir=WorkingDir,SuppressStderr=SuppressStderr,\
        SuppressStdout=SuppressStdout)
    return app(filename)
921357
def align(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(align.__doc__)
    p.add_option("--rnaseq", default=False, action="store_true",
                 help="Input is RNA-seq reads, turn splicing on")
    p.add_option("--native", default=False, action="store_true",
                 help="Convert GSNAP output to NATIVE format")
    p.set_home("eddyyeh")
    p.set_outdir()
    p.set_cpus()
    opts, args = p.parse_args(args)

    if len(args) == 2:
        logging.debug("Single-end alignment")
    elif len(args) == 3:
        logging.debug("Paired-end alignment")
    else:
        sys.exit(not p.print_help())

    dbfile, readfile = args[:2]
    outdir = opts.outdir
    assert op.exists(dbfile) and op.exists(readfile)
    prefix = get_prefix(readfile, dbfile)
    logfile = op.join(outdir, prefix + ".log")
    gsnapfile = op.join(outdir, prefix + ".gsnap")
    nativefile = gsnapfile.rsplit(".", 1)[0] + ".unique.native"

    if not need_update((dbfile, readfile), gsnapfile):
        logging.error("`{0}` exists. `gsnap` already run.".format(gsnapfile))
    else:
        dbdir, dbname = check_index(dbfile)
        cmd = "gsnap -D {0} -d {1}".format(dbdir, dbname)
        cmd += " -B 5 -m 0.1 -i 2 -n 3"  # memory, mismatch, indel penalty, nhits
        if opts.rnaseq:
            cmd += " -N 1"
        cmd += " -t {0}".format(opts.cpus)
        cmd += " --gmap-mode none --nofails"
        if readfile.endswith(".gz"):
            cmd += " --gunzip"
        try:
            offset = "sanger" if guessoffset([readfile]) == 33 else "illumina"
            cmd += " --quality-protocol {0}".format(offset)
        except AssertionError:
            pass
        cmd += " " + " ".join(args[1:])
        sh(cmd, outfile=gsnapfile, errfile=logfile)

    if opts.native:
        EYHOME = opts.eddyyeh_home
        if need_update(gsnapfile, nativefile):
            cmd = op.join(EYHOME, "convert2native.pl")
            cmd += " --gsnap {0} -o {1}".format(gsnapfile, nativefile)
            cmd += " -proc {0}".format(opts.cpus)
            sh(cmd)

    return gsnapfile, logfile
1057671
def _get_seqprop_to_structprop_alignment(self, seqprop, structprop, chain_id):
        """"""
        full_structure_id = '{}-{}'.format(structprop.id, chain_id)
        aln_id = '{}_{}'.format(seqprop.id, full_structure_id)

        if self.sequence_alignments.has_id(aln_id):
            alignment = self.sequence_alignments.get_by_id(aln_id)
            return alignment
        else:
            raise ValueError('{}: structure alignment not found, please run the alignment first'.format(aln_id))
864298
def build_ml_phyml(alignment, outfile, work_dir=".", **kwargs):
    
    phy_file = op.join(work_dir, "work", "aln.phy")
    AlignIO.write(alignment, file(phy_file, "w"), "phylip-relaxed")

    phyml_cl = PhymlCommandline(cmd=PHYML_BIN("phyml"), input=phy_file, **kwargs)
    logging.debug("Building ML tree using PhyML: %s" % phyml_cl)
    stdout, stderr = phyml_cl()

    tree_file = phy_file + "_phyml_tree.txt"
    if not op.exists(tree_file):
        print("***PhyML failed.", file=sys.stderr)
        return None
    sh("cp {0} {1}".format(tree_file, outfile), log=False)

    logging.debug("ML tree printed to %s" % outfile)

    return outfile, phy_file
955413
def _build_index(maf_strm, ref_spec):
  """"""
  idx_strm = StringIO.StringIO()
  bound_iter = functools.partial(genome_alignment_iterator,
                                 reference_species=ref_spec)
  hash_func = JustInTimeGenomeAlignmentBlock.build_hash
  idx = IndexedFile(maf_strm, bound_iter, hash_func)
  idx.write_index(idx_strm)
  idx_strm.seek(0)  # seek to the start
  return idx_strm
1110501
def align_two_alignments(aln1, aln2, params=None):
    
    if not params:
        params = {}

    #create SequenceCollection object from aln1
    aln1_collection = SequenceCollection(aln1)
    #Create mapping between abbreviated IDs and full IDs
    aln1_int_map, aln1_int_keys = aln1_collection.getIntMap(prefix='aln1_')
    #Create SequenceCollection from int_map.
    aln1_int_map = SequenceCollection(aln1_int_map)

    #create SequenceCollection object from aln2
    aln2_collection = SequenceCollection(aln2)
    #Create mapping between abbreviated IDs and full IDs
    aln2_int_map, aln2_int_keys = aln2_collection.getIntMap(prefix='aln2_')
    #Create SequenceCollection from int_map.
    aln2_int_map = SequenceCollection(aln2_int_map)

    #set output and profile options
    params.update({'-out':get_tmp_filename(), '-profile':True})

    #save aln1 to tmp file
    aln1_filename = get_tmp_filename()
    aln1_out = open(aln1_filename,'w')
    aln1_out.write(aln1_int_map.toFasta())
    aln1_out.close()

    #save aln2 to tmp file
    aln2_filename = get_tmp_filename()
    aln2_out = open(aln2_filename, 'w')
    aln2_out.write(aln2_int_map.toFasta())
    aln2_out.close()

    #Create Muscle app and get results
    app = Muscle(InputHandler='_input_as_multifile', params=params,
                 WorkingDir=tempfile.gettempdir())
    res = app((aln1_filename, aln2_filename))

    #Get alignment as dict out of results
    alignment = dict(parse_fasta(res['MuscleOut']))

    #Make new dict mapping original IDs
    new_alignment = {}
    for k,v in alignment.items():
        if k in aln1_int_keys:
            new_alignment[aln1_int_keys[k]] = v
        else:
            new_alignment[aln2_int_keys[k]] = v

    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment)

    #Clean up
    res.cleanUp()
    del(aln1_collection, aln1_int_map, aln1_int_keys)
    del(aln2_collection, aln2_int_map, aln2_int_keys)
    del(app, res, alignment, params)
    remove(aln1_filename)
    remove(aln2_filename)

    return new_alignment
395665
def align(self, alignment_tool = 'clustalw', gap_opening_penalty = 0.2, ignore_bad_chains = False):
        ''''''
        if len(self.pdb_chains) > 1:
            sa = SequenceAligner(alignment_tool = alignment_tool, gap_opening_penalty = gap_opening_penalty)
            for pdb_chain in self.pdb_chains:
                sa.add_sequence('%s_%s' % (pdb_chain['pdb_id'], pdb_chain['chain_id']), pdb_chain['sequence'], ignore_bad_chains = ignore_bad_chains)
            best_matches = sa.align()
            return sa.alignment_output, best_matches
        else:
            raise Exception('Cannot align sequences - less than two chains were specified.')
785335
def check_structure_chain_quality(self, seqprop, structprop, chain_id,
                                      seq_ident_cutoff=0.5, allow_missing_on_termini=0.2,
                                      allow_mutants=True, allow_deletions=False,
                                      allow_insertions=False, allow_unresolved=True):
        """"""
        alignment = self._get_seqprop_to_structprop_alignment(seqprop=seqprop, structprop=structprop, chain_id=chain_id)

        # Compare sequence to structure's sequence using the alignment
        chain_passes_quality_check = ssbio.protein.structure.properties.quality.sequence_checker(reference_seq_aln=alignment[0],
                                                                   structure_seq_aln=alignment[1],
                                                                   seq_ident_cutoff=seq_ident_cutoff,
                                                                   allow_missing_on_termini=allow_missing_on_termini,
                                                                   allow_mutants=allow_mutants,
                                                                   allow_deletions=allow_deletions,
                                                                   allow_insertions=allow_insertions,
                                                                   allow_unresolved=allow_unresolved)
        return chain_passes_quality_check
457146
def build_tree_from_alignment(aln, moltype=DNA, best_tree=False, params=None):
    
    if params is None:
        params = {}

    if moltype == DNA or moltype == RNA:
        params['-nt'] = True
    elif moltype == PROTEIN:
        params['-nt'] = False
    else:
        raise ValueError, \
                "FastTree does not support moltype: %s" % moltype.label

    if best_tree:
        params['-slow'] = True

    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = aln.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)

    app = FastTree(params=params)

    result = app(int_map.toFasta())
    tree = DndParser(result['Tree'].read(), constructor=PhyloNode)
    #remap tip names
    for tip in tree.tips():
        tip.Name = int_keys[tip.Name]

    return tree
294111
def _reader_thread_proc(self, timeout):
    """"""
    for data in self.stream.read_until_close(timeout_ms=timeout):
      if self.stdout is not None:
        self.stdout.write(data)