R/cellranger.R

In RyanYip-Kat/yipCat: Some function like ArchR style in Seurat Object

#' run cellranger reanalysis(only for cellranger version 3.1.0)
#' @param crgPython path of  cellranger python,eg(cellranger-3.1.0/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c10/bin/python)
#' @param targetMatrix target matrix want to transform eg(outs/filtered_feature_bc_matrix)
#' @param id id name
#' @param path2CRG path of cellranger,only support cellranger-3.1.0
#' @examples
#' path2CRG="/home/test/Software/cellranger-3.1.0"
#' @export
generateCloupe<-function(crgPython="python",
                         targetMatrix="outs/filtered_feature_bc_matrix",
                         CRH5="filtered_feature_bc_matrix.h5",
                         id="reanalysis",
                         path2CRG="/home/ye/Software/cellranger-3.1.0"
                         ){

  outH5<-"filtered_feature_bc_matrix_gold.h5"
  .transform10XH5(crgPython=crgPython,
                  targetMatrix = targetMatrix,
                  CRH5 = CRH5,
                  outH5 = outH5,
                  path2CRG = path2CRG
                  )
  cmd<-sprintf("reanalyze --id %s --matrix %s",id,outH5)
  message("reanalysis..")
  cellranger<-file.path(path2CRG,"cellranger-cs/3.1.0/bin/cellranger")
  run <- system2(cellranger, cmd, wait=TRUE, stdout=NULL, stderr=NULL)
}

#' code to transform matrix into 10x h5 matrix
#' @export
scMatrix2CellRangerH5<-"import sys
import argparse
import numpy as np
import h5py
import shutil
import os

parser = argparse.ArgumentParser(description='Turn Seurat output into h5 file that could be re-analyzed by Cellranger.')
parser.add_argument('-m', '--matrix', required=True, help = 'Path to mtx file of Seurat output')
parser.add_argument('-c','--crh5', required=True, help = 'Path of the h5 file generated by Cellranger')
parser.add_argument('-o', '--output', required=True, help = 'H5 file to be handled.')
parser.add_argument('--cellranger_path', default = '/home/zyserver/shizhuoxing/software/cellranger-3.1.0/', help = 'Main path of Cellrager.')
args = parser.parse_args()


sys.path.append(args.cellranger_path+'/cellranger-cs/3.1.0/tenkit/lib/python/')
sys.path.append(args.cellranger_path+'/cellranger-cs/3.1.0/lib/python/')
import cellranger.matrix as cr_matrix

mtx = cr_matrix.CountMatrix.load_mtx(args.matrix)
mtx.save_h5_file(args.output)

CR_H5 = args.matrix + '/' + 'tmp.h5'
shutil.copyfile(args.crh5,CR_H5)

with h5py.File(args.output,'r') as seu, h5py.File(CR_H5,'a') as cel:
        del cel['matrix/barcodes']
        cel.create_dataset('matrix/barcodes', data=seu['matrix/barcodes'])
        del cel['matrix/data']
        cel.create_dataset('matrix/data', data=seu['matrix/data'])
        del cel['matrix/indices']
        cel.create_dataset('matrix/indices', data=seu['matrix/indices'])
        del cel['matrix/indptr']
        cel.create_dataset('matrix/indptr', data=seu['matrix/indptr'])
        del cel['matrix/shape']
        cel.create_dataset('matrix/shape', data=seu['matrix/shape'])

        ftshape = seu['matrix/features/feature_type'].shape[0]
        genome_new = np.repeat(cel['matrix/features/genome'][0],ftshape)

        del cel['matrix/features/genome']
        cel.create_dataset('matrix/features/genome', data=genome_new)
        del cel['matrix/features/feature_type']
        cel.create_dataset('matrix/features/feature_type', data=seu['matrix/features/feature_type'])
        del cel['matrix/features/id']
        cel.create_dataset('matrix/features/id', data=seu['matrix/features/id'])
        del cel['matrix/features/name']
        cel.create_dataset('matrix/features/name', data=seu['matrix/features/name'])

os.remove(args.output)
shutil.copyfile(CR_H5,args.output)
os.remove(CR_H5)
"

#' write strings into script file
#' @export
.writeCode<-function(txt=NULL,file="tmp",ext=".py"){
  stopifnot(ext%in%c(".R",".py",".sh"))
  file=paste0(file,ext)
  write(txt,file=file)
  return(file)
}

#' transform cutsom matrix into 10X h5 matrix
#'
#' @param crgPython path of  cellranger python,eg(cellranger-3.1.0/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c10/bin/python)
#' @param targetMatrix target matrix want to transform eg(outs/filtered_feature_bc_matrix)
#' @param outH5 h5 file to store output
#' @param path2CRG path of cellranger
#' @export
.transform10XH5<-function(crgPython="python",
                          targetMatrix="outs/filtered_feature_bc_matrix",
                          CRH5="filtered_feature_bc_matrix.h5",
                          outH5="filtered_feature_bc_matrix_gold.h5",
                          path2CRG="cellranger"){

  code<-.writeCode(txt = scMatrix2CellRangerH5,file = "h5tranformer",ext = ".py")
  cmd<-sprintf("%s -m %s -c %s -o %s --cellranger_path %s",
               code,targetMatrix,CRH5,outH5,path2CRG)

  message("tranform..")
  run <- system2(crgPython, cmd, wait=TRUE, stdout=NULL, stderr=NULL)
  file.remove(code)
}

#txt<-'cd %s/matrix  && awk '{print $1"\t"$2"\tGene Expression"}' genes.tsv > features.tsv && rm genes.tsv && sed -i 's/\./\-/' barcodes.tsv && gzip *'