%prog frg frgfile

    Extract FASTA sequences from frg reads.
434753
def frg(args):
    
    p = OptionParser(frg.__doc__)
    opts, args = p.parse_args(args)

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

    frgfile, = args
    fastafile = frgfile.rsplit(".", 1)[0] + ".fasta"
    fp = open(frgfile)
    fw = open(fastafile, "w")

    for rec in iter_records(fp):
        if rec.type != "FRG":
            continue
        id = rec.get_field("acc")
        seq = rec.get_field("seq")
        s = SeqRecord(Seq(seq), id=id, description="")
        SeqIO.write([s], fw, "fasta")

    fw.close()
1233138
def summary(args):
    
    p = OptionParser(summary.__doc__)

    opts, args = p.parse_args(args)

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

    fastafile, = args
    f = Fasta(fastafile, index=False)

    halfmaskedseqs = set()
    allmasked = 0
    allbases = 0
    cutoff = 50
    for key, seq in f.iteritems():
        masked = 0
        for base in seq:
            if base not in "AGCT":
                masked += 1
        seqlen = len(seq)
        if masked * 100. / seqlen > cutoff:
            halfmaskedseqs.add(key)
        allmasked += masked
        allbases += seqlen

    seqnum = len(f)
    maskedseqnum = len(halfmaskedseqs)

    print("Total masked bases: {0}".\
            format(percentage(allmasked, allbases)), file=sys.stderr)
    print("Total masked sequences (contain > {0}% masked): {1}".\
            format(cutoff, percentage(maskedseqnum, seqnum)), file=sys.stderr)
966991
def main(args):

  
  sys.stderr.write("reading reference genome\n")
  ref = FASTAData(open(args.genome).read())
  sys.stderr.write("read in "+str(len(ref.keys()))+" chromosomes\n")

  """ Next make the transcriptome """
  txome = {}
  sys.stderr.write("write the transcriptome\n")
  inf = None
  if args.gpd[-3:] == '.gz':
    inf = gzip.open(args.gpd)
  else:
    inf = open(args.gpd)
  stream = GPDStream(inf)
  tof = open(args.tempdir+'/transcriptome.fa','w')
  z = 0
  for gpd in stream:
    z += 1
    if gpd.transcript_name in txome:
      sys.stderr.write("WARNING already have a transcript "+gpd.transcript_name+" ignoring line "+str(z)+" of the gpd\n")
      continue
    txome[gpd.transcript_name] = gpd.gene_name
    tof.write('>'+gpd.transcript_name+"\n"+str(gpd.get_sequence(ref))+"\n")
  tof.close()
  inf.close()
  sys.stderr.write("wrote "+str(len(txome.keys()))+" transcripts\n")

  """Build the salmon index"""
  sys.stderr.write("building a salmon index\n")
  cmd = 'salmon index -p '+str(args.numThreads)+' -t '+args.tempdir+'/transcriptome.fa -i '+args.tempdir+'/salmon_index'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished building the index\n")

  """Use the index to quanitfy"""
  sys.stderr.write("quanitfy reads\n")
  reads = ''
  if args.rU:
    reads = '-r '+args.rU
  else:
    reads = '-1 '+args.r1+' -2 '+args.r2
  cmd = 'salmon quant -p '+str(args.numThreads)+' -i '+args.tempdir+'/salmon_index -l A '+reads+' -o '+args.tempdir+'/output_quant'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished quanitfying\n")

  """Now parse the salmon output to add gene name"""
  salmon = {}
  with open(args.tempdir+'/output_quant/quant.sf') as inf:
    header = inf.readline()
    for line in inf:
      f = line.rstrip().split("\t")
      # by each transcript name hold a data strcture of the information
      salmon[f[0]] = {'name':f[0],'length':int(f[1]),'EffectiveLength':float(f[2]),'TPM':float(f[3]),'NumReads':float(f[4])}
  genes = {}
  for name in salmon:
    gene = txome[name]
    if gene not in genes: genes[gene] = []
    genes[gene].append(salmon[name])
  genetot = {}
  for gene in genes:
    tot = sum([x['TPM'] for x in genes[gene]])
    genetot[gene] = tot
  ordered_gene_names = sorted(genetot.keys(), key=lambda x: genetot[x],reverse=True)
      
  """Collected enough information to make output"""
  sys.stderr.write("generating output\n")
  of = sys.stdout
  if args.output != '-':
    of = open(args.output,'w')
  of.write("geneName\ttranscriptName\tlength\tEffectiveLength\ttxTPM\tNumReads\tgeneTPM\n")
  for gene in ordered_gene_names:
    txs = sorted(genes[gene],key=lambda x: x['TPM'],reverse=True)
    for tx in txs:
      of.write(gene+"\t"+tx['name']+"\t"+str(tx['length'])+"\t"+str(tx['EffectiveLength'])+"\t"+str(tx['TPM'])+"\t"+str(tx['NumReads'])+"\t"+str(genetot[gene])+"\n")
  of.close()

  # Temporary working directory step 3 of 3 - Cleanup
  if not args.specific_tempdir:
    rmtree(args.tempdir)
885254
def parse_args():
    """"""
    parser = argparse.ArgumentParser(
        description='Text Classification with FastText',
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    # Computation options
    group = parser.add_argument_group('Computation arguments')
    group.add_argument('--input', type=str, help='Input file location')
    group.add_argument(
        '--validation', type=str, help='Validation file Location ')
    group.add_argument(
        '--output', type=str, help='Location to save trained model')
    group.add_argument(
        '--ngrams', type=int, default=1, help='NGrams used for training')
    group.add_argument(
        '--batch_size', type=int, default=16, help='Batch size for training.')
    group.add_argument('--epochs', type=int, default=10, help='Epoch limit')
    group.add_argument(
        '--gpu',
        type=int,
        help=('Number (index) of GPU to run on, e.g. 0. '
              'If not specified, uses CPU.'))
    group.add_argument(
        '--no-hybridize',
        action='store_true',
        help='Disable hybridization of gluon HybridBlocks.')

    # Model
    group = parser.add_argument_group('Model arguments')
    group.add_argument(
        '--emsize', type=int, default=100, help='Size of embedding vectors.')

    # Optimization options
    group = parser.add_argument_group('Optimization arguments')
    group.add_argument('--optimizer', type=str, default='adam')
    group.add_argument('--lr', type=float, default=0.05)

    args = parser.parse_args()
    return args
720076
def read_fastas(input_files):
    
    tumor_file = [y for x, y in input_files.items() if x.startswith('T')][0]
    normal_file = [y for x, y in input_files.items() if x.startswith('N')][0]
    output_files = defaultdict(list)
    output_files = _read_fasta(tumor_file, output_files)
    num_entries = len(output_files)
    output_files = _read_fasta(normal_file, output_files)
    assert len(output_files) == num_entries
    return output_files
240081
def fasta(self):
        
        logging.info('Extracting FASTA sequences matching PointFinder database')
        for sample in self.runmetadata.samples:
            # Ensure that there are sequence data to extract from the GenObject
            if GenObject.isattr(sample[self.analysistype], 'sequences'):
                # Set the name of the FASTA file
                sample[self.analysistype].pointfinderfasta = \
                    os.path.join(sample[self.analysistype].outputdir,
                                 '{seqid}_pointfinder.fasta'.format(seqid=sample.name))
                # Create a list to store all the SeqRecords created
                sequences = list()
                with open(sample[self.analysistype].pointfinderfasta, 'w') as fasta:
                    for gene, sequence in sample[self.analysistype].sequences.items():
                        # Create a SeqRecord using a Seq() of the sequence - both SeqRecord and Seq are from BioPython
                        seq = SeqRecord(seq=Seq(sequence),
                                        id=gene,
                                        name=str(),
                                        description=str())
                        sequences.append(seq)
                    # Write all the SeqRecords to file
                    SeqIO.write(sequences, fasta, 'fasta')
779581
def _extract_rsnap_isteps(rproffile):
    """"""
    step_regex = re.compile(r'^\*+step:\s*(\d+) ; time =\s*(\S+)')
    isteps = []  # list of (istep, time, nz)
    rows_to_del = set()
    line = ' '
    with rproffile.open() as stream:
        while line[0] != '*':
            line = stream.readline()
        match = step_regex.match(line)
        istep = int(match.group(1))
        time = float(match.group(2))
        nlines = 0
        iline = 0
        for line in stream:
            if line[0] == '*':
                isteps.append((istep, time, nlines))
                match = step_regex.match(line)
                istep = int(match.group(1))
                time = float(match.group(2))
                nlines = 0
                # remove useless lines produced when run is restarted
                nrows_to_del = 0
                while isteps and istep <= isteps[-1][0]:
                    nrows_to_del += isteps.pop()[-1]
                rows_to_del = rows_to_del.union(
                    range(iline - nrows_to_del, iline))
            else:
                nlines += 1
                iline += 1
        isteps.append((istep, time, nlines))
    return isteps, rows_to_del
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
677722
def gb(args):
    
    from Bio.Alphabet import generic_dna
    try:
        from BCBio import GFF
    except ImportError:
        print("You need to install dep first: $ easy_install bcbio-gff", file=sys.stderr)

    p = OptionParser(gb.__doc__)
    opts, args = p.parse_args(args)

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

    gff_file, fasta_file = args
    pf = op.splitext(gff_file)[0]
    out_file = pf + ".gb"
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")
1242160
def parse_insertion(insertion, gff):
    
    offset = insertion['offset']
    for ins in parse_fasta(insertion['insertion sequence'].split('|')):
        strand = insertion['strand']
        ID = ins[0].split('>')[1].split()[0]
        Start, End = [int(i) for i in ins[0].split('gene-pos=', 1)[1].split()[0].split('-')]
        Start, End = abs(Start + offset), abs(End + offset)
        if strand == '-':
            Start, End = End, Start
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(insertion['source'])
        gff['feature'].append('IVS')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand) # same as rRNA
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s' % (ID))
    return gff
296321
def get_sequences(self, chr, coords):
            
        # Check if we have an index_dir
        if not self.index_dir:
            print("Index dir is not defined!")
            sys.exit()

        # retrieve all information for this specific sequence
        fasta_file = self.fasta_file[chr]
        index_file = self.index_file[chr]
        line_size = self.line_size[chr]
        total_size = self.size[chr]
        index = open(index_file, "rb")
        fasta = open(fasta_file)
        
        seqs = []
        for coordset in coords:
            seq = ""
            for (start,end) in coordset: 
                if start > total_size:
                    raise ValueError("%s: %s, invalid start, greater than sequence length!" % (chr,start))
            
                if start < 0:
                    raise ValueError("Invalid start, < 0!")
                
                if end > total_size:
                    raise ValueError("Invalid end, greater than sequence length!")


                seq += self._read(index, fasta, start, end, line_size)
            seqs.append(seq)
        index.close()
        fasta.close()

        return seqs
1296791
def analyze_fa(fa):
    
    if fa.name == '<stdin>':
        safe = 'temp.id'
    else:
        safe = '%s.id' % (fa.name)
    safe = open(safe, 'w')
    sequences = {} # sequences[id] = sequence
    insertions = {} # insertions[id] = [[start, stop], [start, stop], ...]
    count = 0
    id2name = {}
    names = []
    for seq in parse_fasta(fa):
        id = '%010d' % (count,)
        name = seq[0].split('>', 1)[1]
        id2name[id] = name
        id2name[name] = id
        names.append(name)
        insertions[id] = insertions_from_masked(seq[1])
        sequences[id] = seq
        print('\n'.join(['>%s' % (id), seq[1].upper()]), file=safe)
        count += 1
    safe.close()
    lookup = open('%s.id.lookup' % (fa.name), 'w')
    for i in list(id2name.items()):
        print('\t'.join(i), file=lookup)
    lookup.close()
    return safe.name, sequences, id2name, names, insertions
1277236
def _fasta_slice(fasta, seqid, start, stop, strand):
    
    _strand = 1 if strand == '+' else -1
    return fasta.sequence({'chr': seqid, 'start': start, 'stop': stop, \
        'strand': _strand})
682050
def read_hlas(fasta_fai):
    
    out = []
    with open(fasta_fai) as in_handle:
        for line in in_handle:
            if line.startswith("HLA"):
                out.append(line.split()[0])
    return out
154182
def needle(reads):
    
    from tempfile import mkdtemp
    from shutil import rmtree

    dir = mkdtemp()

    file1 = join(dir, 'file1.fasta')
    with open(file1, 'w') as fp:
        print(reads[0].toString('fasta'), end='', file=fp)

    file2 = join(dir, 'file2.fasta')
    with open(file2, 'w') as fp:
        print(reads[1].toString('fasta'), end='', file=fp)

    out = join(dir, 'result.fasta')

    Executor().execute("needle -asequence '%s' -bsequence '%s' -auto "
                       "-outfile '%s' -aformat fasta" % (
                           file1, file2, out))

    # Use 'list' in the following to force reading the FASTA from disk.
    result = Reads(list(FastaReads(out)))
    rmtree(dir)

    return result
680179
def sequences_from_fasta(path):
    
    from Bio import SeqIO
    return {x.description: x.seq for x in SeqIO.parse(path, 'fasta')}
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