---Integrative scRNA-Seq Interpretation System for Functional Gene Module analysis

Single-cell RNA-Seq data is useful in discovering cell heterogeneity and signature genes in specific cell populations in cancer and other complex diseases. Specifically, the investigation of functional gene modules (FGM) can help to understand gene interactive networks and complex biological processes. QUBIC2 is recognized as one of the most efficient and effective tools for FGM identification from scRNA-Seq data. However, its availability is limited to a C implementation, and its applicative power is affected by only a few downstream analyses functionalities. We developed an R package named IRIS-FGM (integrative scRNA-Seq interpretation system for functional gene module analysis) to support the investigation of FGMs and cell clustering using scRNA-Seq data. Empowered by QUBIC2, IRIS-FGM can identify co-expressed and co-regulated FGMs, predict types/clusters, identify differentially expressed genes, and perform functional enrichment analysis. It is noteworthy that IRIS-FGM also applies Seurat objects that can be easily used in the Seurat vignettes.

1. System environment:

We recommend users to install IRIS-FGM on large memory (32GB) based UNIX/Linux operation system if the user aims at analyzing bicluster-based co-expression analysis. If the user seeks to analyze data by quick mode, we recommend installing IRIS-FGM on a small memory (8GB) based Windows or Linux operation system; IRIS-FGM does not support MAC.

2. R environment: R (equal or greater than 3.5)

1. Install required packages from CRAN:

install.packages(c('BiocManager','devtools', 'AdaptGauss', "pheatmap", 'mixtools','MCL', 'anocva', 'qgraph','Rtools','ggpubr',"ggraph","Seurat"))

2. Install required packages from Bioconductor:

BiocManager::install(c('org.Mm.eg.db','multtest','org.Hs.eg.db','clusterProfiler','DEsingle', 'DrImpute','scater', 'scran'))

3. Install IRIS-FGM from github:

devtools::install_github("BMEngineeR/IRIS-FGM", force = T)

Now the MetaQUBIC is successfully installed and ready for use.

1. Download example data from Yan's RPKM 90 cell embroyonic single cell RNA-Seq data, meta file (cell type information) and paper. The alternative way to download is to use the step 2 method in this section.

2. Set the working directory where you put data and import the IRIS-FGM library.

setwd("./my_working_dir/")
library(IRISFGM)
# this step is for download Yan's expression data and meta file.
download.file(url = "https://bmbl.bmi.osumc.edu/downloadFiles/Yan_expression.txt",destfile = "./my_working_dir/Yan_expression.txt")
download.file(url = "https://bmbl.bmi.osumc.edu/downloadFiles/Yan_cell_label.txt",destfile = "./my_working_dir/Yan_cell_label.txt")

InputMatrix <- read.table("./my_working_dir/Yan_expression.txt",header = T, row.names = 1,check.names = F)
Meta_info <- read.table("./my_working_dir/Yan_cell_label.txt",header = T,row.names = 1)

IRIS-FGM also provides the function to read in 10X scRNA-Seq data format.

Read HDF5 fileReadFrom10X_h5() or read the folder which contain three files (barcode, sparse matrix, genes)ReadFrom10X_folder()

1. Create IRIS-FGM object:

object <- CreateIRISCEMObject(InputMatrix)

2. Adding meta information for your cell:

This step can add the customized cell label by the user. The format of a file passing to meta.info is a data frame of which row name should be cell ID, and column name should be cell type.

object <- AddMeta(object, meta.info = Meta_info)

3. Filtering out low quality data:

Use PlotMeta and SubsetData together to filter out low-quality cells.

PlotMeta(object)

object <- SubsetData(object , nFeature.upper=15000,nFeature.lower=8000,
                         Counts.upper=700000,Counts.lower=400000)

After cells get filtered, we can check the distribution again.

PlotMeta(object)

4. Normalization:

Users can choose to perform normalization based on their needs. The normalization method has two options, one is the simplest CPM normalization (default normalization = 'LibrarySizeNormalization'). The other is from package scran and can be opened by using parameter normalization = 'scran'. Compared to the CPM normalization method, scran will be more accurate but takes more time.

object <- ProcessData(object, normalization = "LibrarySizeNormalization")

5. LTMG modeling:

Here, we will use Left-truncated Mixture Gaussian distribution to model the regulatory signal of each gene. Parameter, 'Gene_use', decides number of top highly variant gene for LTMG modeling, and here we use top 2000 highly variant genes.

object <- RunLTMG(object, Gene_use = "2000", seed = 123)

IRIS-FGM can provide biclustering function, which is based on our in-house novel algorithm, [QUBIC2] (https://github.com/maqin2001/qubic2) to predict functional gene module.

Discretization & biclustering
If you have more cells to analyze functional gene module you can use LTMG or Quantile based discretization (QUBIC). In this step, IRIS-FGM will generate files in the local working directory and do not remove them before finishing your analysis (files include tmp_expression.txt, tmp_expression.txt.chars, tmp_expression.txt.chars.blocks, tmp_expression.txt.chars.chars, and tmp_expression.txt.rules).

1. LTMG based discretization:

You need to run LTMG modeling then binarizing them. It might take a lot of time if you have a large size of cells.

object <- RunLTMG(object, Gene_use = "all", seed = 123)
object <- CalBinaryMultiSignal(object)
object <- RunBicluster(object, DiscretizationModel = "LTMG",OpenDual = FALSE,
                          NumBlockOutput = 100, BlockOverlap = 0.7, BlockCellMin = 15)

2. Quantile based discretization:

It might require users to adjust parameter q for deciding cutoff of binarizing.

object <- RunDiscretization(object, q = 0.06)
object <- RunBicluster(object, DiscretizationModel = "Quantile",OpenDual = FALSE, Extension = 0.90,
                          NumBlockOutput = 1000, BlockOverlap = 0.7, BlockCellMin = 15)

Analyzing functional gene module This section is based on the quantile based discretization and biclustering results.

1. Visualize the gene module-to-module relationship globally:

The figure via the number of overlap genes (controlled edge.by = "gene") shows that bicluster 20, 23, 24, and 25 have similar gene components. Bicluster 21, 22, 25, and until bicluster 35 have similar gene components. In this figure, we can see gene modules in bicluster 20, 23, 24, and 25 may have a similar component and similar functionality, whereas the gene components from this gene module may differ from the other gene modules from the other biclusters (i.e., bicluster 21, 22, 25, and until bicluster 35)

PlotNetwork(object,N.bicluster =c(20:30), edge.by = "gene")

2. Visualize bicluster 20 & bicluster 35 on heatmap.

From step 1 in "Analysis functional gene module," we postulate gene components from bicluster 20, 23, 24, and 25 may differ from gene components from bicluster 21, 22, 25, and until bicluster 35. Therefore, in this section, we will focus on how the difference is. Therefore, we use bicluster 20 and bicluster 35 to generate heatmap and show such a difference.

PlotHeatmap(object,N.bicluster = c(20,35),show.annotation = F)

3. Visualize local co-expression gene module network:

Since we already know the bicluster 20 and bicluster 35 showing the difference at the global level. We then will focus on a local gene module and investigate the co-expression gene network with in the module. Yellow nodes represent the gene module network from bicluster #20. The nodes' size indicates the degree of presence (the number of connected edges for one node). The thickness of edges suggests the value of the correlation coefficient.

• From this figure (bicluster 20, click the figure for zooming in), we can tell the EIFAD gene show a negative correlation (red color edge) to GOSR1 & BBS5, and also show a positive correlation (grey edge) to ZNF394 & POTEM in a group of cells from this bicluster.

PlotModuleNetwork(object, N.bicluster = 20, cutoff=0.6, Node.color = "#E8E504")

• From this figure (bicluster 35, click the figure for zooming in), we can tell the ROBO1 gene shows a negative correlation (red color edge) to NLRP4 & BPGM. GNPDA2 gene shows a positive correlation to BPGM, KIT, and CCDC25 in a group of cells from this bicluster.

PlotModuleNetwork(object, N.bicluster = 35, cutoff=0.6, Node.color = "#E8E504")

4. Functional enrichment analysis.

ISIR-FGM provide a functional enrichment analysis for a selected gene module.

1. For gene module from bicluster 20:

object <- RunPathway(object ,module.number =20, selected.gene.cutoff = 0.05, species = "Human", database = "GO", genes.source = "Bicluster")

2. For gene module from bicluster 35:

object <- RunPathway(object ,module.number =35, selected.gene.cutoff = 0.05, species = "Human", database = "GO", genes.source = "Bicluster")

IRIS-FGM provides the function to use the output biclusters to predict cell type based on the MCL algorithm, which predicted results outperformed the other nine bicluster and cluster methods (Details can be found from figure B in IRIS-FGM paper). In this case, we identify eight cell clusters.

PlotModuleNetwork(object, N.bicluster = 35, cutoff=0.6, Node.color = "#E8E504")
unique(object@MetaInfo$MC_Label)
[1] 1 2 3 4 5 6 7 8

IRIS-FGM integrates the Seurat package for carrying out dimension reduction. Furthermore, IRIS-FGM also integrates the Seurat package cell clustering method as an alternative approach. Before running dimension reduction, you should perform LTMG modeling if you would like to LTMG signal matrix as input rather than the processed expression matrix. The advantage of using the LTMG signal matrix as input is that it can improve dimension reduction and cell clustering. (Detail can be found in this LTMG paper)

1. Use processed expression matrix as input to generate UMAP and cell clusters based on the Seurat package.

• runing program:

object <- IRISFGM::RunDimensionReduction(object,mat.source = "UMImatrix")
object <- IRISFGM::RunClassification(object,resolution = 0.8)

• Visualization by using different cell label: "Cluster" is the benchmark label from paper. "MC_Label" is the label based on MCL prediction. "Seurat0.6" is the label based on the Seurat cell clustering method.

# This function will ask you to choose an index (represent a cluster identity).
# Use benchmark cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
select index of cell condition: 3

# Use MCL predicted cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
select index of cell condition: 4

# Use Seurat predicted cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
select index of cell condition: 5

2. Use LTMG signal matrix as input to generate UMAP and cell clusters based on the Seurat package.

• runing program:

object <- IRISFGM::RunDimensionReduction(object,mat.source = "LTMG")
object <- IRISFGM::RunClassification(object,resolution = 0.8)

• Visualization by using different cell label: "Cluster" is the benchmark label from paper. "MC_Label" is the label based on MCL prediction. "Seurat0.6" is the label based on the Seurat cell clustering method.

# This function will ask you to choose an index (represent a cluster identity).
# Use benchmark cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
6 : Seurat0.8
select index of cell condition: 3

# Use MCL predicted cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
6 : Seurat0.8
select index of cell condition: 4

# Use Seurat predicted cell clusters.
PlotDimension(object)
select condition to present
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
6 : Seurat0.8
select index of cell condition: 6

Cell clusters have been identified by a different method. We can use these cell clusters to analyze DEGs and perform functional enrichment analysis.

1. Find DEG based on MCL prediction.

• Find DEGs.

This function will ask you to provide an index of a cluster label, select a target cluster (group 1), and select a comparing cluster (group 2). Here we select "ML_Label" cell clusters (index is 4), choose the 8th cluster (index is 3, as group 1), and compare to the rest of all clusters (index is 9, as group 2). The DEGs' result can be found in object@BiCluster@MarkerGene .

object <- IRISFGM::FindMarkers(object)
select condition to compare
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
6 : Seurat0.8

select index of cell condition: 4
select index (left) of first group to compare : 
1 : 1
2 : 2
3 : 3
4 : 4
5 : 5
6 : 6
7 : 7
8 : 8
input first group index : 3
select index (left) of second group to compare : 
1 : 1
2 : 2
3 : 3
4 : 4
5 : 5
6 : 6
7 : 7
9 : rest of all

select index of group 2: 9
# DEGs based on MCL predicted label was store in object@BiCluster@MarkerGene
object@BiCluster@MarkerGene[1:5,]

• Functional enrichment analysis.

The result can be found at object@BiCluster@PathwayFromMC.

object <- RunPathway(object, species = "Human", database = "GO", genes.source = "Bicluster")
object@BiCluster@PathwayFromMC[1:5,]                 

2. Find DEG based on Seurat prediction (LTMG as input matrix).

• Find DEGs.

This function will ask you to provide an index of a cluster label, select a target cluster (group 1), and select a comparing cluster (group 2). Here we select "Seurat0.8" (result based on inputing LTMG signal matrix) cell clusters (index is 6), choose the 1th cluster (index is 2, as group 1), and compare to the rest of all clusters (index is 4, as group 2). The DEGs' result can be found in object@LTMG@MarkerGene .

object <- IRISFGM::FindMarkers(object)
select condition to compare
1 : ncount_RNA
2 : nFeature
3 : Cluster
4 : MC_Label
5 : Seurat0.8
6 : Seurat0.8

select index of cell condition: 6
select index of cell condition: 14
select index (left) of first group to compare : 
1 : 0
2 : 1
3 : 2

input first group index : 2
select index (left) of second group to compare : 
1 : 0
3 : 2
4 : rest of all

select index of group 2: 4

# DEGs based on MCL predicted label was store in object@LTMG@MarkerGene
object@LTMG@MarkerGene[1:5,]

• Functional enrichment analysis.

The result can be found at object@LTMG@Pathway.

object <- RunPathway(object, species = "Human", database = "GO", genes.source = "CTS")
object@LTMG@Pathway[1:5,]