reading in phenotype file

    pfile   root of the file containing the phenotypes as NxP matrix
            (N=number of samples, P=number of traits)
546969
def readPhenoFile(pfile,idx=None):
    "

    usecols = None
    if idx!=None:
        """ different traits are comma-seperated """
        usecols = [int(x) for x in idx.split(',')]
        
    phenoFile = pfile+'.phe'
    assert os.path.exists(phenoFile), '%s is missing.'%phenoFile
    Y = SP.loadtxt(phenoFile,usecols=usecols)
    if (usecols is not None) and (len(usecols)==1): Y = Y[:,SP.newaxis]
        
    Y -= Y.mean(0); Y /= Y.std(0)
    return Y
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')
274104
def create_ped_file(samples, base_vcf, out_dir=None):
    
    out_file = "%s.ped" % utils.splitext_plus(base_vcf)[0]
    if out_dir:
        out_file = os.path.join(out_dir, os.path.basename(out_file))
    sample_ped_lines = {}
    header = ["#Family_ID", "Individual_ID", "Paternal_ID", "Maternal_ID", "Sex", "Phenotype", "Ethnicity"]
    for md_ped in list(set([x for x in [tz.get_in(["metadata", "ped"], data)
                                        for data in samples] if x is not None])):
        with open(md_ped) as in_handle:
            reader = csv.reader(in_handle, dialect="excel-tab")
            for parts in reader:
                if parts[0].startswith("#") and len(parts) > len(header):
                    header = header + parts[len(header):]
                else:
                    sample_ped_lines[parts[1]] = parts
    if not utils.file_exists(out_file):
        with file_transaction(samples[0], out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                want_samples = set(vcfutils.get_samples(base_vcf))
                writer = csv.writer(out_handle, dialect="excel-tab")
                writer.writerow(header)
                for data in samples:
                    ped_info = get_ped_info(data, samples)
                    sname = ped_info["individual_id"]
                    if sname in want_samples:
                        want_samples.remove(sname)
                        if sname in sample_ped_lines:
                            writer.writerow(sample_ped_lines[sname])
                        else:
                            writer.writerow([ped_info["family_id"], ped_info["individual_id"],
                                            ped_info["paternal_id"], ped_info["maternal_id"],
                                            ped_info["gender"], ped_info["affected"],
                                            ped_info["ethnicity"]])
    return out_file
1205984
def read_sparse(cls, file_path: str):
        """"""
        
        with open(file_path) as fh:
            next(fh)  # skip header line
            genes = next(fh)[1:-1].split('\t')
            cells = next(fh)[1:-1].split('\t')
            next(fh)
            m, n, nnz = [int(s) for s in next(fh)[:-1].split(' ')]
        
        t = pd.read_csv(file_path, sep=' ', skiprows=5, header=None,
                        dtype={0: np.uint32, 1: np.uint32})
        
        i = t[0].values - 1
        j = t[1].values - 1
        data = t[2].values

        assert data.size == nnz

        X = sparse.coo_matrix((data, (i,j)), shape=[m, n]).todense()

        return cls(X=X, genes=genes, cells=cells)
698179
def printMcmcSamplesToFile(self, profile, numSamples, outFileName):
        

        wmg = profile.getWmg(True)
        V = self.getInitialSample(wmg)

        # Print the number of candidates, phi, and the number of samples.
        outFile = open(outFileName, 'w')
        outFile.write("m," + str(profile.numCands) + '\n')
        outFile.write("phi," + str(self.phi) + '\n')
        outFile.write("numSamples," + str(numSamples))
        
        for i in range(0, numSamples):
            V = self.sampleGenerator.getNextSample(V)
            outFile.write("\n" + json.dumps(V))
        outFile.close()
980300
def _sample_variant_file_in_population(x):
    """
    '"""
    if "population" in x:
        a = _get_project_vcf(x)
        b = _get_variant_file(x, ("vrn_file",))
        decomposed = tz.get_in(("population", "decomposed"), x)
        if (a and b and not decomposed and len(a) > 0 and len(b) > 0 and
              vcfutils.get_samples(a[0]["path"]) == vcfutils.get_samples(b[0]["path"])):
            return True
    return False
223877
def _initGP(self):
        
        if self._inference=='GP2KronSum': 
            signalPos = sp.where(sp.arange(self.n_randEffs)!=self.noisPos)[0][0]
            gp = GP2KronSum(Y=self.Y, F=self.sample_designs, A=self.trait_designs,
                            Cg=self.trait_covars[signalPos], Cn=self.trait_covars[self.noisPos],
                            R=self.sample_covars[signalPos])
        else:
            mean = MeanKronSum(self.Y, self.sample_designs, self.trait_designs) 
            Iok = vec(~sp.isnan(mean.Y))[:,0] 
            if Iok.all():   Iok = None
            covar = SumCov(*[KronCov(self.trait_covars[i], self.sample_covars[i], Iok=Iok) for i in range(self.n_randEffs)])
            gp = GP(covar = covar, mean = mean) 
        self.gp = gp
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
294492
def p_incbin(p):
    
    try:
        fname = zxbpp.search_filename(p[2], p.lineno(2), local_first=True)
        if not fname:
            p[0] = None
            return
        with api.utils.open_file(fname, 'rb') as f:
            filecontent = f.read()
    except IOError:
        error(p.lineno(2), "cannot read file '%s'" % p[2])
        p[0] = None
        return

    p[0] = Asm(p.lineno(1), 'DEFB', filecontent)
1041065
def read_alignments(self, input_dir, file_format, header_grep=None):
        

        compression = ['', 'gz', 'bz2']

        if file_format == 'fasta':
            extensions = ['fa', 'fas', 'fasta']

        elif file_format == 'phylip':
            extensions = ['phy']

        else:
            extensions = []

        extensions = list('.'.join([x,y]) if y else x for x,y in itertools.product(extensions, compression))

        files = fileIO.glob_by_extensions(input_dir, extensions)
        files.sort(key=SORT_KEY)
        self._input_files = files
        records = []

        if self.show_progress:
            pbar = setup_progressbar("Loading files", len(files), simple_progress=True)
            pbar.start()

        for i, f in enumerate(files):
            if f.endswith('.gz') or f.endswith('.bz2'):
                fd, tmpfile = tempfile.mkstemp()
                with fileIO.freader(f, f.endswith('.gz'), f.endswith('.bz2')) as reader,\
                     fileIO.fwriter(tmpfile) as writer:
                    for line in reader:
                        if ISPY3:
                            line = line.decode()
                        writer.write(line)
                try:
                    record = Alignment(tmpfile, file_format, True)
                except ValueError:
                    record = Alignment(tmpfile, file_format, False)
                finally:
                    os.close(fd)
                    os.unlink(tmpfile)

            else:
                try:
                    record = Alignment(f, file_format, True)
                except RuntimeError:
                    record = Alignment(f, file_format, False)

            if header_grep:
                try:
                    datatype = 'dna' if record.is_dna() else 'protein'

                    record = Alignment([(header_grep(x), y) for (x, y) in record.get_sequences()], datatype)

                except TypeError:
                    raise TypeError("Couldn't apply header_grep to header\n"
                                    "alignment number={}, name={}\n"
                                    "header_grep={}".format(i, fileIO.strip_extensions(f), header_grep))
                except RuntimeError:
                    print('RuntimeError occurred processing alignment number={}, name={}'
                          .format(i, fileIO.strip_extensions(f)))
                    raise

            record.name = (fileIO.strip_extensions(f))
            records.append(record)
            if self.show_progress:
                pbar.update(i)
        if self.show_progress:
            pbar.finish()
        return records
1128962
def summary(args):
    
    from jcvi.utils.natsort import natsort_key

    p = OptionParser(summary.__doc__)
    p.add_option("--suffix", default="Mb",
            help="make the base pair counts human readable [default: %default]")
    p.add_option("--ids",
            help="write the ids that have >= 50% N's [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    idsfile = opts.ids
    header = "Seqid Real N's Total %_real".split()
    if idsfile:
        idsfile = open(idsfile, "w")
        nids = 0

    data = []
    for fastafile in args:
        for rec in SeqIO.parse(fastafile, "fasta"):
            seqlen = len(rec)
            nns = rec.seq.count('n') + rec.seq.count('N')
            reals = seqlen - nns
            pct = reals * 100. / seqlen
            pctreal = "{0:.1f}%".format(pct)
            if idsfile and pct < 50:
                nids += 1
                print(rec.id, file=idsfile)

            data.append((rec.id, reals, nns, seqlen, pctreal))

    data.sort(key=natsort_key)
    ids, reals, nns, seqlen, pctreal = zip(*data)
    reals = sum(reals)
    nns = sum(nns)
    seqlen = sum(seqlen)
    pctreal = "{0:.1f}%".format(reals * 100. / seqlen)
    data.append(("Total", reals, nns, seqlen, pctreal))

    write_csv(header, data, sep=" ", filename=opts.outfile, thousands=True)
    if idsfile:
        logging.debug("A total of {0} ids >= 50% N's written to {1}.".\
                      format(nids, idsfile.name))
        idsfile.close()

    return reals, nns, seqlen
33767
def _process_phendesc(self, limit):
        
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        model = Model(graph)
        raw = '/'.join((self.rawdir, 'phendesc'))
        LOG.info("processing G2P")

        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:
                (
                    phendesc_id,
                    genotype_id,
                    environment_id,
                    description,
                    type_id,
                    pub_id) = line
                # 1	2	1	Hemizygous males are wild type, homozygous males are sterile.	60466	209729

                line_counter += 1
                phendesc_key = phendesc_id
                phendesc_id = self._makeInternalIdentifier('phendesc', phendesc_key)

                # for now, just attach the description to the genotype
                genotype_key = genotype_id
                genotype_id = self.idhash['genotype'][genotype_key]
                pub_key = pub_id
                pub_id = self.idhash['publication'][pub_key]

                environment_key = environment_id
                environment_id = self.idhash['environment'][environment_key]

                if self.test_mode and\
                        int(genotype_key) not in self.test_keys['genotype']:
                    continue

                # TODO type id ==> ECO???

                # just make associations with abnormal phenotype
                # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                # note this is not "abnormal phenotype" if that was what is wanted
                # "phenotype": "FBcv:0001347" ~~>  "Phenotype": "UPHENO:0001001"
                # but it is a near collision with an existing term
                phenotype_id = self.localtt['phenotype']
                assoc = G2PAssoc(graph, self.name, genotype_id, phenotype_id)
                assoc.add_source(pub_id)
                assoc.set_description(description)
                assoc.set_environment(environment_id)
                assoc.add_association_to_graph()
                assoc_id = assoc.get_association_id()
                model.addComment(assoc_id, phendesc_id)

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

        return
571758
def read(self, filename, rowprefix=None, colprefix=None, delim=":"):
        """"""

        self.matrix = scipy.io.mmread(filename)
        with open(filename + ".rownames") as in_handle:
            self.rownames = [x.strip() for x in in_handle]
            if rowprefix:
                self.rownames = [rowprefix + delim + x for x in self.rownames]
        with open(filename + ".colnames") as in_handle:
            self.colnames = [x.strip() for x in in_handle]
            if colprefix:
                self.colnames = [colprefix + delim + x for x in self.colnames]
466433
def load(self, read_tuple_name):
        
        self.prefix_width = 0
        self.read_tuple_id_width = 0
        self.genome_id_width = 0
        self.chr_id_width = 0
        self.coor_width = 0

        parts = read_tuple_name.split("__")
        self.prefix_width = len(parts[0])
        self.read_tuple_id_width = len(parts[1])

        segments = parts[2][1:-1].split("),(")
        for segment in segments:
            int_widths = list(map(len, segment.split(",")))
            self.genome_id_width = max(self.genome_id_width, int_widths[0])
            self.chr_id_width = max(self.chr_id_width, int_widths[1])
            self.coor_width = max(self.coor_width, int_widths[2], int_widths[3])
394069
def profiler(self):
        
        # Initialise variables
        profiledata = defaultdict(make_dict)
        profileset = set()
        genedict = {}
        # Find all the unique profiles to use with a set
        for sample in self.metadata.samples:
            if sample[self.analysistype].profile != 'NA':
                profileset.add(sample[self.analysistype].profile[0])
        # Extract the profiles for each set
        for sequenceprofile in profileset:
            # Clear the list of genes
            genelist = []
            for sample in self.metadata.samples:
                if sequenceprofile == sample[self.analysistype].profile[0]:
                    genelist = [os.path.split(x)[1].split('.')[0] for x in sample[self.analysistype].alleles]
            try:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Revert to standard comma separated values
            except KeyError:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Iterate through the rows
            for row in profile:
                # Iterate through the genes
                for gene in genelist:
                    # Add the sequence profile, and type, the gene name and the allele number to the dictionary
                    try:
                        profiledata[sequenceprofile][row['ST']][gene] = row[gene]
                    except KeyError:
                        pass
            # Add the gene list to a dictionary
            genedict[sequenceprofile] = sorted(genelist)
            # Add the profile data, and gene list to each sample
            for sample in self.metadata.samples:
                if sample.general.bestassemblyfile != 'NA':
                    if sequenceprofile == sample[self.analysistype].profile[0]:
                        # Populate the metadata with the profile data
                        sample[self.analysistype].profiledata = profiledata[sample[self.analysistype].profile[0]]
                        # Add the allele directory to a list of directories used in this analysis
                        self.allelefolders.add(sample[self.analysistype].alleledir)
                        dotter()
638711
def read_next_maf( file, species_to_lengths=None, parse_e_rows=False ):
    
    alignment = Alignment(species_to_lengths=species_to_lengths)
    # Attributes line
    line = readline( file, skip_blank=True )
    if not line: return None
    fields = line.split() 
    if fields[0] != 'a': raise Exception("Expected 'a ...' line")
    alignment.attributes = parse_attributes( fields[1:] )
    if 'score' in alignment.attributes:
        alignment.score = alignment.attributes['score']
        del alignment.attributes['score']
    else:
        alignment.score = 0
    # Sequence lines
    last_component = None
    while 1:
        line = readline( file )
        # EOF or Blank line terminates alignment components
        if not line or line.isspace(): break
        if line.isspace(): break 
        # Parse row
        fields = line.split()
        if fields[0] == 's':
            # An 's' row contains sequence for a component
            component = Component()
            component.src = fields[1]
            component.start = int( fields[2] )
            component.size = int( fields[3] )
            component.strand = fields[4]
            component.src_size = int( fields[5] )
            if len(fields) > 6: component.text = fields[6].strip()
            # Add to set
            alignment.add_component( component )
            last_component = component
        elif fields[0] == 'e':
            # An 'e' row, when no bases align for a given species this tells
            # us something about the synteny 
            if parse_e_rows:
                component = Component()
                component.empty = True
                component.src = fields[1]
                component.start = int( fields[2] )
                component.size = int( fields[3] )
                component.strand = fields[4]
                component.src_size = int( fields[5] )
                component.text = None
                synteny = fields[6].strip()
                assert len( synteny ) == 1, \
                    "Synteny status in 'e' rows should be denoted with a single character code"
                component.synteny_empty = synteny
                alignment.add_component( component )
                last_component = component
        elif fields[0] == 'i':
            # An 'i' row, indicates left and right synteny status for the 
            # previous component, we hope ;)
            assert fields[1] == last_component.src, "'i' row does not follow matching 's' row"
            last_component.synteny_left = ( fields[2], int( fields[3] ) )
            last_component.synteny_right = ( fields[4], int( fields[5] ) )
        elif fields[0] == 'q':
            assert fields[1] == last_component.src, "'q' row does not follow matching 's' row"
            # TODO: Should convert this to an integer array?
            last_component.quality = fields[2]
            
    return alignment
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file