<filename>vireoSNP/utils/io_utils.py

import subprocess
import numpy as np
from scipy.io import mmread
from itertools import permutations

from .vireo_base import match
from .vcf_utils import load_VCF, write_VCF, parse_donor_GPb
from .vcf_utils import read_sparse_GeneINFO, GenoINFO_maker


def match_donor_VCF(cell_dat, donor_vcf):
    """Match variants between cell VCF and donor VCF information
    """
    mm_idx = match(cell_dat['variants'], donor_vcf['variants'])
    mm_idx = mm_idx.astype(float)
    if np.sum(mm_idx == mm_idx) == 0 or np.sum(mm_idx >= 0) == 0:
        _cell_id = ["chr" + x for x in cell_dat['variants']]
        mm_idx = match(_cell_id, donor_vcf['variants'])
    if np.sum(mm_idx == mm_idx) == 0 or np.sum(mm_idx >= 0) == 0:
        _donor_id = ["chr" + x for x in donor_vcf['variants']]
        mm_idx = match(cell_dat['variants'], _donor_id)

    idx1 = np.where(mm_idx == mm_idx)[0] #remove None
    # TODO: check when chr is not compatible! given warning.
    if len(idx1) == 0:
        print("[vireo] warning: no variants matched to donor VCF, " + 
                "please check chr format!")
    else:
        print("[vireo] %d out %d variants matched to donor VCF" 
                %(len(idx1), len(cell_dat['variants'])))
    idx2 = mm_idx[idx1].astype(int)

    cell_dat['AD'] = cell_dat['AD'][idx1, :]
    cell_dat['DP'] = cell_dat['DP'][idx1, :]
    cell_dat["variants"]  = [cell_dat["variants"][x] for x in idx1]
    for _key in cell_dat["FixedINFO"].keys():
        cell_dat["FixedINFO"][_key] = [
            cell_dat["FixedINFO"][_key][x] for x in idx1]

    donor_vcf["variants"]  = [donor_vcf["variants"][x] for x in idx2]
    for _key in donor_vcf["FixedINFO"].keys():
        donor_vcf["FixedINFO"][_key] = [
            donor_vcf["FixedINFO"][_key][x] for x in idx2]
    for _key in donor_vcf["GenoINFO"].keys():
        donor_vcf["GenoINFO"][_key] = [
            donor_vcf["GenoINFO"][_key][x] for x in idx2]  

    return cell_dat, donor_vcf


def read_cellSNP(dir_name, layers=['AD', 'DP']):
    """Read data from the cellSNP output directory

    Parameters
    ----------
    dir_name:
        directory full path name for cellSNP output
    
    Return
    ------
    A disctionary containing AD, DP, cells and variants
    """
    cell_dat = load_VCF(dir_name + "/cellSNP.base.vcf.gz", load_sample=False,
                        biallelic_only=False)
    for _layer in layers:
        cell_dat[_layer] = mmread(dir_name + "/cellSNP.tag.%s.mtx" %(_layer)).tocsc()
    cell_dat['samples'] = np.genfromtxt(dir_name + "/cellSNP.samples.tsv", dtype=str)
    return cell_dat


def read_vartrix(alt_mtx, ref_mtx, cell_file, vcf_file=None):
    """Read data from VarTrix

    Parameters
    ----------
    alt_mtx:
        sparse matrix file for alternative alleles
    ref_mtx:
        sparse matrix file for reference alleles
    cell_file:
        file for cell barcodes, each per line
    vcf_file:
        the vcf file used for fetch variants in VarTrix
    
    Return
    ------
    A disctionary containing AD, DP, cells and optionally variants
    """
    if vcf_file is not None:
        cell_dat = load_VCF(vcf_file, load_sample=False, biallelic_only=False)
        cell_dat['variants'] = np.array(cell_dat['variants'])
    else:
        cell_dat = {}
    cell_dat['AD'] = mmread(alt_mtx).tocsc()
    cell_dat['DP'] = mmread(ref_mtx).tocsc() + cell_dat['AD']
    cell_dat['samples'] = np.genfromtxt(cell_file, dtype=str)
    return cell_dat


def write_donor_id(out_dir, donor_names, cell_names, n_vars, res_vireo):
    """
    Write the results of donor id into files.
    """
    ID_prob, doublet_prob = res_vireo['ID_prob'], res_vireo['doublet_prob']

    prob_max = np.max(ID_prob, axis=1)
    prob_doublet_out = np.max(doublet_prob, axis=1)
    donor_singlet = np.array(donor_names, "U100")[np.argmax(ID_prob, axis=1)]

    doublet_names = [",".join(x) for x in permutations(donor_names, 2)]
    donor_doublet = np.array(doublet_names, "U100")[np.argmax(doublet_prob, 
                                                              axis=1)]

    donor_ids = donor_singlet.copy()
    donor_ids[prob_max < 0.9] = "unassigned"
    donor_ids[prob_doublet_out >= 0.9] = "doublet"
    donor_ids[n_vars < 10] = "unassigned"


    ## save log file
    fid = open(out_dir + "/_log.txt", "w")
    fid.writelines("logLik: %.3e\n" %(res_vireo['LB_doublet']))
    fid.writelines("thetas: \n%s\n" %(res_vireo['theta_shapes']))
    fid.close()

    ## save summary file
    fid = open(out_dir + "/summary.tsv", "w")
    fid.writelines("Var1\tFreq\n")
    donor_ids_uniq, donor_ids_count = np.unique(donor_ids, return_counts=True)
    for i in range(len(donor_ids_uniq)):
        fid.writelines("%s\t%d\n" %(donor_ids_uniq[i], donor_ids_count[i]))
    fid.close()
    print("[vireo] final donor size:")
    print("\t".join([str(x) for x in donor_ids_uniq]))
    print("\t".join([str(x) for x in donor_ids_count]))

    ## save donor_ids file
    fid = open(out_dir + "/donor_ids.tsv", "w")
    header = ["cell", "donor_id", "prob_max", "prob_doublet", "n_vars", 
              "best_singlet", "best_doublet"]
    fid.writelines("\t".join(header) + "\n")
    for i in range(len(cell_names)):
        line = [cell_names[i], donor_ids[i], "%.2e" %prob_max[i],
                "%.2e" %prob_doublet_out[i], "%d" %n_vars[i],
                donor_singlet[i], donor_doublet[i]]
        fid.writelines("\t".join(line) + "\n") 
    fid.close()

    ## save singlet probability file
    fid = open(out_dir + "/prob_singlet.tsv", "w")
    fid.writelines("\t".join(["cell"] + donor_names) + "\n")
    for i in range(len(cell_names)):
        line = ["%.2e" %x for x in ID_prob[i, :]]
        fid.writelines("\t".join([cell_names[i]] + line) + "\n") 
    fid.close()

    ## save doublet probability file
    fid = open(out_dir + "/prob_doublet.tsv", "w")
    fid.writelines("\t".join(["cell"] + doublet_names) + "\n")
    for i in range(len(cell_names)):
        line = ["%.2e" %x for x in doublet_prob[i, :]]
        fid.writelines("\t".join([cell_names[i]] + line) + "\n")
    fid.close()

    ## save ambient RNA file
    if res_vireo['ambient_Psi'] is not None:
        fid = open(out_dir + "/prop_ambient.tsv", "w")
        fid.writelines("\t".join(["cell"] + donor_names + ['logLik_ratio']) + "\n")
        for i in range(len(cell_names)):
            line = ["%.4e" %x for x in res_vireo['ambient_Psi'][i, :]]
            line += ['%.2f' %res_vireo['Psi_LLRatio'][i]]
            fid.writelines("\t".join([cell_names[i]] + line) + "\n")
        fid.close()

    bashCommand = "gzip -f %s %s" %(out_dir + "/prob_singlet.tsv", 
        out_dir + "/prob_doublet.tsv")
    pro = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
    pro.communicate()[0]

def make_whitelists(donor_id_file, out_prefix):
    """Generate whitelist for each donor as input for umi_tools extract
    """
    table = np.genfromtxt(donor_id_file, dtype='str', delimiter='\t')[1:, :]
    table = table[table[:, 1] != 'unassigned', :]
    table = table[table[:, 1] != 'doublet', :]
    
    for _donor in np.unique(table[:, 1]):
        idx = table[:, 1] == _donor
        barcodes = table[idx, 0]
        fid = open(out_prefix + "_%s.txt" %_donor, "w")
        for _line in barcodes:
            fid.writelines(_line.split('-')[0] + '\n')
        fid.close()

<filename>vireoSNP/utils/io_utils.py

import subprocess
import numpy as np
from scipy.io import mmread
from itertools import permutations

from .vireo_base import match
from .vcf_utils import load_VCF, write_VCF, parse_donor_GPb
from .vcf_utils import read_sparse_GeneINFO, GenoINFO_maker


def match_donor_VCF(cell_dat, donor_vcf):
    """Match variants between cell VCF and donor VCF information
    """
    mm_idx = match(cell_dat['variants'], donor_vcf['variants'])
    mm_idx = mm_idx.astype(float)
    if np.sum(mm_idx == mm_idx) == 0 or np.sum(mm_idx >= 0) == 0:
        _cell_id = ["chr" + x for x in cell_dat['variants']]
        mm_idx = match(_cell_id, donor_vcf['variants'])
    if np.sum(mm_idx == mm_idx) == 0 or np.sum(mm_idx >= 0) == 0:
        _donor_id = ["chr" + x for x in donor_vcf['variants']]
        mm_idx = match(cell_dat['variants'], _donor_id)

    idx1 = np.where(mm_idx == mm_idx)[0] #remove None
    # TODO: check when chr is not compatible! given warning.
    if len(idx1) == 0:
        print("[vireo] warning: no variants matched to donor VCF, " + 
                "please check chr format!")
    else:
        print("[vireo] %d out %d variants matched to donor VCF" 
                %(len(idx1), len(cell_dat['variants'])))
    idx2 = mm_idx[idx1].astype(int)

    cell_dat['AD'] = cell_dat['AD'][idx1, :]
    cell_dat['DP'] = cell_dat['DP'][idx1, :]
    cell_dat["variants"]  = [cell_dat["variants"][x] for x in idx1]
    for _key in cell_dat["FixedINFO"].keys():
        cell_dat["FixedINFO"][_key] = [
            cell_dat["FixedINFO"][_key][x] for x in idx1]

    donor_vcf["variants"]  = [donor_vcf["variants"][x] for x in idx2]
    for _key in donor_vcf["FixedINFO"].keys():
        donor_vcf["FixedINFO"][_key] = [
            donor_vcf["FixedINFO"][_key][x] for x in idx2]
    for _key in donor_vcf["GenoINFO"].keys():
        donor_vcf["GenoINFO"][_key] = [
            donor_vcf["GenoINFO"][_key][x] for x in idx2]  

    return cell_dat, donor_vcf


def read_cellSNP(dir_name, layers=['AD', 'DP']):
    """Read data from the cellSNP output directory

    Parameters
    ----------
    dir_name:
        directory full path name for cellSNP output
    
    Return
    ------
    A disctionary containing AD, DP, cells and variants
    """
    cell_dat = load_VCF(dir_name + "/cellSNP.base.vcf.gz", load_sample=False,
                        biallelic_only=False)
    for _layer in layers:
        cell_dat[_layer] = mmread(dir_name + "/cellSNP.tag.%s.mtx" %(_layer)).tocsc()
    cell_dat['samples'] = np.genfromtxt(dir_name + "/cellSNP.samples.tsv", dtype=str)
    return cell_dat


def read_vartrix(alt_mtx, ref_mtx, cell_file, vcf_file=None):
    """Read data from VarTrix

    Parameters
    ----------
    alt_mtx:
        sparse matrix file for alternative alleles
    ref_mtx:
        sparse matrix file for reference alleles
    cell_file:
        file for cell barcodes, each per line
    vcf_file:
        the vcf file used for fetch variants in VarTrix
    
    Return
    ------
    A disctionary containing AD, DP, cells and optionally variants
    """
    if vcf_file is not None:
        cell_dat = load_VCF(vcf_file, load_sample=False, biallelic_only=False)
        cell_dat['variants'] = np.array(cell_dat['variants'])
    else:
        cell_dat = {}
    cell_dat['AD'] = mmread(alt_mtx).tocsc()
    cell_dat['DP'] = mmread(ref_mtx).tocsc() + cell_dat['AD']
    cell_dat['samples'] = np.genfromtxt(cell_file, dtype=str)
    return cell_dat


def write_donor_id(out_dir, donor_names, cell_names, n_vars, res_vireo):
    """
    Write the results of donor id into files.
    """
    ID_prob, doublet_prob = res_vireo['ID_prob'], res_vireo['doublet_prob']

    prob_max = np.max(ID_prob, axis=1)
    prob_doublet_out = np.max(doublet_prob, axis=1)
    donor_singlet = np.array(donor_names, "U100")[np.argmax(ID_prob, axis=1)]

    doublet_names = [",".join(x) for x in permutations(donor_names, 2)]
    donor_doublet = np.array(doublet_names, "U100")[np.argmax(doublet_prob, 
                                                              axis=1)]

    donor_ids = donor_singlet.copy()
    donor_ids[prob_max < 0.9] = "unassigned"
    donor_ids[prob_doublet_out >= 0.9] = "doublet"
    donor_ids[n_vars < 10] = "unassigned"


    ## save log file
    fid = open(out_dir + "/_log.txt", "w")
    fid.writelines("logLik: %.3e\n" %(res_vireo['LB_doublet']))
    fid.writelines("thetas: \n%s\n" %(res_vireo['theta_shapes']))
    fid.close()

    ## save summary file
    fid = open(out_dir + "/summary.tsv", "w")
    fid.writelines("Var1\tFreq\n")
    donor_ids_uniq, donor_ids_count = np.unique(donor_ids, return_counts=True)
    for i in range(len(donor_ids_uniq)):
        fid.writelines("%s\t%d\n" %(donor_ids_uniq[i], donor_ids_count[i]))
    fid.close()
    print("[vireo] final donor size:")
    print("\t".join([str(x) for x in donor_ids_uniq]))
    print("\t".join([str(x) for x in donor_ids_count]))

    ## save donor_ids file
    fid = open(out_dir + "/donor_ids.tsv", "w")
    header = ["cell", "donor_id", "prob_max", "prob_doublet", "n_vars", 
              "best_singlet", "best_doublet"]
    fid.writelines("\t".join(header) + "\n")
    for i in range(len(cell_names)):
        line = [cell_names[i], donor_ids[i], "%.2e" %prob_max[i],
                "%.2e" %prob_doublet_out[i], "%d" %n_vars[i],
                donor_singlet[i], donor_doublet[i]]
        fid.writelines("\t".join(line) + "\n") 
    fid.close()

    ## save singlet probability file
    fid = open(out_dir + "/prob_singlet.tsv", "w")
    fid.writelines("\t".join(["cell"] + donor_names) + "\n")
    for i in range(len(cell_names)):
        line = ["%.2e" %x for x in ID_prob[i, :]]
        fid.writelines("\t".join([cell_names[i]] + line) + "\n") 
    fid.close()

    ## save doublet probability file
    fid = open(out_dir + "/prob_doublet.tsv", "w")
    fid.writelines("\t".join(["cell"] + doublet_names) + "\n")
    for i in range(len(cell_names)):
        line = ["%.2e" %x for x in doublet_prob[i, :]]
        fid.writelines("\t".join([cell_names[i]] + line) + "\n")
    fid.close()

    ## save ambient RNA file
    if res_vireo['ambient_Psi'] is not None:
        fid = open(out_dir + "/prop_ambient.tsv", "w")
        fid.writelines("\t".join(["cell"] + donor_names + ['logLik_ratio']) + "\n")
        for i in range(len(cell_names)):
            line = ["%.4e" %x for x in res_vireo['ambient_Psi'][i, :]]
            line += ['%.2f' %res_vireo['Psi_LLRatio'][i]]
            fid.writelines("\t".join([cell_names[i]] + line) + "\n")
        fid.close()

    bashCommand = "gzip -f %s %s" %(out_dir + "/prob_singlet.tsv", 
        out_dir + "/prob_doublet.tsv")
    pro = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
    pro.communicate()[0]

def make_whitelists(donor_id_file, out_prefix):
    """Generate whitelist for each donor as input for umi_tools extract
    """
    table = np.genfromtxt(donor_id_file, dtype='str', delimiter='\t')[1:, :]
    table = table[table[:, 1] != 'unassigned', :]
    table = table[table[:, 1] != 'doublet', :]
    
    for _donor in np.unique(table[:, 1]):
        idx = table[:, 1] == _donor
        barcodes = table[idx, 0]
        fid = open(out_prefix + "_%s.txt" %_donor, "w")
        for _line in barcodes:
            fid.writelines(_line.split('-')[0] + '\n')
        fid.close()

Match variants between cell VCF and donor VCF information #remove None # TODO: check when chr is not compatible! given warning. Read data from the cellSNP output directory      Parameters     ----------     dir_name:         directory full path name for cellSNP output          Return     ------     A disctionary containing AD, DP, cells and variants Read data from VarTrix      Parameters     ----------     alt_mtx:         sparse matrix file for alternative alleles     ref_mtx:         sparse matrix file for reference alleles     cell_file:         file for cell barcodes, each per line     vcf_file:         the vcf file used for fetch variants in VarTrix          Return     ------     A disctionary containing AD, DP, cells and optionally variants Write the results of donor id into files. ## save log file ## save summary file ## save donor_ids file ## save singlet probability file ## save doublet probability file ## save ambient RNA file Generate whitelist for each donor as input for umi_tools extract