Convert dropEst velocyto output to anndata

dropEst2adata.py

#!/usr/bin/env python
import sys
import warnings
warnings.filterwarnings("ignore")
import logging
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
from anndata2ri.scipy2ri._r2py import rmat_to_spmat
import anndata as ann
import pandas as pd
# Include standard modules
import argparse

help = 'Convert mtx derived from dropEst to anndata objects.'

parser = argparse.ArgumentParser(description=help)
parser.add_argument("--file", "-f", help="dropEst rds containing matrices for spliced, unspliced and spanning.")
parser.add_argument("--output", "-o", help="output path and file name (wihout extension).")
parser.add_argument("--h5ad", help="save the anndata file as h5ad (default is loom).", action='store_true')
logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.INFO)


def dropEst2adata(rds):
    readRDS = robjects.r['readRDS']

    logging.info("loading rds file")
    ccm=readRDS(rds)
    colnames = robjects.r("colnames")
    rownames = robjects.r("rownames")
    names = robjects.r("names")

    def make_adata(mat):
        return ann.AnnData(rmat_to_spmat(ccm.rx2(mat)).T,
               obs=pd.DataFrame(index=list(colnames(ccm.rx2(mat)))),
               var=pd.DataFrame(index=list(rownames(ccm.rx2(mat)))))

    logging.info("generate anndata for each matrices")
    adatas=[make_adata(mat) for mat in list(names(ccm))]
    logging.info("concatenating")
    adata=ann.AnnData.concatenate(*adatas,join="outer")
    adata=ann.AnnData(adata[adatas[0].obs_names+"-0"].X,
               obs=adatas[0].obs,
               var=adata.var,
               layers={"spliced":adata[adatas[0].obs_names+"-0"].X,
                       "unspliced":adata[adatas[0].obs_names+"-1"].X,
                       "ambiguous":adata[adatas[0].obs_names+"-2"].X})

    return adata


if __name__ == "__main__":
    args = parser.parse_args()
    if args.file is None:
        parser.print_help(sys.stderr)
        sys.exit(1)
    if args.output is not None:
        out = args.output
    else:
        out = "dropest_anndata"

    adata = dropEst2adata(args.file)
    if args.h5ad:
        logging.info("writing to %s.h5ad"%out)
        adata.write_h5ad(out+".h5ad")
    else:
        logging.info("writing %s.loom"%out)
        adata.write_loom(out+".loom")

    print(adata)