| 24 | |
| 25 | |
| 26 | def initialFiltering( args ): |
| 27 | |
| 28 | sc.logging.print_versions() |
| 29 | |
| 30 | ################################################## |
| 31 | # read unfiltered data from a loom file |
| 32 | adata = sc.read_loom(args.loom_input) |
| 33 | |
| 34 | ################################################## |
| 35 | # basic filtering / stats |
| 36 | |
| 37 | nCountsPerGene = np.sum(adata.X, axis=0) |
| 38 | nCellsPerGene = np.sum(adata.X>0, axis=0) |
| 39 | |
| 40 | # Show info |
| 41 | print("Number of counts (in the dataset units) per gene:", nCountsPerGene.min(), " - " ,nCountsPerGene.max()) |
| 42 | print("Number of cells in which each gene is detected:", nCellsPerGene.min(), " - " ,nCellsPerGene.max()) |
| 43 | |
| 44 | nCells=adata.X.shape[0] |
| 45 | |
| 46 | # pySCENIC thresholds |
| 47 | minCountsPerGene=3*.01*nCells # 3 counts in 1% of cells |
| 48 | print("minCountsPerGene: ", minCountsPerGene) |
| 49 | |
| 50 | minSamples=.01*nCells # 1% of cells |
| 51 | print("minSamples: ", minSamples) |
| 52 | |
| 53 | #################### |
| 54 | # initial cuts |
| 55 | sc.pp.filter_cells(adata, min_genes=args.thr_min_genes) |
| 56 | sc.pp.filter_genes(adata, min_cells=args.thr_min_cells) |
| 57 | |
| 58 | #################### |
| 59 | # mito and genes/counts cuts |
| 60 | mito_genes = adata.var_names.str.startswith('MT-') |
| 61 | # for each cell compute fraction of counts in mito genes vs. all genes |
| 62 | if( sum(mito_genes)==0 ): |
| 63 | adata.obs['percent_mito'] = 0.0 |
| 64 | else: |
| 65 | adata.obs['percent_mito'] = np.ravel(np.sum(np.asmatrix(adata[:, mito_genes].X.todense()), axis=1)) / np.ravel(np.sum(adata.X, axis=1)) |
| 66 | # add the total counts per cell as observations-annotation to adata |
| 67 | adata.obs['n_counts'] = np.ravel(adata.X.sum(axis=1)) |
| 68 | |
| 69 | adata = adata[adata.obs['n_genes'] < args.thr_n_genes, :] |
| 70 | adata = adata[adata.obs['percent_mito'] < args.thr_pct_mito, :] |
| 71 | |
| 72 | ################################################## |
| 73 | # output to loom file: |
| 74 | row_attrs = { |
| 75 | "Gene": np.array(adata.var_names) , |
| 76 | } |
| 77 | col_attrs = { |
| 78 | "CellID": np.array(adata.obs_names) , |
| 79 | "nGene": np.array( np.sum(adata.X.transpose()>0 , axis=0)).flatten() , |
| 80 | "nUMI": np.array( np.sum(adata.X.transpose() , axis=0)).flatten() , |
| 81 | } |
| 82 | |
| 83 | lp.create(args.loom_filtered, adata.X.transpose(), row_attrs, col_attrs) |