count reads mapping to genes using featureCounts
    http://subread.sourceforge.net
4817
def count(data):
    
    in_bam = dd.get_work_bam(data) or dd.get_align_bam(data)
    out_dir = os.path.join(dd.get_work_dir(data), "align", dd.get_sample_name(data))
    if dd.get_aligner(data) == "star":
        out_dir = os.path.join(out_dir, "%s_%s" % (dd.get_sample_name(data), dd.get_aligner(data)))
    sorted_bam = bam.sort(in_bam, dd.get_config(data), order="queryname", out_dir=safe_makedir(out_dir))
    gtf_file = dd.get_gtf_file(data)
    work_dir = dd.get_work_dir(data)
    out_dir = os.path.join(work_dir, "htseq-count")
    safe_makedir(out_dir)
    count_file = os.path.join(out_dir, dd.get_sample_name(data)) + ".counts"
    summary_file = os.path.join(out_dir, dd.get_sample_name(data)) + ".counts.summary"
    if file_exists(count_file):
        return count_file

    featureCounts = config_utils.get_program("featureCounts", dd.get_config(data))
    paired_flag = _paired_flag(in_bam)
    strand_flag = _strand_flag(data)

    filtered_bam = bam.filter_primary(sorted_bam, data)

    cmd = ("{featureCounts} -a {gtf_file} -o {tx_count_file} -s {strand_flag} "
           "{paired_flag} {filtered_bam}")

    message = ("Count reads in {tx_count_file} mapping to {gtf_file} using "
               "featureCounts")
    with file_transaction(data, [count_file, summary_file]) as tx_files:
        tx_count_file, tx_summary_file = tx_files
        do.run(cmd.format(**locals()), message.format(**locals()))
    fixed_count_file = _format_count_file(count_file, data)
    fixed_summary_file = _change_sample_name(
        summary_file, dd.get_sample_name(data), data=data)
    shutil.move(fixed_count_file, count_file)
    shutil.move(fixed_summary_file, summary_file)

    return count_file
780980
def iTable2GFF(iTable, fa, contig = False):
    
    columns = ['#seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame', 'attribute']
    gff = {c:[] for c in columns}
    for insertion in iTable.iterrows():
        insertion = insertion[1]
        if insertion['ID'] not in fa:
            continue
        # rRNA strand
        strand = insertion['sequence'].split('strand=', 1)[1].split()[0]
        # set rRNA positions for reporting features on contig or extracted sequence
        if contig is True:
            gene = [int(i) for i in insertion['sequence'].split('pos=', 1)[1].split()[0].split('-')]
            if strand == '-':
                offset = -1 * (gene[1])
            else:
                offset = gene[0]
        else:
            strand = '+'
            gene = [1, int(insertion['sequence'].split('total-len=', 1)[1].split()[0])]
            offset = gene[0]
        insertion['strand'] = strand
        insertion['offset'] = offset
        # source for prediction
        source = insertion['sequence'].split('::model', 1)[0].rsplit(' ', 1)[-1]
        insertion['source'] = source
        # rRNA gene
        geneAnnot = '%s rRNA gene' % (source.split('from', 1)[0])
        geneNum = insertion['sequence'].split('seq=', 1)[1].split()[0]
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(source)
        gff['feature'].append('Gene')
        gff['start'].append(gene[0])
        gff['end'].append(gene[1])
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (geneNum, geneAnnot))
        # rRNA
        gff = parse_rRNA(insertion, fa[insertion['ID']], gff)
        # insertions
        gff = parse_insertion(insertion, gff)
        # orfs
        gff = parse_orf(insertion, gff)
        # catalytic RNAs
        gff = parse_catalytic(insertion, gff)
    return pd.DataFrame(gff)[columns].drop_duplicates()
1053360
def annotated_reads(txome,sorted_read_stream,output,args):
   if not os.path.exists(args.tempdir+'/'+output):
      os.makedirs(args.tempdir+'/'+output)
   ts = OrderedStream(iter(txome.transcript_stream()))
   rs = OrderedStream(GPDStream(sorted_read_stream))
   mls = MultiLocusStream([ts,rs])
   
   aof = gzip.open(args.tempdir+'/'+output+'/annotated.txt.gz','w')
   of1 = gzip.open(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz','w')
   of2 = gzip.open(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz','w')
   of3 = gzip.open(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz','w')
   z = 0
   for ml in mls:
      refs, reads = ml.payload
      if len(refs) == 0: continue
      """Check and see if we have single exons annotated"""
      single_exon_refs = [x for x in refs if x.get_exon_count()==1]
      single_exon_reads = [x for x in reads if x.get_exon_count()==1]
      #print str(len(single_exon_refs))+" "+str(len(single_exon_reads))
      unannotated_single_exon_reads = []
      for seread in single_exon_reads:
         ovs  = [(x.exons[0].length,seread.exons[0].length,x.exons[0].overlap_size(seread.exons[0]),x) for x in single_exon_refs if x.exons[0].overlaps(seread.exons[0])]
         #check for minimum overlap
         ovs = [x for x in ovs if x[2] >= args.single_exon_minimum_overlap]
         ovs = [x for x in ovs if float(x[2])/float(max(x[0],x[1])) >= args.single_exon_mutual_overlap]
         ovs = sorted(ovs,key=lambda x: float(x[2])/float(max(x[0],x[1])),reverse=True)
         if len(ovs) == 0: 
            unannotated_single_exon_reads.append(seread)
            continue
         best_ref = ovs[0][3]
         aof.write(seread.name+"\t"+best_ref.name+"\tSE"+"\n")

      """May want an optional check for any better matches among exons"""
      #now we can look for best matches among multi-exon transcripts
      reads = [x for x in reads if x.get_exon_count() > 1]
      multiexon_refs = [x for x in refs if x.get_exon_count() > 1]
      unannotated_multi_exon_reads = []
      partial_annotated_multi_exon_reads = []
      for read in reads:
         # we dont' need to have multiple exons matched. one is enough to call
         ovs = [y for y in [(x,x.exon_overlap(read,
                                              multi_minover=args.multi_exon_minimum_overlap,
                                              multi_endfrac=args.multi_exon_end_frac,
                                              multi_midfrac=args.multi_exon_mid_frac,
                                              multi_consec=False)) for x in multiexon_refs] if y[1]]
         for o in ovs: o[1].analyze_overs()
         full =  sorted([x for x in ovs if x[1].is_subset()==1],
                        key = lambda y: float(y[1].overlap_size())/float(max(y[1].tx_obj1.length,y[1].tx_obj2.length)),
                        reverse=True
                       )
         if len(full) > 0:
            aof.write(read.name+"\t"+full[0][0].name+"\tfull"+"\n")
            continue
         subset =  sorted([x for x in ovs if x[1].is_subset()==2],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(subset) > 0:
            aof.write(read.name+"\t"+subset[0][0].name+"\tpartial"+"\n")
            continue
         #check for supersets
         superset = sorted([x for x in ovs if x[1].is_subset()==3],
                          key = lambda y: (y[1].match_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         if len(superset) > 0:
            partial_annotated_multi_exon_reads.append((read,superset[0][0]))
            #print read.name+"\t"+superset[0][0].name+"\tsuper"
            continue
         #check for noncompatible overlaps
         overset = sorted([x for x in ovs if x[1].match_exon_count > 0],
                          key = lambda y: (y[1].consecutive_exon_count(),
                                           y[1].min_overlap_fraction()),
                          reverse = True
                         )
         #print [(x[1].consecutive_exon_count(), x[1].min_overlap_fraction()) for x in overset]
         if len(overset) > 0:
            partial_annotated_multi_exon_reads.append((read,overset[0][0]))
            #print read.name+"\t"+overset[0][0].name+"\tover"
            continue
         unannotated_multi_exon_reads.append(read)
      """Now we have partially annotated multi and unannotated multi and unannotated single"""
      if len(unannotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(unannotated_multi_exon_reads))+"   \r")
         d = Deconvolution(downsample(unannotated_multi_exon_reads,args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample)
         for tx in groups:
            z+=1
            of1.write(tx.get_gpd_line()+"\n")
      if len(partial_annotated_multi_exon_reads) > 0:
         sys.stderr.write(str(z)+" "+str(len(partial_annotated_multi_exon_reads))+"   \r")
         ### set the direction of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_strand(v[1].direction)
         ### set the gene name of the transcript
         for v in partial_annotated_multi_exon_reads:
            v[0].set_gene_name(v[1].gene_name)         
         d = Deconvolution(downsample([x[0] for x in partial_annotated_multi_exon_reads],args.downsample_locus))
         groups = d.parse(tolerance=20,downsample=args.downsample,use_gene_names=True)
         for tx in groups:
            z += 1
            of2.write(tx.get_gpd_line()+"\n")
      """Do the unannotated single exon reads"""
      g = Graph()
      for r in unannotated_single_exon_reads:
         if len([x for x in partial_annotated_multi_exon_reads if x[0].overlaps(r)]) > 0: continue
         if len([x for x in unannotated_multi_exon_reads if x.overlaps(r)]) > 0: continue
         if len([x for x in txome.transcripts if x.overlaps(r)]) > 0: continue
         g.add_node(Node(r))
      for i in range(0,len(g.nodes)):
         for j in range(0,len(g.nodes)):
            if i == j: continue
            if g.nodes[i].payload.overlaps(g.nodes[j].payload):
               g.add_edge(Edge(g.nodes[i],g.nodes[j]))
      g.merge_cycles()
      for r in g.roots:
         se = []
         se += r.payload_list
         children = g.get_children(r)
         for child in children:  se += child.payload_list
         rng = GenomicRange(se[0].exons[0].chr,
                            min([x.exons[0].start for x in se]),
                            max([x.exons[0].end for x in se]))
         tx = Transcript([rng],Transcript.Options(direction='+'))
         of3.write(tx.get_gpd_line()+"\n")

   of1.close()
   of2.close()
   of3.close()
   sys.stderr.write("\n")      
   """"""
   sys.stderr.write("sort transcriptome\n")
   sort_gpd(args.tempdir+'/'+output+'/partial_annotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/partial_annotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_multiexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_multiexon.sorted.gpd.gz',args)
   sort_gpd(args.tempdir+'/'+output+'/unannotated_singleexon.gpd.gz',args.tempdir+'/'+output+'/unannotated_singleexon.sorted.gpd.gz',args)
   """We still have the unannotated single exon reads to deal with"""
360726
def _dict_seq_locus(list_c, loci_obj, seq_obj):
    
    seqs = defaultdict(set)
    # n = len(list_c.keys())
    for c in list_c.values():
        for l in c.loci2seq:
            [seqs[s].add(c.id) for s in c.loci2seq[l]]

    common = [s for s in seqs if len(seqs[s]) > 1]
    seqs_in_c = defaultdict(float)
    for c in list_c.values():
        for l in c.loci2seq:
            # total = sum([v for v in loci_obj[l].coverage.values()])
            for s in c.loci2seq[l]:
                if s in common:
                    pos = seq_obj[s].pos[l]
                    # cov = 1.0 * loci_obj[l].coverage[pos] / total
                    cov = 1.0 * loci_obj[l].coverage[pos]
                    if seqs_in_c[(s, c.id)] < cov:
                        seqs_in_c[(s, c.id)] = cov
    seqs_in_c = _transform(seqs_in_c)
    return seqs_in_c
1158286
def count_genes(model):
    """"""
    genes = set()
    for reaction in model.reactions:
        if reaction.genes is None:
            continue

        if isinstance(reaction.genes, boolean.Expression):
            genes.update(v.symbol for v in reaction.genes.variables)
        else:
            genes.update(reaction.genes)

    return len(genes)
314743
def reads(args):
    
    p = OptionParser(reads.__doc__)
    p.add_option("-p", dest="prefix_length", default=4, type="int",
            help="group the reads based on the first N chars [default: %default]")
    opts, args = p.parse_args(args)

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

    frgscffile, = args
    prefix_length = opts.prefix_length

    fp = open(frgscffile)
    keyfn = lambda: defaultdict(int)
    counts = defaultdict(keyfn)
    for row in fp:
        f = FrgScfLine(row)
        fi = f.fragmentID[:prefix_length]
        counts[f.scaffoldID][fi] += 1

    for scf, count in sorted(counts.items()):
        print("{0}\t{1}".format(scf,
                ", ".join("{0}:{1}".format(*x) for x in sorted(count.items()))))
784287
def _process_feature_pub(self, limit):
        

        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph

        raw = '/'.join((self.rawdir, 'feature_pub'))
        LOG.info("processing feature_pub")

        line_counter = 0

        with open(raw, 'r') as f:
            filereader = csv.reader(f, delimiter='\t', quotechar='\"')
            f.readline()  # read the header row; skip
            for line in filereader:
                (feature_pub_id, feature_id, pub_id) = line
                # 1440    3175682 62137
                # 2       3160606 99159

                feature_key = feature_id
                if self.test_mode and not (
                        int(feature_key)
                        in self.test_keys['gene'] + self.test_keys['allele'] and
                        int(pub_id) in self.test_keys['pub']):
                    continue
                if feature_key not in self.idhash['feature']:
                    continue
                feature_id = self.idhash['feature'][feature_key]
                pub_key = pub_id
                pub_id = self.idhash['publication'][pub_key]

                graph.addTriple(pub_id, self.globaltt['mentions'], feature_id)

                line_counter += 1

                if not self.test_mode and limit is not None and line_counter > limit:
                    break

        return
876452
def find_features(seqs, locus_tag="all", utr_len=200):
    """"""
    found_features = []

    for seq_i in seqs:
        for feature in seq_i.features:
            if feature.type == "CDS" and (locus_tag == "all" or \
                    ('locus_tag' in feature.qualifiers and \
                    feature.qualifiers['locus_tag'][0] == locus_tag)):
                start = max(0, feature.location.nofuzzy_start - utr_len)
                stop  = max(0, feature.location.nofuzzy_end + utr_len)
                feature_seq = seq_i.seq[start:stop]
                f_match = FeatureMatch(feature, feature_seq, feature.strand,
                                       utr_len)
                found_features.append(f_match)

    return found_features
371118
static public void addGeneIDGFF2GeneFeatures(LinkedHashMap<String, ChromosomeSequence> chromosomeSequenceList, FeatureList listGenes) throws Exception {
		Collection<String> geneIds = listGenes.attributeValues("gene_id");
		for (String geneid : geneIds) {
			FeatureList gene = listGenes.selectByAttribute("gene_id", geneid);
			FeatureI geneFeature = gene.get(0);
			ChromosomeSequence seq = chromosomeSequenceList.get(geneFeature.seqname());
			geneid = geneid.replaceAll("_", ".G");
			AccessionID geneAccessionID = new AccessionID(geneid);
			GeneSequence geneSequence = null;
			Collection<String> transcriptids = gene.attributeValues("gene_id");
			for (String transcriptid : transcriptids) {
				// get all the individual features (exons, CDS regions, etc.) of this gene
				FeatureList transcriptFeature = listGenes.selectByAttribute("gene_id", transcriptid);
				transcriptid = transcriptid.replaceAll("_", ".G");




				//      String seqName = feature.seqname();
				//FeatureI startCodon = null;
				//FeatureI stopCodon = null;
				Integer startCodonBegin = null;
				Integer stopCodonEnd = null;
				//String startCodonName = "";
				//String stopCodonName = "";


				// now select only the coding regions of this gene
				FeatureList firstFeatures = transcriptFeature.selectByType("First");
				FeatureList terminalFeatures = transcriptFeature.selectByType("Terminal");
				FeatureList internalFeatures = transcriptFeature.selectByType("Internal");
				FeatureList singleFeatures = transcriptFeature.selectByType("Single");
				FeatureList cdsFeatures = new FeatureList();
				cdsFeatures.add(firstFeatures);
				cdsFeatures.add(terminalFeatures);
				cdsFeatures.add(internalFeatures);
				cdsFeatures.add(singleFeatures);
				// sort them
				cdsFeatures = cdsFeatures.sortByStart();
				Strand strand = Strand.POSITIVE;
				FeatureI feature = cdsFeatures.get(0);
				if (feature.location().isNegative()) {
					strand = Strand.NEGATIVE;
				}
				if (startCodonBegin == null) {
					FeatureI firstFeature = cdsFeatures.get(0);
					if (strand == Strand.NEGATIVE) {
						startCodonBegin = firstFeature.location().bioEnd();
					} else {
						startCodonBegin = firstFeature.location().bioStart();
					}
				}

				if (stopCodonEnd == null) {

					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					if (strand == Strand.NEGATIVE) {
						stopCodonEnd = lastFeature.location().bioStart();
					} else {
						stopCodonEnd = lastFeature.location().bioEnd();
					}
				}
				//for gtf ordering can be strand based so first is last and last is first
				if (startCodonBegin > stopCodonEnd) {
					int temp = startCodonBegin;
					startCodonBegin = stopCodonEnd;
					stopCodonEnd = temp;
				}

				AccessionID transcriptAccessionID = new AccessionID(transcriptid);
				if (geneSequence == null) {
					geneSequence = seq.addGene(geneAccessionID, startCodonBegin, stopCodonEnd, strand);
					geneSequence.setSource(((Feature) feature).source());
				} else {
					//if multiple transcripts for one gene make sure the gene is defined as the min and max start/end

					if (startCodonBegin < geneSequence.getBioBegin()) {
						geneSequence.setBioBegin(startCodonBegin);
					}
					if (stopCodonEnd > geneSequence.getBioBegin()) {
						geneSequence.setBioEnd(stopCodonEnd);
					}

				}
				TranscriptSequence transcriptSequence = geneSequence.addTranscript(transcriptAccessionID, startCodonBegin, stopCodonEnd);
				/*
				if (startCodon != null) {
					if (startCodonName == null || startCodonName.length() == 0) {
						startCodonName = transcriptid + "-start_codon-" + startCodon.location().bioStart() + "-" + startCodon.location().bioEnd();
					}
					transcriptSequence.addStartCodonSequence(new AccessionID(startCodonName), startCodon.location().bioStart(), startCodon.location().bioEnd());
				}
				if (stopCodon != null) {
					if (stopCodonName == null || stopCodonName.length() == 0) {
						stopCodonName = transcriptid + "-stop_codon-" + stopCodon.location().bioStart() + "-" + stopCodon.location().bioEnd();
					}
					transcriptSequence.addStopCodonSequence(new AccessionID(stopCodonName), stopCodon.location().bioStart(), stopCodon.location().bioEnd());
				}
				*/

				for (FeatureI cdsFeature : cdsFeatures) {
					Feature cds = (Feature) cdsFeature;
					String cdsName = cds.getAttribute("transcript_name");
					if (cdsName == null || cdsName.length() == 0) {
						cdsName = transcriptid + "-cds-" + cds.location().bioStart() + "-" + cds.location().bioEnd();
					}
					AccessionID cdsAccessionID = new AccessionID(cdsName);
					//ExonSequence exonSequence = geneSequence.addExon(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd());
					CDSSequence cdsSequence = transcriptSequence.addCDS(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd(), cds.frame());
					cdsSequence.setSequenceScore(cds.score());
				}
			}
		}

	}
136325
def _splitGenoGeneWindow(self,annotation_file=None,cis=1e4,funct='protein_coding',minSnps=1.,maxSnps=SP.inf):
        
         #1. load annotation
        assert annotation_file is not None, 'Splitter:: specify annotation file'
        try:
            f = h5py.File(annotation_file,'r')
            geneID = f['geneID'][:]
            gene_chrom = f['chrom'][:]
            gene_start = f['start'][:]
            gene_end = f['end'][:]
            gene_strand = f['strand'][:]
            gene_function = f['function'][:]
            f.close()
        except:
            print('Splitter:: format annotation file not valid')

        # if funct is not None, it has to be a list
        if funct is not None and funct!=list:   funct=[funct]

        windows = []
        nSnps   = []
        Igene   = []
        #2. calculates windows 
        for gene_i in range(geneID.shape[0]):
            if funct is not None:
                if gene_function[gene_i] not in funct:
                    Igene.append(False)
                    continue
            wnd = [gene_chrom[gene_i],gene_start[gene_i]-cis,gene_end[gene_i]+cis]
            Ir = (self.chrom==wnd[0])*(self.pos>=wnd[1])*(self.pos<=wnd[2])
            _nSnps = Ir.sum()
            if _nSnps>=minSnps and _nSnps<=maxSnps:
                windows.append(wnd)
                nSnps.append(_nSnps)
                Igene.append(True)
            else:
                Igene.append(False)
        Igene = SP.array(Igene)
        self.info['nSnps'] = SP.array(nSnps)
        self.info['geneID']        = geneID[Igene]
        self.info['gene_start']    = gene_start[Igene]
        self.info['gene_end']      = gene_end[Igene]
        self.info['gene_strand']   = gene_strand[Igene]
        self.info['gene_function'] = gene_function[Igene]
        return SP.array(windows)
1323543
def run_all2(protgroup, memornot, subsequences, base_outdir,
            protgroup_dict, protein_feathers_dir, date, errfile, impute_counts=True,
            cutoff_num_proteins=0, core_only_genes=None,
            length_filter_pid=.8, remove_correlated_feats=True,
            force_rerun_counts=False, force_rerun_percentages=False, force_rerun_pca=False):
    """"""
    import ssbio.utils

    # Need to set multiprocessing limit for scipy/numpy stuff if parallelizing anything
    import os
    os.environ['OMP_NUM_THREADS'] = '1'

    # First, filter down the protein group to the membrane/nonmembrane definition
    prots_filtered_feathers = get_protein_feather_paths(protgroup=protgroup, memornot=memornot,
                                                        protgroup_dict=protgroup_dict,
                                                        protein_feathers_dir=protein_feathers_dir,
                                                        core_only_genes=core_only_genes)
    num_proteins = len(prots_filtered_feathers)
    if num_proteins <= cutoff_num_proteins:
        return

    # Make output directories
    protscale = 'proteome_unscaled'
    outdir_d0 = ssbio.utils.make_dir(op.join(base_outdir, protscale))
    outdir_d1 = ssbio.utils.make_dir(op.join(outdir_d0, '-'.join(memornot)))
    outdir_final = ssbio.utils.make_dir(op.join(outdir_d1, '-'.join(protgroup)))

    if impute_counts:
        big_strain_counts_df = get_proteome_counts_impute_missing(prots_filtered_feathers=prots_filtered_feathers,
                                                                  outpath=op.join(outdir_final,
                                                                                  '{}-subsequence_proteome_IMP.fthr'.format(
                                                                                      date)),
                                                                  length_filter_pid=length_filter_pid,
                                                                  force_rerun=force_rerun_counts)

        big_strain_percents_df = get_proteome_percentages(counts_df=big_strain_counts_df,
                                                          outpath=op.join(outdir_final,
                                                                          '{}-subsequence_proteome_perc_IMP.fthr'.format(
                                                                                  date)),
                                                          force_rerun=force_rerun_percentages)
        pca_pickle = op.join(outdir_final, '{}-subsequence_pca.pckl'.format(date))
    # Divide by totals to get percentages in a new dataframe
    else:
        try:
            big_strain_percents_df = get_proteome_correct_percentages(prots_filtered_feathers=prots_filtered_feathers,
                                                                      outpath=op.join(outdir_final,
                                                                                      '{}-subsequence_proteome_perc_AVG.fthr'.format(
                                                                                              date)),
                                                                      length_filter_pid=length_filter_pid,
                                                                      force_rerun=force_rerun_percentages)
            pca_pickle = op.join(outdir_final, '{}-subsequence_pca_AVG.pckl'.format(date))
        except:
            with open(errfile, "a") as myfile:
                myfile.write('PERCENTAGES ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return



    if ssbio.utils.force_rerun(flag=force_rerun_pca, outfile=pca_pickle):

        # Then, get filters for rows of the loaded feathers for interested subsequences
        keep_subsequences = get_interested_subsequences(subsequences=subsequences)

        # Some numbers: number of features
        num_feats = len(big_strain_percents_df)

        # Make an unwieldy title
        big_title = 'LOC={0}; PROTGROUP={1};\n' \
                    'NUMPROTS={2}; NUMFEATS={3}'.format('-'.join(memornot),
                                                                        '-'.join(protgroup),
                                                                        num_proteins,
                                                                        num_feats)

        # Run PCA and make plots
        runner = PCAMultiROS(features_df=big_strain_percents_df, observations_df=pd.DataFrame(), plot_title=big_title)
        try:
            runner.clean_data(keep_features=keep_subsequences, remove_correlated_feats=remove_correlated_feats)
        except:
            with open(errfile, "a") as myfile:
                myfile.write(
                    'CLEAN ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return
        # try:
        runner.run_pca()
        # except:
        #     with open(errfile, "a") as myfile:
        #         myfile.write(
        #             'PCA ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
        #     return
        with open(pca_pickle, 'wb') as f:
            pickle.dump(runner, f)
    else:
        with open(pca_pickle, 'rb') as f:
            runner = pickle.load(f)
161104
def read_10xgenomics(cls, tarball_fpath: str, prefix: str,
                         use_ensembl_ids: bool = False):
        """"""

        _LOGGER.info('Reading file: %s', tarball_fpath)

        with tarfile.open(tarball_fpath, mode='r:gz') as tf:
            ti = tf.getmember('%smatrix.mtx' % prefix)
            with tf.extractfile(ti) as fh:
                mtx = scipy.io.mmread(fh)

            ti = tf.getmember('%sgenes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                i = 1
                if use_ensembl_ids:
                    i = 0
                gene_names = \
                        [row[i] for row in csv.reader(wrapper, delimiter='\t')]

            ti = tf.getmember('%sbarcodes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                barcodes = \
                        [row[0] for row in csv.reader(wrapper, delimiter='\t')]


            assert mtx.shape[0] == len(gene_names)
            assert mtx.shape[1] == len(barcodes)
        
        _LOGGER.info('Matrix dimensions: %s', str(mtx.shape))
        X = mtx.todense()
        matrix = cls(X=X, genes=gene_names, cells=barcodes)
        
        return matrix
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()
164845
def parse_genome_results(self, f):
    
    regexes = {
        'bam_file': r"bam file = (.+)",
        'total_reads': r"number of reads = ([\d,]+)",
        'mapped_reads': r"number of mapped reads = ([\d,]+)",
        'mapped_bases': r"number of mapped bases = ([\d,]+)",
        'sequenced_bases': r"number of sequenced bases = ([\d,]+)",
        'mean_insert_size': r"mean insert size = ([\d,\.]+)",
        'median_insert_size': r"median insert size = ([\d,\.]+)",
        'mean_mapping_quality': r"mean mapping quality = ([\d,\.]+)",
        'general_error_rate': r"general error rate = ([\d,\.]+)",
    }
    d = dict()
    for k, r in regexes.items():
        r_search = re.search(r, f['f'], re.MULTILINE)
        if r_search:
            try:
                d[k] = float(r_search.group(1).replace(',',''))
            except ValueError:
                d[k] = r_search.group(1)
    # Check we have an input filename
    if 'bam_file' not in d:
        log.debug("Couldn't find an input filename in genome_results file {}".format(f['fn']))
        return None

    # Get a nice sample name
    s_name = self.clean_s_name(d['bam_file'], f['root'])

    # Add to general stats table & calculate a nice % aligned
    try:
        self.general_stats_data[s_name]['total_reads'] = d['total_reads']
        self.general_stats_data[s_name]['mapped_reads'] = d['mapped_reads']
        d['percentage_aligned'] = (d['mapped_reads'] / d['total_reads'])*100
        self.general_stats_data[s_name]['percentage_aligned'] = d['percentage_aligned']
        self.general_stats_data[s_name]['general_error_rate'] = d['general_error_rate']*100
    except KeyError:
        pass

    # Save results
    if s_name in self.qualimap_bamqc_genome_results:
        log.debug("Duplicate genome results sample name found! Overwriting: {}".format(s_name))
    self.qualimap_bamqc_genome_results[s_name] = d
    self.add_data_source(f, s_name=s_name, section='genome_results')
511246
def _read_v3_10x_h5(filename):
    
    with tables.open_file(str(filename), 'r') as f:
        try:
            dsets = {}
            for node in f.walk_nodes('/matrix', 'Array'):
                dsets[node.name] = node.read()
            from scipy.sparse import csr_matrix
            M, N = dsets['shape']
            data = dsets['data']
            if dsets['data'].dtype == np.dtype('int32'):
                data = dsets['data'].view('float32')
                data[:] = dsets['data']
            matrix = csr_matrix((data, dsets['indices'], dsets['indptr']),
                                shape=(N, M))
            adata = AnnData(matrix,
                            {'obs_names': dsets['barcodes'].astype(str)},
                            {'var_names': dsets['name'].astype(str),
                             'gene_ids': dsets['id'].astype(str),
                             'feature_types': dsets['feature_type'].astype(str),
                             'genome': dsets['genome'].astype(str)})
            logg.info(t=True)
            return adata
        except KeyError:
            raise Exception('File is missing one or more required datasets.')
485728
def parse_featurecounts_report (self, f):
        

        file_names = list()
        parsed_data = dict()
        for l in f['f'].splitlines():
            thisrow = list()
            s = l.split("\t")
            if len(s) < 2:
                continue
            if s[0] == 'Status':
                for f_name in s[1:]:
                    file_names.append(f_name)
            else:
                k = s[0]
                if k not in self.featurecounts_keys:
                    self.featurecounts_keys.append(k)
                for val in s[1:]:
                    try:
                        thisrow.append(int(val))
                    except ValueError:
                        pass
            if len(thisrow) > 0:
                parsed_data[k] = thisrow
        # Check that this actually is a featureCounts file, as format and parsing is quite general
        if 'Assigned' not in parsed_data.keys():
            return None
        for idx, f_name in enumerate(file_names):

            # Clean up sample name
            s_name = self.clean_s_name(f_name, f['root'])

            # Reorganised parsed data for this sample
            # Collect total count number
            data = dict()
            data['Total'] = 0
            for k in parsed_data:
                data[k] = parsed_data[k][idx]
                data['Total'] += parsed_data[k][idx]

            # Calculate the percent aligned if we can
            try:
                data['percent_assigned'] = (float(data['Assigned'])/float(data['Total'])) * 100.0
            except (KeyError, ZeroDivisionError):
                pass

            # Add to the main dictionary
            if len(data) > 1:
                if s_name in self.featurecounts_data:
                    log.debug("Duplicate sample name found! Overwriting: {}".format(s_name))
                self.add_data_source(f, s_name)
                self.featurecounts_data[s_name] = data
668170
def mdl(ll, k, data):
    

    """
    N - no. of genes
    n - no. of cells 
    k - no. of cell types
    R - sum(Dataset) i.e. total no. of reads

    function TotCost = TotBits(N,m,p,R,C)
        # C is the cost from the cost function
        TotCost = C + (N*m + m*p)*(log(R/(N*p)));
    """
    N, m = data.shape
    cost = ll + (N*m + m*k)*(np.log(data.sum()/(N*k)))
    return cost