scRICA: an R package for multiple-sample single-cell RNA-Seq Integrative Comparative Analysis

5.4.1 Cellular compositions summary and visualization

1. These can be done via getClusterSummaryReplot(). In our demonstrated data named as ‘scRICA_demo.rds’, where ‘expCond1’ specify the patient donor IDs, and ‘expCond2’ specify fallopian tube anatomical sites (I/A/F). If you would like to update the cells annotation, just run on newAnnotation=T, and provide the corresponding annotation file with the option newAnnotationRscriptName=/full/path/to/file, an example of new annotation file named as anno_example.R can be found from ‘misc’ folder, where all cell types need to be updated can be included.

   An example code to summarize the cellular composition based on samples attribute specified in ‘expCond1’ with respect to different donors is shown below:

rds1 <- readRDS('~/Desktop/757_scRNA-seq/RDS/scRNAseq_FT_finalized_wAnno.rds')
library(scRICA)
## seperation on donor: expCond1
getClusterSummaryReplot(rds = rds1, resDir = 'scRICA_results', newAnnotation = F, expCondCheck = 'expCond1', expCondCheckFname = 'donors_summary')
## seperation on fallopian tissue types (I/A/F): expCond2 togher with new annotation
getClusterSummaryReplot(rds = rds1, resDir = 'scRICA_results', newAnnotation = T, newAnnotationRscriptName = 'anno_example.R', expCondCheck = 'expCond1', expCondCheckFname = 'donors_summary')

2. Depending on whether newAnnotation is True or False, all related results will be saved into scRICA_results/results_wOrgClusterAnnotation/expCond_donors_summary or scRICA_results/results_wOrgClusterAnnotation/expCond_donors_summary respectively. An example of data structure of result with newAnnotation=F can be seen as:

   -- scRICA_results
     |__results_wOrgClusterAnnotation
   |expCond_donors_summary
      |__cellNo_summary_orgClusterAnnotation_expCond_donors_summary.pdf
      |__cellNo_summary_orgClusterAnnotation_expCond_donors_summary.txt
      |__new_tSNE_plot_expCond_donors_summary
         |__tsne_plot_noLabel_integrate_orgAnnotation.pdf
         |__tsne_plot_wLabel_integrate_orgAnnotation.pdf
         |__tsne_plot_wLabel_orgAnnotation_expCond_donors_summary.pdf
      |__new_UMAP_plot_expCond_donors_summary
         |__UMAP_plot_noLabel_integrate_orgAnnotation.pdf
         |__UMAP_plot_wLabel_integrate_orgAnnotation.pdf
         |__UMAP_plot_wLabel_orgAnnotation_expCond_donors_summary.pdf

• where ‘cellNo_summary_orgClusterAnnotation_expCond_donors_summary.pdf’: is a cellular composition plot separated by the specified sample’s attribute by the option expCondCheck = 'expCond1'.
• ‘cellNo_summary_orgClusterAnnotation_expCond_donors_summary.txt’: provides a summary of corresponding cell number and percentage of cells normalized by the total number of cells from each attribute.
• cellular composition plot separated by patient donors (expCondCheck = 'expCond1') is:

• cellular composition plot separated by fallopian tube anatomical sites (expCondCheck = 'expCond2') with new cell annotation specified in newAnnotation = T, newAnnotationRscriptName = 'anno_example.R' is:

• UMAP plot separated into 3 different anatomical sites (I/A/F) together with new cell annotations is:

  

3. There are several options which can be used to alter the cellular composition percentage bar-plot, such as

   • stack =T: make stacked or un-stacked cellular composition bar-plot by turning on/off this option

   • perPlotVertical = F, either in the vertical or horizontal position for this cellular composition bar-plot

   • fpRemake = F, perRemake = T: only make cellular composition bar-plot by turn on perRemake without making corresponding UMAP/tSNE plots (perPlotVertical = F)

   • perPlotVertical = F, perPlotHeight = 3, perPlotWidth = 7, fpRemake = F, stack =T,

   • expCondNameOrder = c('I', 'A', 'F'): specify the order of attributes specified from expCondCheck = 'expCond2', where ‘expCond2’ has 3 levels, which can be re-ordered as ‘I->A->F’.

   • An example of un-stacked re-ordered cellular composition plot shown below can be made by:

   getClusterSummaryReplot(rds = rds1, resDir = 'misc/visualization_results', 
                      perRemake = T, fpRemake = F, stack =F,  
                      perPlotVertical = T, perPlotHeight = 3, perPlotWidth = 7, 
                      expCondNameOrder = c('I', 'A', 'F'),
                      newAnnotation = F, 
                      expCondCheck = 'expCond2', expCondCheckFname = 'tissue_sep1')

4. Other options can be found from the help document via ?getClusterSummaryReplot

getGoiHeatmap(heatmap.view = 'expCond.idents', resDir = 'misc/visualization_results', rds = rds1,
            newAnnotation = F, geneNames = c('CD44', 'ITGA6', 'EPCAM', 'KRT5', 'KRT7', 'TUBB4A', 'TUBB4B', 'PAX8'),
            expCondCheck = 'expCond2',
            expCondCheckFname = 'heatmap_tissue',
            plotFnamePrefix = 'geneset1_expCond_idents_merge_SECEonly_IAonly',
            expCondReorderLevels = c('I CE', 'A CE', 'I SE', 'A SE'),
            display = 'sample.merge', draw.lines = F,
            cellcluster = c('CE', 'SE'),
            expCond = c('I', 'A'),
            fontsize.top = 2, fontsize.y = 12, barHeight.top = 0.05,
            plotWidth = 5, plotHeight = 3, fontangle.top = 90, debug = T)

5.5.1 Differential expression analysis between different samples attributes for cell types

• DE analysis on cells

  scRICA allow users to conduct DE analysis on any 2 group of cells with different sample’s attributes defined in any specific metadata table column (expCond*). For example, the demonstration RDs data has 2 columns, expCond1 and expCond2, defining samples attributes of donors and fallopian tissue types respectively. If we would like to check DEGs between tissue types I and F for a specific cell type with option cellcluster or all cell types without specifying cellcluster. Once you are running on all cell clusters, all cell types DE results are saved into the excel file spreadsheet, which can cause memory deficiency for excel file storage. commands to run these 2 analysis can be seen as below:

  de.res1 <- getClusterExpCondDe(resDir = 'misc/', rds = rds1,
                             newAnnotation = F, 
                             cellcluster = 'MP', deMethod = 'wilcox',
                             expCondCheck = 'expCond2', expCondCheckFname = 'IFcomp_wilcox', 
                             compGroup = 'I/F')

  de.res2 <- getClusterExpCondDe(resDir = 'misc/', rds = rds1,
                             newAnnotation = F, deMethod = 'MAST',
                             outputExcel = as.logical(F),
                             expCondCheck = 'expCond2', expCondCheckFname = 'IFcomp_allCells_MAST',
                             compGroup = 'I/F')

  All analysis results will be saved into the folder specified by resDir= 'misc' with newAnnotation = F into misc/results_wOrgClusterAnnotation_DEGs or misc/results_wNewAnnotation_DEGs when newAnnotation = T. For each different analysis, results will be saved into specific folder named via option expCondCheckFname, as shown below, above 2 analysis are saved into folder misc/results_wNewAnnotation_DEGs/IFcomp_wilcox and misc/results_wNewAnnotation_DEGs/IFcomp_allCells_MAST respectively, including 1) full list of DEGs (differentially expressed genes) results; 2) significantly DEGs at the specified pAdjValCutoff level, by default it is 0.05; 3) significantly up- and down-regulated DEGs at the specified pAdjValCutoff level; 4) summary of number of DEGs; 5) top number of DEGs list. When outputExcel = F, full list of DE results and significantly expressed DEGs are saved into the spreadsheet of the same excel file, otherwise, they are saved into the sperate txt file with respect to each cell type.

  -- results_wOrgClusterAnnotation_DEGs
     |__IFcomp_allCells_MAST
    |__expCondCompDeMarkers_I-F_adjSig_allClusters_cluster*.txt
    |__expCondCompDeMarkers_I-F_adjSig_down_allClusters_cluster*.txt
    |__expCondCompDeMarkers_I-F_adjSig_up_allClusters_cluster*.txt
    |__expCondCompDeMarkers_I-F_full_allClusters_cluster*.txt
    |__expCondCompDeMarkers_I-F_NoDEmarkers_summary.xlsx
    |__expCondCompDeMarkers_I-F_top10_downDe.xlsx
    |__expCondCompDeMarkers_I-F_top10_upDe.xlsx
     |__IFcomp_wilcox
    |__expCondCompDeMarkers_I-F_adjSig_1SelClusters.xlsx
    |__expCondCompDeMarkers_I-F_adjSig_down_1SelClusters.xlsx
    |__expCondCompDeMarkers_I-F_adjSig_up_1SelClusters.xlsx
    |__expCondCompDeMarkers_I-F_full_1SelClusters.xlsx
    |__expCondCompDeMarkers_I-F_NoDEmarkers_summary.xlsx
    |__expCondCompDeMarkers_I-F_top10_downDe.xlsx
    |__expCondCompDeMarkers_I-F_top10_upDe.xlsx

  Above returned DE analysis results are listed item, where you can access all DE analysis results with below command with de.res1 as an example:

  head(de.res1$clusterDeMarkers$MP)
  de.res1$clusterDeResSummar
  head(de.res1$clusterTopDeMarkers$MP)

• DE analysis on sample’s pseudo-bulk counts

  scRICA also provides DE analysis based on pseudo-bulk counts aggregated on samples defined in the initial metadata table, this information is inherited into the integrated RDS metadata table column expCond (rds1$expCond). For example, in our demonstration RDS file, there are totally 18 samples from different donor fallopian tube tissues, then we can use option deMethod = 'DESeq2.bulk' to conduct pseudo-bulk DE analysis, where counts from each sample are aggregrated and normalized based on the total number of cells. Additionally, users can also use an addtional metadata table to specify certain confounding covariates for model to correct those potential confounding effects for DE analysis with DESeq2. In this confounding factor defining metadata table, 1) column bulksamp should include all samples listed in rds1$expCond or initial metadata table for integration analysis; 2) other confounding factors to be corrected should be included in the columns of this table, such as we make an example confounding factor metadata table for our RDS demonstration data as below, here we would like to address potential confounding factors effects from gender and bmi.

  meta1 <- data.frame('bulksamp' = names(table(rds1$expCond)),
                  'gender' = c('F', 'F', 'F', 'M', 'M', 'M', 'M', 'M', 'F', 'F', 'M', 'M', 'M', 'M', 'M', 'M', 'F', 'F' ),
                  'bmi' = c(23.5, 23.5, 23.5, 18, 18, 28.4, 28.4, 28.4, 32, 32, 20, 20, 20, 30, 30, 30, 21, 21))

  With this confounding factor metadata table, we can use option deseq2bulk.metaCovariateInput to remove all potential confounding factor effects from pseudo-bulk DE analysis performed via DESeq2. Users also can use option norm.method to adjust normalization methods implemented across samples, including TMM, upper-quantile (UQ), and default DESeq2 normalization methods. An example of running command can be seen as below:

  de.res3 <- getClusterExpCondDe(resDir = 'misc/', rds = rds1,
                             newAnnotation = F, cellcluster = 'MP',
                             deMethod = 'DESeq2.bulk',
                             deseq2bulk.metaCovariateInput = meta1,
                             expCondCheck = 'expCond2', expCondCheckFname = 'IFcomp_DEseq2bulk',
                             compGroup = 'I/F')

  In addition to above analysis results, scRICA also return original and cell number normalized count table, which users can use below command to access them:

  head(de.res3$normCounts$MP)
  head(de.res3$orgCounts$MP)

5.5.2 Over expressed marker genes across comparison cell types with respect to different sample’s attributes levels.

scRICA allow users to check over expressed marker genes across defined cell types via function getExpCondClusterMarkers(). For example, users can check over-expressed genes among 5 selected immune cells, including T/NK1, T/NK2, T/NK3, MP, and B/P with respect to samples tissue types defined in expCond2 with below command:

oe.res1 <- getExpCondClusterMarkers(resDir = 'misc/', rds = rds1,
                                    deMethod = 'MAST',
                                    cellcluster = c('T/NK1', 'T/NK2', 'T/NK3', 'MP', 'B/P'),
                                    newAnnotation = F,
                                    expCondCheck = 'expCond2',
                                    expCondCheckFname = 'tissue_IR_overExpressedGenes')

Where, over-expressed marker genes are identified from each cell type in comparison with all other cell types with respect to each fallopian tube tissue type (I/A/F). Similarly, as other DE results, all analysis results are also saved into the folder based on defined options resDir, newAnnotation, and expCondCheckFname. In our example, results are saved into the foler named as misc/clusterMarkerGenes_results_wOrgClusterAnnotation with below data structure, where 1) full list of marker genes analysis results with respect to each fallopian tube tissue type are saved into txt file named as *_resFull.txt; 2) significant over-expressed genes list are also saved into different txt file named as *_UpSig.txt; 3) the number of over-expressed marker genes are also summarized into *_UpSig_NoSummary.txt; 4) top expressed marker genes heatmap is also made with respect to different levels of samples attributes as *_UpSig_heatmap.pdf.

-- clusterMarkerGenes_results_wOrgClusterAnnotation
   |__tissue_IR_overExpressedGenes
      |__expCond_tissue_IR_overExpressedGenes_A_allCluster_pos_markers_resFull.txt
      |__expCond_tissue_IR_overExpressedGenes_A_allCluster_pos_markers_UpSig_heatmap.pdf
      |__expCond_tissue_IR_overExpressedGenes_A_allCluster_pos_markers_UpSig_NoSummary.txt
      |__expCond_tissue_IR_overExpressedGenes_A_allCluster_pos_markers_UpSig.txt
      |__expCond_tissue_IR_overExpressedGenes_F_allCluster_pos_markers_resFull.txt
      |__expCond_tissue_IR_overExpressedGenes_F_allCluster_pos_markers_UpSig_heatmap.pdf
      |__expCond_tissue_IR_overExpressedGenes_F_allCluster_pos_markers_UpSig_NoSummary.txt
      |__expCond_tissue_IR_overExpressedGenes_F_allCluster_pos_markers_UpSig.txt
      |__expCond_tissue_IR_overExpressedGenes_I_allCluster_pos_markers_resFull.txt
      |__expCond_tissue_IR_overExpressedGenes_I_allCluster_pos_markers_UpSig_heatmap.pdf
      |__expCond_tissue_IR_overExpressedGenes_I_allCluster_pos_markers_UpSig_NoSummary.txt
      |__expCond_tissue_IR_overExpressedGenes_I_allCluster_pos_markers_UpSig.txt

Similarly, some analysis results can also be accessed from returned analysis object, including the full list of over-expressed genes and significantly over-express gene list, such as oe.res1 from our above analysis example.

head(oe.res1$fullRes$A)
head(oe.res1$sigUp$A)
table(oe.res1$sigUp$A$cluster)

5.5.3 HIPPO cells sub-clustering analysis

scRICA allow users to conduct sub-clustering HIPPO analysis for any specific cell types (cellcluster) with respect to cell from all samples or any specific sample’s attribute defined via expCondCheck and expCond together. This analyis can be done inside R with R function getHippoRes() or by command line exections Rscript run_scRICA_hippo.R [Options]....

• R function exectution: ?getHippoRes()

  For example, you can sub-cluster 3 types of T/NK cells (T/NK1, T/NK2 and T/NK3) into 5 sub-clusters via noClusters=5 based on zero-inflation clustering method for all cells regardless of theirs samples attributes, the command is as below:

  getHippoRes(resDir = 'misc/', rds = rds1, newAnnotation = F,
          hippoResNamePrefix = 'Allcells_TNK',
          cellcluster = c('T/NK1', 'T/NK2', 'T/NK3'), noClusters = 4)

  Users also can conduct HIPPO sub-clustering analysis on cells defined from a specific sample’s attribute level, such as you can conduct HIPPO sub-clustering analysis on 2 groups of T/NK cells from fallopian tissue F only via below command:

  getHippoRes(resDir = 'misc/', rds = rds1, newAnnotation = F,
          expCondCheck = 'expCond2', expCond = 'F',
          hippoResNamePrefix = 'Fonly_TNK23',
          cellcluster = c('T/NK2', 'T/NK3'), noClusters = 5)

  If too many cells need to be performed, you also can turn on sparseMatrix = T to avoid memory deficiency. The HIPPO analysis results are saved into the defined folder named via option hippoResNamePrefix inside resDir direcotry sub-folder hippo_results, an example of above analysis results data structure can be seen as below:

  -- hippo_results
     |__Allcells_TNK
    |__Allcells_TNK_k4_clustersTree.pdf
    |__diagnosticPlot_cellCluster_Allcells_TNK_k4.pdf
    |__lightHIPPOres_Allcells_TNK_k4.Rdata
    |__Per_cellCluster_Allcells_TNK_k4.xlsx
    |__topFeatures_cellCluster_Allcells_TNK_k4_dotplot.pdf
     |__Fonly_TNK23
    |__...

  where includes 1) hierarchy clustering tree plot; 2) diagnostic clustering plot; 3) all analysis results saved in Rdata; 4) number of cells in each HIPPO identified sub-clusters in excel file; 5) top expressed genes with respect to each HIPPO identified sub-clusters.

• Command line execution: Rscript run_scRICA_hippo.R --help

  Users can also conduct above analysis with command line executions as below:

  Rscript run_scRICA_hippo.R --outputDir misc/ --hippoResNamePrefix Allcells_TNK_commandLine --expCondCheck 'all' --rdsFname ~/Desktop/scRICA_demo.rds --cellclusters 'T/NK1, T/NK2, T/NK3' --noCluster 4 > $PWD/run_scRICA_hippo_tnk_allcells_k4.log 2>&1

  Rscript run_scRICA_hippo.R --outputDir misc/ --hippoResNamePrefix Fonly_TNK23_commandLine --expCondCheck 'expCond2' --expCond 'F' --rdsFname ~/Desktop/scRICA_demo.rds --cellclusters 'T/NK2, T/NK3' --noCluster 5 > $PWD/run_scRICA_hippo_tnk23_Fonly_k5.log 2>&1

• Update RDS file with HIPPO results

  As shown in above examples, 3 types of T/NK cells were conducted with sub-clustering analysis with HIPPO for a k=4 sub-clusters identification. Below command can be used to update the cells clustering information saved in Idents() based on HIPPO results saved in Rdata format. If list names are provided in the option hippoResList, updated cells idents names are based on them. Otherwise, users also can use option hippoResNames for updated HIIPO sub-clusters names. More detailed information regarding to these options, please refer ?updateHippoIdents.

  rds.hippo1 <- updateHippoIdents(resDir = 'misc/', rds = rds1, resSaveFname = 'hippo_tnk.Rds',
                           hippoResList = list('TNK_hippo' = 'misc/hippo_results/Allcells_TNK_k4/lightHIPPOres_Allcells_TNK_k4.Rdata'),
                           hippoResK = c(4),
                           newAnnotation = F )
  table(Idents(rds.hippo1))

5.5.4 pseudo-time trajectory analysis with slingshot

Similarly, scRICA also allow users to conduct pseudo-time trajectory analysis with slingshot for any specific cell types (cellcluster) with respect to cells from all samples or any specific sample’s attribute defined via expCondCheck and expCond together. This analyis can be done inside R with R function getExpCondClusterPseudotime() or by command line exections Rscript run_scRICA_pseudoTime.R [Options]....

• R function exectution: ?getExpCondClusterPseudotime(). For example, you can run slingshot pseudo-time trajectory analysis on SE and CE cell types with respect to 3 different fallopian tube anatomical sites seperations with below command:

  psdt.res1 <- getExpCondClusterPseudotime(resDir = 'misc/visualization_results', rds = rds1,
                                     expCondCheck = 'expCond2', 
                                     expCondCheckFname = 'tissueSep_SECE',
                                     cellcluster = c('SE1', 'CE'),
                                     newAnnotation = F)

• The above analysis results are saved inside folder named via options resDir and expCondCheckFname. For example, the above analysis results are saved into the folder named as ./misc/visualization_results/results_wOrgClusterAnnotation_pseudoTime/tissueSep_SECE with below data structure:

  -- results_wOrgClusterAnnotation_pseudoTime
   |__tissueSep_SECE
    |__expCondLevel_A_cluster_SE_fnPseudoTimeRes.Rdata
    |__expCondLevel_A_cluster_SE_ptGMM_GAMheatmap.pdf
    |__expCondLevel_A_cluster_SE_ptGMM_GAMheatmapGenenames.txt
    |__expCondLevel_A_cluster_SE_ptGMM_lineage.pdf
    |__expCondLevel_A_cluster_SE_ptGMM_plots.pdf
    |__...

  where includes results with respect to each samples attributes condition and cell types, such as for fallpian site ampulla, SE cell, 1) slingshot linage plots (expCondLevel_A_cluster_SE_ptGMM_lineage.pdf); 2) PCA and diffusion map seperation plots (expCondLevel_A_cluster_SE_ptGMM_plots.pdf); 3) heatmap of identified genes seperating lingeages (expCondLevel_A_cluster_SE_ptGMM_GAMheatmap.pdf); 4) corresponding gene list presented in heatmap (expCondLevel_A_cluster_SE_ptGMM_GAMheatmapGenenames.txt).

• Command line execution: Rscript run_scRICA_pseudoTime.R --help

  For above analysis, it can also be done with command line option shown as below:

  Rscript run_scRICA_pseudoTime.R --outputDir misc/ --hippoResNamePrefix tissueSep_SECE --expCondCheck 'expCond2' --rdsFname ~/Desktop/scRICA_demo.rds --cellclusters 'SE, CE' > $PWD/run_scRICA_pseudoTime_CESEonly.log 2>&1