Modern developments in single cell sequencing technologies enable broad insights into cellular state. Single cell RNA sequencing (scRNA-seq) can be used to explore cell types, states, and developmental trajectories to broaden understanding of cell heterogeneity in tissues and organs. Analysis of these sparse, high-dimensional experimental results requires dimension reduction. Several methods have been developed to estimate low-dimensional embeddings for filtered and normalized single cell data. However, methods have yet to be developed for unfiltered and unnormalized count data. We present a nonlinear latent variable model with robust, heavy-tailed error and adaptive kernel learning to estimate low-dimensional nonlinear structure in scRNA-seq data. Gene expression in a single cell is modeled as a noisy draw from a Gaussian process in high dimensions from low-dimensional latent positions. This model is called the Gaussian process latent variable model (GPLVM). We model residual errors with a heavy-tailed Student's t-distribution to estimate a manifold that is robust to technical and biological noise. We compare our approach to common dimension reduction tools to highlight our model's ability to enable important downstream tasks, including clustering and inferring cell developmental trajectories, on available experimental data. We show that our robust nonlinear manifold is well suited for raw, unfiltered gene counts from high throughput sequencing technologies for visualization and exploration of cell states. 62 factor analysis (ZIFA) (Pierson and Yau), t-SNE (Van Der Maaten and Hinton) (perplex-63 ity set to default 30), and PCA (Hotelling, 1933) were tested as comparison methods. To 64 evaluate the robust adaptations of the tGPLVM model, we fit tGPLVM with only an SE 65 kernel or SE and Matérn 1/2 kernel as well as tGPLVM with normally distributed error. 66 Trajectory building with minimum spanning trees 67 tGPLVM was used to fit a two dimensional latent mappings for the Lonnberg (Lönnberg 68 et al., 2017) developmental data. The minimum spanning tree was found on Euclidean 69 distance matrix of the posterior means of the low dimensional embedding and compared 70 to sequencing time to verify correct ordering. The same analysis was performed with 71 ZIFA (Pierson and Yau), t-SNE (Van Der Maaten and Hinton) (perplexity set to default 72 30), and PCA (Hotelling, 1933). 73 Visualization of sparse, raw count matrices 74 tGPLVM was used to fit a three-dimensional mapping for the two 10x data sets, CD34+ 75 cells and mice brain cells. Pearson correlation between latent position posterior mean and 76 expression counts was used to identify genes associated with latent dimensions. ZIFA (Pier-77 son and Yau), t-SNE (Van Der Maaten and Hinton) (perplexity set to default 30), and 78 PCA (truncated SVD (Halko et al.)) were also fit to the CD34+ cell data to compare 79 computational time.