Mpath-package: Mpath: an analysis algorithm that maps multi-branching...

In Zouter/MPath_gh: Mpath: a single-cell RNAseq data analysis algorithm that maps multi-branching single-cell trajectories revealing cell progression during development

Description

This package provides a new algorithm that maps multi-branching cell developmental pathways and aligns individual cells along the continuum of developmental trajectories. Mpath computationally reconstructs cell developmental pathways as a multi-destination journey on a map of connected landmarks wherein individual cells are placed in order along the paths connecting the landmarks. To achieve that, it first identifies clusters of cells and designates landmark clusters each defines a discrete cellular state. Subsequently it identifies and counts cells that are potentially transitioning from one landmark state to the next based on transcriptional similarities. It then uses the cell counts to infer putative transitions between landmark states giving rise to a state transition network. After that, Mpath sorts individual cells according to their various stages during transition to resemble the landmark-to-landmark continuum. Lastly, Mpath detects genes that were differentially expressed along the single-cell trajectories and identifies candidate regulatory markers.

Examples

## Not run: 
#### Install and load Mpath package

install.packages("Mpath_1.0.tar.gz",repos = NULL, type="source")
library(Mpath)

#################################################
###### Analysis of mouse DC dataset GSE60783 ####
#################################################

path <- getwd()
setwd(paste(path,"/GSE60783",sep=""))

##### remove low detection rate genes

rpkmFile="TPM_GSE60783_noOutlier.txt";
rpkmQCFile="TPM_GSE60783_noOutlier_geneQC0.05anyGroup.txt";
sampleFile="sample_GSE60783_noOutlier.txt";
QC_gene(rpkmFile=rpkmFile,
        rpkmQCFile=rpkmQCFile,
        sampleFile=sampleFile,threshold=0.05,method="any")

### Mpath using spleenic CD4 vs CD8 DEGs 

rpkmFile = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG.txt";
sampleFile = "sample_GSE60783_noOutlier.txt";
find_optimal_cluster_number(rpkmFile = rpkmFile,
                            sampleFile = sampleFile,
                            min_cluster_num = 7, max_cluster_num = 15,
                            diversity_cut = 0.6, size_cut = 0.05)
 
### Landmark designation

rpkmFile = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG.txt";
baseName = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG";
sampleFile = "sample_GSE60783_noOutlier.txt";
landmark_cluster <- landmark_designation(rpkmFile = rpkmFile,
                                         baseName = baseName,
                                         sampleFile = sampleFile,
                                         method = "diversity_size",
                                         numcluster = 11, diversity_cut=0.6,
                                         size_cut=0.05)

### Plot hierachical clustering 

dataFile = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG.txt";
SC_hc_colorCode(dataFile = dataFile,
                cuttree_k = 11,
                sampleFile= "sample_GSE60783_noOutlier.txt",
                width = 22, height = 10, iflog2 = TRUE,
                colorPalette = c("red","green","blue"))

### Construct weighted neighborhood network

exprs = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG.txt";
baseName = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG";
neighbor_network <- build_network(exprs = exprs,
                                  landmark_cluster = landmark_cluster,
                                  baseName = baseName)

### TrimNet: trim edges of lower weights

baseName = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG";
trimmed_net <- trim_net(neighbor_network,textSize=30,
                        baseName = baseName,
                        method = "mst")

### plot trimmed net and color-code the nodes by gene expression

rpkmFile="TPM_GSE60783_noOutlier.txt";
lmFile="TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG_landmark_cluster.txt";
color_code_node_2(networkFile=trimmed_net,
                  rpkmFile=rpkmFile,
                  lmFile=lmFile,
                  geneName=c("Irf8","Id2","Batf3"),
                  baseName="cDC1_marker",
                  seed=NULL)

### Re-order the cells on the path connecting landmark 
### "CDP_2","CDP_1","PreDC_9","PreDC_3"

exprs = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG.txt";
ccFile = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup_CD4vsCD8DEG_landmark_cluster.txt";
order <- nbor_order(exprs = exprs,
                    ccFile = ccFile, 
                    lm_order = c("CDP_2","CDP_1","preDC_9","preDC_3"),
                    if_bb_only=FALSE,
                    method=1)

### identify genes that changed along the cell re-ordering

deg <- vgam_deg(exprs = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup.txt",
                order = order,
                lm_order = c("CDP_2","CDP_1","preDC_9","preDC_3"),
                min_expr=1,
                p_threshold=0.05)

### plot heatmap of genes that changed along the cell re-ordering

heatmap_nbor(exprs = "TPM_GSE60783_noOutlier_geneQC0.05anyGroup.txt",                        
             cell_order = "CDP_2_CDP_1_preDC_9_preDC_3_order.txt", 
             plot_genes = "CDP_2_CDP_1_preDC_9_preDC_3_vgam_deg0.05.txt",
             cell_annotation = "sample_GSE60783_noOutlier.txt", 
             num_gene_cluster = 6,
             hm_height = 15, hm_width = 10,
             baseName = "CDP_2_CDP_1_preDC_9_preDC_3_order_vgam_deg0.05")
             
###################################################
### Analysis of human myoblast dataset GSE52529 ###
###################################################

setwd(paste(path,"/GSE52529",sep=""))

### remove low detection rate genes

QC_gene(rpkmFile = "GSE52529_fpkm_matrix_nooutliers.txt",
        rpkmQCFile = "GSE52529_fpkm_matrix_nooutliers_geneQC0.05anyGroup.txt",
        sampleFile = "sample_nooutlier.txt",threshold=0.05,method="any")

rpkmFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt";
sampleFile = "sample_nooutlier.txt";
find_optimal_cluster_number(rpkmFile = rpkmFile,
                            sampleFile = sampleFile,
                            min_cluster_num = 7, max_cluster_num = 18,
                            diversity_cut = 0.9, size_cut = 0.05)

rpkmFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt";
baseName = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG";
landmark_cluster <- landmark_designation(rpkmFile = rpkmFile,
                                         baseName = baseName,
                                         sampleFile = "sample_nooutlier.txt",
                                         method = "diversity_size",
                                         numcluster = 14, diversity_cut=0.9, 
                                         size_cut=0.05)

SC_hc_colorCode(dataFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",
                cuttree_k = 14,
                sampleFile= "sample_nooutlier.txt",
                width = 22, height = 10, iflog2 = TRUE)

### Construct weighted neighborhood network

exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt";
neighbor_network <- build_network(exprs = exprs,
                                  landmark_cluster = landmark_cluster)

### TrimNet: trim edges of lower weights

trimmed_net <- trim_net(neighbor_network,textSize=30,
                       baseName = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG",
                       method = "mst")

### Color code the landmarks with marker expression

networkFile="GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_state_transition_mst.txt";
lmFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_landmark_cluster.txt";
rpkmFile = "GSE52529_fpkm_matrix_nooutliers_geneSymbol.txt";
geneName=c("SPHK1","PBX1","XBP1","ZIC1","MZF1","CUX1","ARID5B","POU2F1","CDK1","MYOG");
color_code_node_2(networkFile = networkFile,
                  rpkmFile = rpkmFile,
                  lmFile = lmFile,
                  geneName = geneName,
                  baseName = "Marker",
                  seed=3)


### path 1: "T0_2","T0_1","T24_8","T48_10","T72_13"

ccFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_landmark_cluster.txt";
order <- nbor_order(exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",
                    ccFile = ccFile, 
                    lm_order = c("T0_2","T0_1","T24_8","T48_10","T72_13"),
                    if_bb_only=TRUE,
                    method=1)

deg <- vgam_deg(exprs = "GSE52529_fpkm_matrix_nooutliers_geneQC0.05anyGroup.txt",
               order = order,
               lm_order = c("T0_2","T0_1","T24_8","T48_10","T72_13"),
               min_expr=1,
               p_threshold=0.05)

heatmap_nbor(exprs = "GSE52529_fpkm_matrix_nooutliers_geneQC0.05anyGroup.txt",                        
            cell_order = order,
            plot_genes = row.names(deg),
            cell_annotation = "sample_nooutlier.txt", 
            num_gene_cluster = 6,
            hm_height = 15, hm_width = 10,
            baseName = "GSE52529_path1_order_vgam_deg0.05")

### path 2: "T0_2","T0_1","T24_7","T48_4","T72_14"
ccFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_landmark_cluster.txt";
order <- nbor_order(exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",
                    ccFile = ccFile, 
                    lm_order = c("T0_2","T0_1","T24_7","T48_4","T72_14"),
                    if_bb_only = TRUE,
                    method=1)

deg <- vgam_deg(exprs = "GSE52529_fpkm_matrix_nooutliers_geneQC0.05anyGroup.txt",
               order = order,
               lm_order = c("T0_2","T0_1","T24_7","T48_4","T72_14"),
               min_expr=1,
               p_threshold=0.05)

heatmap_nbor(exprs = "GSE52529_fpkm_matrix_nooutliers_geneQC0.05anyGroup.txt",                        
            cell_order = order,
            plot_genes = row.names(deg),
            cell_annotation = "sample_nooutlier.txt", 
            num_gene_cluster = 6,
            hm_height = 15, hm_width = 10,
            baseName = "GSE52529_path2_order_vgam_deg0.05")

### Generate Figure 5

ccFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_landmark_cluster.txt";
order1 <- nbor_order(exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",
                     ccFile = ccFile, 
                     lm_order = c("T0_2","T0_1","T24_8","T48_10","T72_13"),
                     if_bb_only=TRUE,
                     method=1)

ccFile = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG_landmark_cluster.txt";
order2 <- nbor_order(exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",
                     ccFile = ccFile, 
                     lm_order = c("T0_2","T0_1","T24_7","T48_4","T72_14"),
                     if_bb_only=TRUE,
                     method=1)

deg1 <- read.table("T0_2_T0_1_T24_8_T48_10_T72_13_vgam_deg0.05.txt",sep="\t",header=T)
deg2 <- read.table("T0_2_T0_1_T24_7_T48_4_T72_14_vgam_deg0.05.txt",sep="\t",header=T)
deg <- unique(c(as.character(deg1[,1]),as.character(deg2[,1])))

heatmap_nbor(exprs = "GSE52529_fpkm_matrix_nooutliers_ANOVA_p0.05_DEG.txt",                        
            cell_order = c(order1,order2), 
            plot_genes = deg,
            cell_annotation = "sample_nooutlier.txt", 
            num_gene_cluster = 7,
            hm_height = 15, hm_width = 10,
            baseName = "Path12_method1orderedbackbone_progression_heatmap", 
            n_linechart = list(order1,order2))

## End(Not run)
             

Zouter/MPath_gh documentation built on May 14, 2019, 8:58 a.m.