diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml
new file mode 100644
index 0000000..105ce2d
--- /dev/null
+++ b/.idea/inspectionProfiles/profiles_settings.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
\ No newline at end of file
diff --git a/.idea/misc.xml b/.idea/misc.xml
new file mode 100644
index 0000000..ac8300d
--- /dev/null
+++ b/.idea/misc.xml
@@ -0,0 +1,4 @@
+
+
+
+
\ No newline at end of file
diff --git a/.idea/modules.xml b/.idea/modules.xml
new file mode 100644
index 0000000..8a3509e
--- /dev/null
+++ b/.idea/modules.xml
@@ -0,0 +1,8 @@
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/.idea/scGCN.iml b/.idea/scGCN.iml
new file mode 100644
index 0000000..38f76ec
--- /dev/null
+++ b/.idea/scGCN.iml
@@ -0,0 +1,12 @@
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/.idea/vcs.xml b/.idea/vcs.xml
new file mode 100644
index 0000000..94a25f7
--- /dev/null
+++ b/.idea/vcs.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
\ No newline at end of file
diff --git a/.idea/workspace.xml b/.idea/workspace.xml
new file mode 100644
index 0000000..5466975
--- /dev/null
+++ b/.idea/workspace.xml
@@ -0,0 +1,63 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 1647938476993
+
+
+ 1647938476993
+
+
+
+
\ No newline at end of file
diff --git a/README.md b/README.md
index dc85a41..86f0459 100644
--- a/README.md
+++ b/README.md
@@ -1,75 +1,56 @@
-scGCN is a Graph Convolutional Networks Algorithm for Knowledge Transfer in Single Cell Omics
+scRNA_GCN流程是在scGCN流程的基础上进行的升级和完善，支持单细胞RNA基因表达矩阵作为输入，用以跨物种之间的细胞类型比较。
+## 概述
-This is a TensorFlow implementation of scGCN for leveraging and label transfer across differnt single cell datasets.
+基于卷积神经网络（Convolutional Neural Networks）的跨物种单细胞卷积神经网络细胞类型注释软件是一个linux平台的python软件，旨在帮助没有过多的编程知识的用户在下游分析单细胞转录组测序（scRNA-seq）及相关组学数据。软件提供了一种方便的算法来对不同物种的基因表达相似性进行分析，从而在无需知道跨物种之间的物种特异性marker的情况下，根据参考物种的细胞类型注释对待测物种的细胞类型进行注释。并且，在细胞类型注释完成后，还能简单地比较分析不同物种细胞类型的基因表达差异。软件操作难度低，利于没有编程基础的研究人员学习并使用。
-[](https://zenodo.org/badge/latestdoi/294531199)
-
-## Overview
-
-Single-cell omics represent the fastest-growing genomics data type in the literature and the public genomics repositories. Leveraging the growing repository of labeled datasets and transferring labels from existing datasets to newly generated datasets will empower the exploration of the single-cell omics. The current label transfer methods have limited performance, largely due to the intrinsic heterogeneity among cell populations and extrinsic differences between datasets. Here, we present a robust graph artificial intelligence model, single-cell Graph Convolutional Network (scGCN), to achieve effective knowledge transfer across disparate datasets. Benchmarked with other label transfer methods on different single cell omics datasets, scGCN has consistently demonstrated superior accuracy on leveraging cells from different tissues, platforms, and species, as well as cells profiled at different molecular layers. scGCN is implemented as an integrated workflow and provided here.
-
-## Requirements
+## 需要安装的包
* setuptools >= 40.6.3
* numpy >= 1.15.4
* tensorflow >= 1.15.0
* networkx >= 2.2
* scipy >= 1.1.0
-## Installation
+## 流程的安装
-Download scGCN:
+下载流程：
```
-git clone https://github.com/QSong-github/scGCN
+git clone https://github.com/Dee-chen/scGCN/tree/cdy
```
-Install requirements and scGCN:
+安装依赖:
```bash
python setup.py install
```
-The general installation time is less than 10 seconds, and have been tested on mac OS and linux system.
-## Run the demo
-load the example data using the data_preprocess.R script
-In the example data, we include the data from Mouse (reference) and Human (query) of GSE84133 dataset. The reference dataset contains 1,841 cells and the query dataset contains more cells (N=7,264) and 12,182 genes.
+## 运行示例
+
+在示例文档中，有两个物种的相关基因表达矩阵数据
```bash
cd scGCN
-Rscript data_preprocess.R # load example data
-python train.py # run scGCN
+Rscript data_preprocess.R # 加载和预处理示例数据
+python train.py # 运行主程序
```
-All output will be shown in the output_log.txt file. Performance will be shown at the bottom.
-We also provide the Seurat performance on this reference-qeury set (as in Figure 4), by run
-
```
Rscript Seurat_result.R
```
-## Input data
-
-When using your own data, you have to provide
-* the raw data matrix of reference data and cell labels
-* the raw data matrix of query data
-
-## Output
-
-The output files with scGCN predicted labels will be stored in the results folder.
-
-## Model options
-
-We also provide other GCN models includidng GAT (Veličković et al., ICLR 2018), HyperGCN (Chami et al., NIPS 2019) and GWNN (Xu et al., ICLR 2019) for optional use.
-
-## Detecting unknown cells
+## 输入数据
-For the query data that have cell types not appearing in reference data, we provide a screening step in our scGCN model using two statistical metrics, entropy score and enrichment score. If certain cells in query data have higher entropy and lower enrichment, these cells should be assigned as unknown cells. Specifically, choose check_unknown=TRUE in the function 'save_processed_data' to detect unknown cells.
+当使用自定义数据时，需要提供一个待测物种的基因表达矩阵，和一个参考物种的基因表达矩阵以及对应的细胞类型标签数据
+* 基因表达矩阵是数据框结构，行是细胞，列是基因
+* 细胞类型标签是单列数据框结构，以"type"作为单列数据的表头
-## Reproduction instructions
+## 输出
-The above scripts can reproduce the quantitative results in our manuscript based on our provided data.
+输出的数据信息会在results文件夹中，可视化结果则在当前运行的主程序目录下
-## Cite
+## 差异分析
+在原始scGCN的基础上，增加了初步的差异分析功能，用以可视化跨物种同一细胞类型的主要差异基因
-Please cite our paper and the related GCN papers if you use this code in your own work:
+## 引用
+原始scGCN的文章引用信息：
```
Song, Q., Su, J., & Zhang, W. (2020). scGCN: a Graph Convolutional Networks Algorithm for Knowledge Transfer in Single Cell Omics. bioRxiv.
```
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..7466fe9
--- /dev/null
+++ b/main.py
@@ -0,0 +1,98 @@
+import argparse
+import threading
+import time
+import os
+import configparser
+import logging
+import coloredlogs
+import sys
+import subprocess
+import multiprocessing
+import pandas as pd
+
+class NewConfigParser(configparser.ConfigParser):
+ def optionxform(self, optionstr):
+ return optionstr
+
+def err_exit():
+ sys.exit('\033[1;31;47m!!The program exited abnormally, please check the log file !!\033[0m')
+
+def main():
+ pwd= os.getcwd()
+ parser = argparse.ArgumentParser(description='scGCN pipline')
+ parser.add_argument('-i', type=str, required=True, metavar='input_list',help='The count matix and cell type.')
+ parser.add_argument('-o', type=str, metavar='outputdir', default=os.sep.join([pwd, 'output']),
+ help='The output dir.default:' + os.sep.join([pwd, 'output']))
+ parser.add_argument('--log', type=str, metavar='logfile',help='log file name,or log will print on stdout')
+ parser.add_argument('--debug', action='store_true', default=False, help='The log file will contain the output of each software itself, which is convenient for finding errors (-log is required)')
+
+ args = parser.parse_args()
+
+
+ if not args.log:
+ coloredlogs.install(
+ fmt='%(asctime)s: %(levelname)s\t%(message)s',
+ evel='info'.upper(), stream=sys.stdout
+ )
+ else:
+ if args.debug:
+ l='debug'
+ else:
+ l='info'
+ print('Logs will be written to {}'.format(args.log))
+ if os.path.exists(args.log):
+ os.remove(args.log)
+ logging.basicConfig(
+ filename=args.log,
+ filemode='a',
+ format='%(asctime)s: %(levelname)s\t%(message)s',
+ datefmt='%H:%M:%S',
+ level=l.upper()
+ )
+ home_dir=os.path.dirname(os.path.abspath(__file__))
+
+ logging.info("Start checking the input file..")
+ try:
+ input_list=pd.read_csv(args.i, sep='\t', header=0)
+ except FileNotFoundError:
+ logging.error('Checking input file :[ERR] -- The input file :%s does not exist' %(args.input_list))
+ err_exit()
+ if input_list.columns.tolist()!=["Species","count","type","dbclass"]:
+ logging.error('The input file :%s format error.' % (args.i))
+ err_exit()
+ for f in input_list['count']:
+ if not os.path.exists(f):
+ logging.error('The input count file :%s does not exist.'%f)
+ err_exit()
+ for f in input_list['type']:
+ if not os.path.exists(f):
+ logging.error('The cell type count file :%s does not exist.'%f)
+ err_exit()
+
+ outpath=os.path.abspath(args.o)
+ os.mkdir(outpath)
+
+ for ref in input_list[input_list.dbclass=="ref"].itertuples():
+ for query in input_list[input_list.dbclass=="query"].itertuples():
+
+ os.mkdir(os.sep.join([outpath,ref.Species+"_"+query.Species]))
+ logging.info('Start processing %s to %s type predictions.' % (ref.Species, query.Species))
+ logging.info('Data pre-processing...')
+ os.chdir(os.sep.join([outpath,ref.Species+"_"+query.Species]))
+ logging.info(['Rscript',os.sep.join([home_dir,'data_preprocess.R']),ref.count,ref.type,query.count,query.type])
+ stdout=subprocess.run(['Rscript',os.sep.join([home_dir,'data_preprocess.R']),ref.count,ref.type,query.count,query.type],stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+ logging.debug('{} stdout:\n'.format('Data pre-processing') +stdout.stdout.decode('utf-8'))
+ logging.debug('{} stderr:\n'.format('Data pre-processing') +stdout.stderr.decode('utf-8'))
+ logging.info("Data pre-processing finish.")
+
+ logging.info("Cell type predictions...")
+ stdout = subprocess.run(['python',os.sep.join([home_dir,'train.py'])],stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+ logging.debug('{} stdout:\n'.format('Data pre-processing') +stdout.stdout.decode('utf-8'))
+ logging.debug('{} stderr:\n'.format('Data pre-processing') +stdout.stderr.decode('utf-8'))
+ logging.info("Cell type predictions finish.")
+
+ logging.info("ALL finish.")
+
+if __name__ == '__main__':
+ main()
+
diff --git a/scGCN/.Rhistory b/scGCN/.Rhistory
new file mode 100644
index 0000000..e69de29
diff --git a/scGCN/__pycache__/data.cpython-38.pyc b/scGCN/__pycache__/data.cpython-38.pyc
new file mode 100644
index 0000000..60915b6
Binary files /dev/null and b/scGCN/__pycache__/data.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/graph.cpython-38.pyc b/scGCN/__pycache__/graph.cpython-38.pyc
new file mode 100644
index 0000000..3539653
Binary files /dev/null and b/scGCN/__pycache__/graph.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/layers.cpython-38.pyc b/scGCN/__pycache__/layers.cpython-38.pyc
new file mode 100644
index 0000000..235de3a
Binary files /dev/null and b/scGCN/__pycache__/layers.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/models.cpython-38.pyc b/scGCN/__pycache__/models.cpython-38.pyc
new file mode 100644
index 0000000..a8580e8
Binary files /dev/null and b/scGCN/__pycache__/models.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/sankey.cpython-38.pyc b/scGCN/__pycache__/sankey.cpython-38.pyc
new file mode 100644
index 0000000..ae7db46
Binary files /dev/null and b/scGCN/__pycache__/sankey.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/utility.cpython-38.pyc b/scGCN/__pycache__/utility.cpython-38.pyc
new file mode 100644
index 0000000..0450c4b
Binary files /dev/null and b/scGCN/__pycache__/utility.cpython-38.pyc differ
diff --git a/scGCN/__pycache__/utils.cpython-38.pyc b/scGCN/__pycache__/utils.cpython-38.pyc
new file mode 100644
index 0000000..646a5f7
Binary files /dev/null and b/scGCN/__pycache__/utils.cpython-38.pyc differ
diff --git a/scGCN/checkpoints/best_validation.data-00000-of-00001 b/scGCN/checkpoints/best_validation.data-00000-of-00001
new file mode 100644
index 0000000..f655b60
Binary files /dev/null and b/scGCN/checkpoints/best_validation.data-00000-of-00001 differ
diff --git a/scGCN/checkpoints/best_validation.index b/scGCN/checkpoints/best_validation.index
new file mode 100644
index 0000000..5a01eed
Binary files /dev/null and b/scGCN/checkpoints/best_validation.index differ
diff --git a/scGCN/checkpoints/best_validation.meta b/scGCN/checkpoints/best_validation.meta
new file mode 100644
index 0000000..cc86f3b
Binary files /dev/null and b/scGCN/checkpoints/best_validation.meta differ
diff --git a/scGCN/checkpoints/checkpoint b/scGCN/checkpoints/checkpoint
new file mode 100644
index 0000000..fce2a8c
--- /dev/null
+++ b/scGCN/checkpoints/checkpoint
@@ -0,0 +1,2 @@
+model_checkpoint_path: "best_validation"
+all_model_checkpoint_paths: "best_validation"
diff --git a/scGCN/data.py b/scGCN/data.py
index a1e1163..705f3c9 100644
--- a/scGCN/data.py
+++ b/scGCN/data.py
@@ -29,7 +29,7 @@ def input_data(DataDir,Rgraph=True):
lab_data2 = data2.reset_index(drop=True) #.transpose()
lab_label1.columns = ['type']
lab_label2.columns = ['type']
-
+
types = np.unique(lab_label1['type']).tolist()
random.seed(123)
@@ -83,14 +83,15 @@ def input_data(DataDir,Rgraph=True):
lab_train2 = pd.concat([label_train1, lab_label2])
#' save objects
-
+ types_all = np.unique([*lab_label1['type'],*lab_label2['type']]).tolist()
PIK = "{}/datasets.dat".format(DataDir)
res = [
data_train1, data_test1, data_val1, label_train1, label_test1,
- label_val1, lab_data2, lab_label2, types
+ label_val1, lab_data2, lab_label2, types_all
]
with open(PIK, "wb") as f:
pkl.dump(res, f)
-
+ with open("data.txt","w") as out:
+ out.write(str(res))
print('load data succesfully....')
diff --git a/scGCN/data.txt b/scGCN/data.txt
new file mode 100644
index 0000000..8f34b83
--- /dev/null
+++ b/scGCN/data.txt
@@ -0,0 +1,105 @@
+[ ABCC3 ABI3 ACADVL ACKR3 ACTB ... PEPD VIL1 TMEM37 CA2 PGRMC2
+722 0 0 1 0 24 ... 0 4 0 0 2
+241 0 0 0 0 22 ... 0 0 0 0 3
+5925 0 0 0 0 3 ... 0 0 0 0 0
+2037 0 0 1 0 9 ... 0 0 0 0 0
+3496 0 0 2 0 24 ... 0 0 0 0 2
+... ... ... ... ... ... ... ... ... ... ... ...
+1052 0 0 2 0 0 ... 0 0 0 0 0
+4638 0 0 0 0 0 ... 0 0 0 0 0
+2228 0 0 0 0 18 ... 1 0 0 0 0
+402 0 0 0 0 4 ... 0 0 0 0 0
+975 0 0 0 0 14 ... 0 0 0 1 2
+
+[5877 rows x 2000 columns], ABCC3 ABI3 ACADVL ACKR3 ACTB ... PEPD VIL1 TMEM37 CA2 PGRMC2
+4436 0 0 2 0 1 ... 0 0 0 0 0
+5086 0 0 0 0 14 ... 1 0 0 0 2
+386 0 0 1 0 1 ... 0 0 0 0 1
+4657 0 0 0 0 5 ... 0 0 0 0 0
+4323 0 0 3 1 7 ... 0 0 0 0 0
+... ... ... ... ... ... ... ... ... ... ... ...
+2806 0 0 2 0 7 ... 0 0 0 0 1
+1073 0 0 0 0 1 ... 0 0 0 0 0
+1081 0 0 2 0 2 ... 1 0 0 0 1
+1050 0 0 1 0 3 ... 0 0 0 0 0
+1059 0 0 3 0 1 ... 0 0 0 0 0
+
+[731 rows x 2000 columns], ABCC3 ABI3 ACADVL ACKR3 ACTB ... PEPD VIL1 TMEM37 CA2 PGRMC2
+2786 0 0 1 0 3 ... 0 1 0 0 0
+4371 0 0 0 0 13 ... 0 0 0 0 0
+3963 0 0 0 0 6 ... 0 0 0 0 0
+3690 0 0 0 0 32 ... 0 0 0 0 0
+1993 0 0 0 0 2 ... 0 2 0 0 1
+... ... ... ... ... ... ... ... ... ... ... ...
+230 0 0 4 0 5 ... 0 0 0 0 2
+2789 0 0 1 0 12 ... 0 0 0 0 0
+4452 0 0 1 0 5 ... 0 0 0 0 0
+1058 0 0 0 0 2 ... 0 0 0 0 0
+6569 0 0 0 0 11 ... 0 0 0 1 0
+
+[656 rows x 2000 columns], type
+722 alpha
+241 alpha
+5925 alpha
+2037 alpha
+3496 alpha
+... ...
+1052 quiescent_stellate
+4638 quiescent_stellate
+2228 quiescent_stellate
+402 quiescent_stellate
+975 quiescent_stellate
+
+[5877 rows x 1 columns], type
+4436 alpha
+5086 alpha
+386 alpha
+4657 alpha
+4323 alpha
+... ...
+2806 quiescent_stellate
+1073 quiescent_stellate
+1081 quiescent_stellate
+1050 quiescent_stellate
+1059 quiescent_stellate
+
+[731 rows x 1 columns], type
+2786 alpha
+4371 alpha
+3963 alpha
+3690 alpha
+1993 alpha
+... ...
+230 quiescent_stellate
+2789 quiescent_stellate
+4452 quiescent_stellate
+1058 quiescent_stellate
+6569 quiescent_stellate
+
+[656 rows x 1 columns], ABCC3 ABI3 ACADVL ACKR3 ACTB ... PEPD VIL1 TMEM37 CA2 PGRMC2
+0 0 0 0 0 17 ... 0 0 0 0 10
+1 1 0 0 0 82 ... 1 0 0 0 1
+2 0 0 0 0 8 ... 0 0 0 0 1
+3 0 0 0 0 2 ... 0 0 0 0 1
+4 0 0 0 0 8 ... 0 0 0 0 8
+... ... ... ... ... ... ... ... ... ... ... ...
+1836 0 0 0 0 0 ... 0 0 0 0 3
+1837 0 0 0 0 0 ... 0 0 0 0 1
+1838 0 0 0 0 9 ... 0 0 0 0 0
+1839 0 0 0 0 22 ... 0 0 0 0 0
+1840 0 0 0 0 10 ... 0 0 0 0 0
+
+[1841 rows x 2000 columns], type
+0 beta
+1 ductal
+2 delta
+3 delta
+4 beta
+... ...
+1836 beta
+1837 beta
+1838 ductal
+1839 quiescent_stellate
+1840 ductal
+
+[1841 rows x 1 columns], ['alpha', 'beta', 'delta', 'ductal', 'endothelial', 'gamma', 'macrophage', 'quiescent_stellate']]
\ No newline at end of file
diff --git a/scGCN/data_preprocess.R b/scGCN/data_preprocess.R
index 3f2c41b..e4e1f4e 100644
--- a/scGCN/data_preprocess.R
+++ b/scGCN/data_preprocess.R
@@ -1,3 +1,230 @@
+#' This functions provides the optional way to generate scGCN labels
+#' @param Dat1 reference data; rows are genes and columns are cells
+#' @param Dat2 query data; rows are genes and columns are cells
+#' @param Lab1 label of reference data; rows are genes and columns are cells
+#' @param Lab2 query data; rows are genes and columns are cells
+#' @param inter_graph inter-dataset graph in input folder
+
+#' @return This function returns a list of two statistical scores: entropy score and enrichment score that unknown cells have higher entropy score and lower enrichment score
+#' @examples: first load count.list and label.list from folder "example_data" and run 'data_preprocess.R'
+#' then run 'metrics(dat1,lab1,dat2,lab2,inter_graph,cluster)'
+
+suppressMessages(library(Seurat)); suppressMessages(library(entropy))
+
+
+args=commandArgs(T)
+
+ref_c=read.table(args[1],header=T)
+qur_c=read.table(args[3],header=T)
+count.list=list(ref_c,qur_c)
+ref_t=read.table(args[2],header=T)
+qur_t=read.table(args[4],header=T)
+label.list=list(ref_t,qur_t)
+
+
+#save_processed_data(count.list,label.list)
+
+
+
+metrics <- function(lab1,inter_graph,clusters){
+ inter.graph <- inter_graph; hc <- clusters
+
+ #' ------------ check cluster enrichment in cell types ------------
+ #' Examine in clusters
+
+ #' check enrichment
+ eval1 <- sapply(names(table(hc)),function(name){
+ ps1 <- which(hc==name)
+ foo <- table(lab1[inter_graph[which(inter_graph[,2]%in%ps1),1],1])
+ tem <- table(lab1)/sum(table(lab1))
+ if (length(foo)>0){
+ temp <- unlist(sapply(1:length(foo),function(i){
+ ind <- grep(names(foo)[i],names(tem))
+ r1 <- foo[i]/sum(foo)
+ return (r1/tem[ind])}))
+ return (max(temp)/sum(temp))
+ } else {return (0)}
+ })
+
+ #' check mixtureness
+ eval2 <- sapply(names(table(hc)),function(name){
+ ps1 <- which(hc==name)
+ foo <- table(lab1[inter_graph[which(inter_graph[,2]%in%ps1),1],1])
+ tem <- table(lab1)/sum(table(lab1))
+ if (length(foo)>0){
+ temp <- unlist(sapply(1:length(foo),function(i){
+ ind <- grep(names(foo)[i],names(tem))
+ r1 <- foo[i]/sum(foo)
+ return (r1/tem[ind])}))
+ temp2 <- entropy.empirical(temp,unit='log')
+ return (temp2)
+ } else {return (0) }
+ })
+ return (list(eval1=eval1,eval2=eval2))
+}
+
+
+#' This functions provides the optional way to generate scGCN labels
+#' @param Dat1 reference data; rows are genes and columns are cells
+#' @param Dat2 query data; rows are genes and columns are cells
+#' @param Lab1 label of reference data; rows are genes and columns are cells
+#' @param Lab2 query data; rows are genes and columns are cells
+
+#' @return This function returns files saved in folders "input" & "process_data"
+#' @export: all files are saved in current path
+#' @examples: load count.list and label.list from folder "example_data"
+#' save_processed_data(count.list,label.list)
+
+GenerateGraph <- function(Dat1,Dat2,Lab1,K,check.unknown){
+ object1 <- CreateSeuratObject(counts=Dat1,project = "1",assay = "Data1",
+ min.cells = 0,min.features = 0,
+ names.field = 1,names.delim = "_")
+
+ object2 <- CreateSeuratObject(counts=Dat2,project = "2",assay = "Data2",
+ min.cells = 0,min.features =0,names.field = 1,
+ names.delim = "_")
+
+ objects <- list(object1,object2)
+ objects1 <- lapply(objects,function(obj){
+ obj <- NormalizeData(obj,verbose=F)
+ obj <- FindVariableFeatures(obj,
+ selection.method = "vst",
+ nfeatures = 2000,verbose=F)
+ obj <- ScaleData(obj,features=rownames(obj),verbose=FALSE)
+ obj <- RunPCA(obj, features=rownames(obj), verbose = FALSE)
+ return(obj)})
+ #' Inter-data graph
+ object.nn <- FindIntegrationAnchors(object.list = objects1,k.anchor=K,verbose=F)
+ arc=object.nn@anchors
+ d1.arc1=cbind(arc[arc[,4]==1,1],arc[arc[,4]==1,2],arc[arc[,4]==1,3])
+ grp1=d1.arc1[d1.arc1[,3]>0,1:2]-1
+
+ if (check.unknown){
+ obj <- objects1[[2]]
+ obj <- RunPCA(obj, features = VariableFeatures(object = obj),npcs=30,verbose=F)
+ obj <- FindNeighbors(obj,verbose=F)
+ obj <- FindClusters(obj, resolution = 0.5,verbose=F)
+ hc <- Idents(obj); inter.graph=grp1+1
+ scores <- metrics(lab1=Lab1,inter_graph=inter.graph,clusters=hc)
+ saveRDS(scores,file='./input/statistical_scores.RDS')
+ }
+ #' Intra-data graph
+ d2.list <- list(objects1[[2]],objects1[[2]])
+ d2.nn <- FindIntegrationAnchors(object.list =d2.list,k.anchor=K,verbose=F)
+ d2.arc=d2.nn@anchors
+ d2.arc1=cbind(d2.arc[d2.arc[,4]==1,1],d2.arc[d2.arc[,4]==1,2],d2.arc[d2.arc[,4]==1,3])
+ d2.grp=d2.arc1[d2.arc1[,3]>0,1:2]-1
+ final <- list(inteG=grp1,intraG=d2.grp)
+ return (final)
+}
+
+#' This function normalize count data matrix from refrence and query set
+#' @param count.list list of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @return This function returns normalized data list
+#' @export
+#' @examples
+#' normalize_data(count.list)
+normalize_data <- function(count.list){
+ norm.list <- vector('list')
+ var.features <- vector('list')
+ for ( i in 1:length(count.list)){
+ norm.list[[i]] <- as.matrix(Seurat:::NormalizeData.default(count.list[[i]],verbose=F))
+ #' select variable features
+ hvf.info <- Seurat:::FindVariableFeatures.default(count.list[[i]],selection.method='vst',verbose=F)
+ hvf.info <- hvf.info[which(x = hvf.info[, 1, drop = TRUE] != 0), ]
+ hvf.info <- hvf.info[order(hvf.info$vst.variance.standardized, decreasing = TRUE), , drop = FALSE]
+ var.features[[i]] <- head(rownames(hvf.info), n = 2000)
+ }
+ #' select variable features
+ sel.features <- selectHVFeature(count.list,var.features)
+ return (list(norm.list,sel.features))}
+
+#' This function returns scaled data matrix from refrence and query set
+#' @param count.list list of raw counts of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param norm.list list of normalized (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param hvg.features variable features of data list
+#' @return This function returns scaled data list
+#' @export
+#' @examples
+#' scale_data(count.list,norm.list,hvg.features)
+scale_data <- function(count.list,norm.list,hvg.features){
+ scale.list <- lapply(norm.list,function(mat){
+ Seurat:::ScaleData.default(object = mat, features = hvg.features,verbose=F)})
+ scale.list <- lapply(1:length(count.list),function(i){
+ return (scale.list[[i]][na.omit(match(rownames(count.list[[i]]),rownames(scale.list[[i]]))),])})
+ return (scale.list)}
+
+#' This function returns a set of highly variable features
+#' @param count.list list of raw counts of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param var.features variable features of data list
+#' @return This function returns a common set of highly variable features
+#' @export
+#' @examples
+#' selectHVFeature(count.list,var.features)
+selectHVFeature <- function(count.list,var.features,nfeatures = 2000){
+ var.features1 <- unname(unlist(var.features))
+ var.features2 <- sort(table(var.features1), decreasing = TRUE)
+ for (i in 1:length(count.list)) {
+ var.features3 <- var.features2[names(var.features2) %in% rownames(count.list[[i]])]}
+ tie.val <- var.features3[min(nfeatures, length(var.features3))]
+ features <- names(var.features3[which(var.features3 > tie.val)])
+ if (length(features) > 0) {
+ feature.ranks <- sapply(features, function(x) {
+ ranks <- sapply(var.features, function(y) {
+ if (x %in% y) {
+ return(which(x == y))
+ }
+ return(NULL)
+ })
+ median(unlist(ranks))
+ })
+ features <- names(sort(feature.ranks))
+ }
+ features.tie <- var.features3[which(var.features3 == tie.val)]
+ tie.ranks <- sapply(names(features.tie), function(x) {
+ ranks <- sapply(var.features, function(y) {
+ if (x %in% y) {return(which(x == y))}
+ return(NULL)
+ })
+ median(unlist(ranks))
+ })
+ features <- c(features, names(head(sort(tie.ranks), nfeatures - length(features))))
+ return(features)
+}
+
+#' This functions takes refrence data and labels to identify variable gene features
+#' @param data reference data; rows are genes and columns are cells
+#' @param label data frame with rownames identical with colnames of data; the first column is cell type
+#' @param nf number of variable features
+#' @return This function returns variable gene features using ANOVA
+#' @export
+#' @examples
+#' select_feature(data=reference.data,label=reference.label)
+select_feature <- function(data,label,nf=2000){
+ M <- nrow(data); new.label <- label[,1]
+ pv1 <- sapply(1:M, function(i){
+ mydataframe <- data.frame(y=as.numeric(data[i,]), ig=new.label)
+ fit <- aov(y ~ ig, data=mydataframe)
+ summary(fit)[[1]][["Pr(>F)"]][1]})
+ names(pv1) <- rownames(data)
+ pv1.sig <- names(pv1)[order(pv1)[1:nf]]
+ egen <- unique(pv1.sig)
+ return (egen)
+}
+
+#' This functions takes refrence data and labels to identify variable gene features
+#' @param count_list list of reference data and query data; rows are genes and columns are cells
+#' @param label_list list of reference label and query label (if any), both are data frames with rownames identical with colnames of data; the first column is cell type
+#' @return This function returns four elements, including the normalized data, scaled data, hvg features, and selected features
+#' @export
+#' @examples
+#' pre_process(count_list,label_list)
+pre_process <- function(count_list,label_list){
+ sel.features <- select_feature(count_list[[1]],label_list[[1]])
+ count_list_new <- list(count_list[[1]][sel.features,],count_list[[2]][sel.features,])
+ return (count_list_new)
+}
+
#' This functions takes raw counts and labels of reference/query set to generate scGCN training input
#' @param count.list list of reference data and query data; rows are genes and columns are cells
#' @param label.list list of reference label and query label (if any), both are data frames with rownames identical with colnames of data; the first column is cell type
@@ -5,11 +232,65 @@
#' @export: all files are saved in current path
#' @examples: load count.list and label.list from folder "example_data"
#' save_processed_data(count.list,label.list)
+save_processed_data <- function(count.list,label.list,Rgraph=TRUE,check_unknown=FALSE,outdir="input"){
+ count.list <- pre_process(count_list=count.list,label_list=label.list)
+ #' save counts data to certain path: 'input'
+ dir.create(outdir);
+ write.csv(t(count.list[[1]]),file='input/Data1.csv',quote=F,row.names=T)
+ write.csv(t(count.list[[2]]),file='input/Data2.csv',quote=F,row.names=T)
-source('data_preprocess_utility.R')
-
-count.list <- readRDS('example_data/count.list.RDS')
-label.list <- readRDS('example_data/label.list.RDS')
-
+ #' optional graph: R genreated graph has minor differnce with python, user can choose the one with better performance
+ if (Rgraph){
+ #' use R generated graph
+ new.dat1 <- count.list[[1]]; new.dat2 <- count.list[[2]]
+ new.lab1 <- label.list[[1]]; new.lab2 <- label.list[[2]]
+ graphs <- suppressWarnings(GenerateGraph(Dat1=new.dat1,Dat2=new.dat2,
+ Lab1=new.lab1,K=5,
+ check.unknown=check_unknown))
+ write.csv(graphs[[1]],file='input/inter_graph.csv',quote=F,row.names=T)
+ write.csv(graphs[[2]],file='input/intra_graph.csv',quote=F,row.names=T)
+ write.csv(new.lab1,file='input/Label1.csv',quote=F,row.names=F)
+ write.csv(new.lab2,file='input/Label2.csv',quote=F,row.names=F)
+ } else {
+ #' use python generated graph
+ dir.create('results')
+ #' @param norm.list normalized data
+ res1 <- normalize_data(count.list)
+ norm.list <- res1[[1]]; hvg.features <- res1[[2]];
+ #' @param scale.list scaled data
+ scale.list <- scale_data(count.list,norm.list,hvg.features)
+ outputdir <- 'process_data'; dir.create(outputdir)
+ write.csv(hvg.features,file=paste0(outputdir,'/sel_features.csv'),quote=F,row.names=F)
+ N <- length(count.list)
+ for (i in 1:N){
+ df = count.list[[i]]
+ if (!dir.exists(paste0(outputdir,'/count_data'))){dir.create(paste0(outputdir,'/count_data'))}
+ file.name=paste0(outputdir,'/count_data/count_data_',i,'.csv')
+ write.csv(df,file=file.name,quote=F)
+ }
+
+ for (i in 1:N){
+ df = label.list[[i]]
+ if (!dir.exists(paste0(outputdir,'/label_data'))){dir.create(paste0(outputdir,'/label_data'))}
+ file.name=paste0(outputdir,'/label_data/label_data_',i,'.csv')
+ write.csv(df,file=file.name,quote=F)
+ }
+
+ for (i in 1:N){
+ df = norm.list[[i]]
+ if (!dir.exists(paste0(outputdir,'/norm_data'))){dir.create(paste0(outputdir,'/norm_data'))}
+ file.name=paste0(outputdir,'/norm_data/norm_data_',i,'.csv')
+ write.csv(df,file=file.name,quote=F)
+ }
+
+ for (i in 1:N){
+ df = scale.list[[i]]
+ if (!dir.exists(paste0(outputdir,'/scale_data'))){dir.create(paste0(outputdir,'/scale_data'))}
+ file.name=paste0(outputdir,'/scale_data/scale_data_',i,'.csv')
+ write.csv(df,file=file.name,quote=F)
+ }
+ }
+}
save_processed_data(count.list,label.list)
+
diff --git a/scGCN/data_preprocess.raw.R b/scGCN/data_preprocess.raw.R
new file mode 100644
index 0000000..9fb4716
--- /dev/null
+++ b/scGCN/data_preprocess.raw.R
@@ -0,0 +1,25 @@
+#' This functions takes raw counts and labels of reference/query set to generate scGCN training input
+#' @param count.list list of reference data and query data; rows are genes and columns are cells
+#' @param label.list list of reference label and query label (if any), both are data frames with rownames identical with colnames of data; the first column is cell type
+#' @return This function returns files saved in folders "input" & "process_data"
+#' @export: all files are saved in current path
+#' @examples: load count.list and label.list from folder "example_data"
+#' save_processed_data(count.list,label.list)
+
+source('data_preprocess_utility.R')
+
+count.list <- readRDS('example_data/count.list.RDS')
+label.list <- readRDS('example_data/label.list.RDS')
+
+save_processed_data(count.list,label.list)
+#四个文件，一个参考基因表达矩阵，一个待测基因表达矩阵，一个参考标签数据框，一个待测标签数据框
+#> source('data_preprocess_utility.R')
+#> ref_c=read.table("参考基因表达矩阵",header=T)
+#> qur_c=read.table("待测基因表达矩阵",header=T)
+#> count.list=list(ref_c,qur_c)
+#> ref_t=read.table("ref.class.txt参考的标签数据框",header=T)
+#> qur_t=read.table("qur.class.txt待测的标签数据框",header=T)
+#> label.list=list(ref_t,qur_t)
+#> save_processed_data(count.list,label.list)
+
+
diff --git a/scGCN/data_preprocess_utility.R b/scGCN/data_preprocess_utility.R
index c206395..0d02f45 100644
--- a/scGCN/data_preprocess_utility.R
+++ b/scGCN/data_preprocess_utility.R
@@ -217,10 +217,10 @@ pre_process <- function(count_list,label_list){
#' @export: all files are saved in current path
#' @examples: load count.list and label.list from folder "example_data"
#' save_processed_data(count.list,label.list)
-save_processed_data <- function(count.list,label.list,Rgraph=TRUE,check_unknown=FALSE){
+save_processed_data <- function(count.list,label.list,Rgraph=TRUE,check_unknown=FALSE,outdir="input"){
count.list <- pre_process(count_list=count.list,label_list=label.list)
#' save counts data to certain path: 'input'
- dir.create('input');
+ dir.create(outdir);
write.csv(t(count.list[[1]]),file='input/Data1.csv',quote=F,row.names=T)
write.csv(t(count.list[[2]]),file='input/Data2.csv',quote=F,row.names=T)
@@ -234,8 +234,15 @@ save_processed_data <- function(count.list,label.list,Rgraph=TRUE,check_unknown=
check.unknown=check_unknown))
write.csv(graphs[[1]],file='input/inter_graph.csv',quote=F,row.names=T)
write.csv(graphs[[2]],file='input/intra_graph.csv',quote=F,row.names=T)
- write.csv(new.lab1,file='input/label1.csv',quote=F,row.names=F)
- write.csv(new.lab2,file='input/label2.csv',quote=F,row.names=F)
+<<<<<<< Updated upstream
+ write.csv(new.lab1,file='input/Label1.csv',quote=F,row.names=F)
+ write.csv(new.lab2,file='input/Label2.csv',quote=F,row.names=F)
+=======
+ #数据集之间锚点
+ write.csv(new.lab1,file='input/Label1.csv',quote=F,row.names=F)
+ write.csv(new.lab2,file='input/Label2.csv',quote=F,row.names=F)
+ #到此未对label2做处理
+>>>>>>> Stashed changes
} else {
#' use python generated graph
dir.create('results')
diff --git a/scGCN/data_preprocess_utility.raw.R b/scGCN/data_preprocess_utility.raw.R
new file mode 100644
index 0000000..8df9789
--- /dev/null
+++ b/scGCN/data_preprocess_utility.raw.R
@@ -0,0 +1,323 @@
+#' This functions provides the optional way to generate scGCN labels
+#' @param Dat1 reference data; rows are genes and columns are cells
+#' @param Dat2 query data; rows are genes and columns are cells
+#' @param Lab1 label of reference data; rows are genes and columns are cells
+#' @param Lab2 query data; rows are genes and columns are cells
+#' @param inter_graph inter-dataset graph in input folder
+
+#' @return This function returns a list of two statistical scores: entropy score and enrichment score that unknown cells have higher entropy score and lower enrichment score
+#熵值和富集值：未知类型的细胞具有较高的熵和较低的富集度；
+#' @examples: first load count.list and label.list from folder "example_data" and run 'data_preprocess.R'
+#' then run 'metrics(dat1,lab1,dat2,lab2,inter_graph,cluster)'
+
+suppressMessages(library(Seurat)); suppressMessages(library(entropy))
+
+metrics <- function(lab1,inter_graph,clusters){
+ #调用的时候是处理待测数据时
+ #scores <- metrics(lab1=Lab1,inter_graph=inter.graph,clusters=hc)
+ #inter_graph和细胞整合锚点有关；clusters和细胞标识有关
+ inter.graph <- inter_graph; hc <- clusters
+
+ #' ------------ check cluster enrichment in cell types ------------
+ #' Examine in clusters
+
+ #' check enrichment 细胞富集度
+ eval1 <- sapply(names(table(hc)),function(name){
+ ps1 <- which(hc==name)
+ foo <- table(lab1[inter_graph[which(inter_graph[,2]%in%ps1),1],1])
+ tem <- table(lab1)/sum(table(lab1))
+ if (length(foo)>0){
+ temp <- unlist(sapply(1:length(foo),function(i){
+ ind <- grep(names(foo)[i],names(tem))
+ r1 <- foo[i]/sum(foo)
+ return (r1/tem[ind])}))
+ return (max(temp)/sum(temp))
+ } else {return (0)}
+ })
+
+ #' check mixtureness
+ eval2 <- sapply(names(table(hc)),function(name){
+ ps1 <- which(hc==name)
+ foo <- table(lab1[inter_graph[which(inter_graph[,2]%in%ps1),1],1])
+ tem <- table(lab1)/sum(table(lab1))
+ if (length(foo)>0){
+ temp <- unlist(sapply(1:length(foo),function(i){
+ ind <- grep(names(foo)[i],names(tem))
+ r1 <- foo[i]/sum(foo)
+ return (r1/tem[ind])}))
+ #以上和之前一样？
+ temp2 <- entropy.empirical(temp,unit='log')
+ return (temp2)
+ } else {return (0) }
+ })
+ return (list(eval1=eval1,eval2=eval2))
+}
+
+
+#' This functions provides the optional way to generate scGCN labels
+#' @param Dat1 reference data; rows are genes and columns are cells
+#' @param Dat2 query data; rows are genes and columns are cells
+#' @param Lab1 label of reference data; rows are genes and columns are cells
+#' @param Lab2 query data; rows are genes and columns are cells
+
+#' @return This function returns files saved in folders "input" & "process_data"
+#' @export: all files are saved in current path
+#' @examples: load count.list and label.list from folder "example_data"
+#' save_processed_data(count.list,label.list)
+
+GenerateGraph <- function(Dat1,Dat2,Lab1,K,check.unknown){
+ #调用的时候k的赋值为5
+ object1 <- CreateSeuratObject(counts=Dat1,project = "1",assay = "Data1",
+ min.cells = 0,min.features = 0,
+ names.field = 1,names.delim = "_")
+
+ object2 <- CreateSeuratObject(counts=Dat2,project = "2",assay = "Data2",
+ min.cells = 0,min.features =0,names.field = 1,
+ names.delim = "_")
+ #分别对ref和qur创建seurat对象
+ objects <- list(object1,object2)
+ objects1 <- lapply(objects,function(obj){
+ obj <- NormalizeData(obj,verbose=F)
+ obj <- FindVariableFeatures(obj,
+ selection.method = "vst",
+ nfeatures = 2000,verbose=F)
+ #这里也有一个基因数目的处理2000
+ obj <- ScaleData(obj,features=rownames(obj),verbose=FALSE)
+ obj <- RunPCA(obj, features=rownames(obj), verbose = FALSE)
+ return(obj)})
+ #用seurat做一些常规的单细胞数据预处理
+
+ #' Inter-data graph
+ object.nn <- FindIntegrationAnchors(object.list = objects1,k.anchor=K,verbose=F)
+ #寻找锚点：Find a set of anchors between a list of Seurat objects, K=5
+ #These anchors can later be used to integrate the objects using the IntegrateData function
+ #k.anchor:How many neighbors (k) to use when picking anchors
+ arc=object.nn@anchors
+ #将整体细胞的锚点信息赋值给变量
+ d1.arc1=cbind(arc[arc[,4]==1,1],arc[arc[,4]==1,2],arc[arc[,4]==1,3])
+ grp1=d1.arc1[d1.arc1[,3]>0,1:2]-1
+ #这里是？个性化的处理；最终得到的grp1是一个数据框，第一列是cell1的标识-1，第二列是cell2的标识-1，保留标准是score>0
+
+ if (check.unknown){
+ obj <- objects1[[2]]
+ obj <- RunPCA(obj, features = VariableFeatures(object = obj),npcs=30,verbose=F)
+ obj <- FindNeighbors(obj,verbose=F)
+ #Computes the k.param nearest neighbors for a given dataset. Can also optionally (via compute.SNN),
+ #construct a shared nearest neighbor graph by calculating the neighborhood overlap (Jaccard index)
+ # between every cell and its k.param nearest neighbors.
+ obj <- FindClusters(obj, resolution = 0.5,verbose=F)
+ hc <- Idents(obj); inter.graph=grp1+1
+ #hc是用Idents函数获取细胞标识，那么seurat对象的标识是？
+ #inter.graph和细胞整合的锚点有关
+ scores <- metrics(lab1=Lab1,inter_graph=inter.graph,clusters=hc)
+ #这里是对待测数据处理，用的label1？
+ saveRDS(scores,file='./input/statistical_scores.RDS')
+ #检验待测细胞类型的熵和富集度
+ }
+ #' Intra-data graph
+ d2.list <- list(objects1[[2]],objects1[[2]])
+ d2.nn <- FindIntegrationAnchors(object.list =d2.list,k.anchor=K,verbose=F)
+
+ d2.arc=d2.nn@anchors
+ d2.arc1=cbind(d2.arc[d2.arc[,4]==1,1],d2.arc[d2.arc[,4]==1,2],d2.arc[d2.arc[,4]==1,3])
+ d2.grp=d2.arc1[d2.arc1[,3]>0,1:2]-1
+ #个性化操作锚点
+ final <- list(inteG=grp1,intraG=d2.grp)
+ #返回两个有关锚点的值;一个是参考数据集和待测数据集之间的相似度，一个是数据集内部的相似度
+ return (final)
+}
+
+#' This function normalize count data matrix from refrence and query set
+#' @param count.list list of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @return This function returns normalized data list
+#' @export
+#' @examples
+#' normalize_data(count.list)
+normalize_data <- function(count.list){
+ norm.list <- vector('list')
+ var.features <- vector('list')
+ for ( i in 1:length(count.list)){
+ norm.list[[i]] <- as.matrix(Seurat:::NormalizeData.default(count.list[[i]],verbose=F))
+ #' select variable features
+ hvf.info <- Seurat:::FindVariableFeatures.default(count.list[[i]],selection.method='vst',verbose=F)
+ hvf.info <- hvf.info[which(x = hvf.info[, 1, drop = TRUE] != 0), ]
+ hvf.info <- hvf.info[order(hvf.info$vst.variance.standardized, decreasing = TRUE), , drop = FALSE]
+ var.features[[i]] <- head(rownames(hvf.info), n = 2000)
+ }
+ #' select variable features
+ sel.features <- selectHVFeature(count.list,var.features)
+ return (list(norm.list,sel.features))}
+
+#' This function returns scaled data matrix from refrence and query set
+#' @param count.list list of raw counts of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param norm.list list of normalized (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param hvg.features variable features of data list
+#' @return This function returns scaled data list
+#' @export
+#' @examples
+#' scale_data(count.list,norm.list,hvg.features)
+scale_data <- function(count.list,norm.list,hvg.features){
+ scale.list <- lapply(norm.list,function(mat){
+ Seurat:::ScaleData.default(object = mat, features = hvg.features,verbose=F)})
+ scale.list <- lapply(1:length(count.list),function(i){
+ return (scale.list[[i]][na.omit(match(rownames(count.list[[i]]),rownames(scale.list[[i]]))),])})
+ return (scale.list)}
+
+#' This function returns a set of highly variable features
+#' @param count.list list of raw counts of (1) reference data and (2) query data, rows are genes and columns are cells
+#' @param var.features variable features of data list
+#' @return This function returns a common set of highly variable features
+#' @export
+#' @examples
+#' selectHVFeature(count.list,var.features)
+selectHVFeature <- function(count.list,var.features,nfeatures = 2000){
+ var.features1 <- unname(unlist(var.features))
+ var.features2 <- sort(table(var.features1), decreasing = TRUE)
+ for (i in 1:length(count.list)) {
+ var.features3 <- var.features2[names(var.features2) %in% rownames(count.list[[i]])]}
+ #确保基因既在ref中也在qur中
+ tie.val <- var.features3[min(nfeatures, length(var.features3))]
+ #反正只保留2000个基因数目
+ features <- names(var.features3[which(var.features3 > tie.val)])
+ if (length(features) > 0) {
+ feature.ranks <- sapply(features, function(x) {
+ ranks <- sapply(var.features, function(y) {
+ if (x %in% y) {
+ return(which(x == y))
+ }
+ return(NULL)
+ })
+ median(unlist(ranks))
+ })
+ features <- names(sort(feature.ranks))
+ }
+ features.tie <- var.features3[which(var.features3 == tie.val)]
+ tie.ranks <- sapply(names(features.tie), function(x) {
+ ranks <- sapply(var.features, function(y) {
+ if (x %in% y) {return(which(x == y))}
+ return(NULL)
+ })
+ median(unlist(ranks))
+ })
+ features <- c(features, names(head(sort(tie.ranks), nfeatures - length(features))))
+ return(features)
+}
+
+#' This functions takes refrence data and labels to identify variable gene features
+#筛选高变基因的过程
+#' @param data reference data; rows are genes and columns are cells
+#' @param label data frame with rownames identical with colnames of data; the first column is cell type
+#' @param nf number of variable features
+#' @return This function returns variable gene features using ANOVA
+#' @export
+#' @examples
+#' select_feature(data=reference.data,label=reference.label)
+select_feature <- function(data,label,nf=2000){
+ M <- nrow(data); new.label <- label[,1]
+ #M是行数，参考基因表达矩阵的行数；new.label是细胞类型标签
+ pv1 <- sapply(1:M, function(i){
+ mydataframe <- data.frame(y=as.numeric(data[i,]), ig=new.label)
+ fit <- aov(y ~ ig, data=mydataframe)
+ #aov函数是在做方差分析；分析同一个基因的不同细胞的关系？
+ summary(fit)[[1]][["Pr(>F)"]][1]})
+ names(pv1) <- rownames(data)
+ #rownames是基因名？
+ pv1.sig <- names(pv1)[order(pv1)[1:nf]]
+ #nf是挑选的高变基因的数目，在此是2000个
+ egen <- unique(pv1.sig)
+ return (egen)
+ #返回的是高变基因特征
+}
+
+#' This functions takes refrence data and labels to identify variable gene features
+#' @param count_list list of reference data and query data; rows are genes and columns are cells
+#' @param label_list list of reference label and query label (if any), both are data frames with rownames identical with colnames of data; the first column is cell type
+#' @return This function returns four elements, including the normalized data, scaled data, hvg features, and selected features
+#' @export
+#' @examples
+#' pre_process(count_list,label_list)
+pre_process <- function(count_list,label_list){
+ sel.features <- select_feature(count_list[[1]],label_list[[1]])
+ #提取参考基因表达矩阵和参考的细胞标签；进行高变基因的挑选
+ count_list_new <- list(count_list[[1]][sel.features,],count_list[[2]][sel.features,])
+ #对参考基因表达矩阵做了处理，保留2000个高变基因的数据
+ return (count_list_new)
+}
+
+#' This functions takes raw counts and labels of reference/query set to generate scGCN training input
+#' @param count.list list of reference data and query data; rows are genes and columns are cells
+#' @param label.list list of reference label and query label (if any), both are data frames with rownames identical with colnames of data; the first column is cell type
+#' @return This function returns files saved in folders "input" & "process_data"
+#' @export: all files are saved in current path
+#' @examples: load count.list and label.list from folder "example_data"
+#' save_processed_data(count.list,label.list)
+save_processed_data <- function(count.list,label.list,Rgraph=TRUE,check_unknown=FALSE){
+ count.list <- pre_process(count_list=count.list,label_list=label.list)
+ #对参考数据处理，这一步得到的是做了2000个高变基因提取的基因表达矩阵和标签
+ #' save counts data to certain path: 'input'
+ dir.create('input');
+ write.csv(t(count.list[[1]]),file='input/Data1.csv',quote=F,row.names=T)
+ write.csv(t(count.list[[2]]),file='input/Data2.csv',quote=F,row.names=T)
+ #分别保存两个基因表达矩阵：（1）对参考做了2000个高变基因处理（2）对待测未处理
+
+ #' optional graph: R genreated graph has minor differnce with python, user can choose the one with better performance
+ #使用R或者python两种方法来生成graph
+ if (Rgraph){
+ #' use R generated graph
+ new.dat1 <- count.list[[1]]; new.dat2 <- count.list[[2]]
+ new.lab1 <- label.list[[1]]; new.lab2 <- label.list[[2]]
+ graphs <- suppressWarnings(GenerateGraph(Dat1=new.dat1,Dat2=new.dat2,
+ Lab1=new.lab1,K=5,
+ check.unknown=check_unknown))
+ #返回两个锚点信息
+ write.csv(graphs[[1]],file='input/inter_graph.csv',quote=F,row.names=T)
+ #数据集内部锚点
+ write.csv(graphs[[2]],file='input/intra_graph.csv',quote=F,row.names=T)
+ #数据集之间锚点
+ write.csv(new.lab1,file='input/label1.csv',quote=F,row.names=F)
+ write.csv(new.lab2,file='input/label2.csv',quote=F,row.names=F)
+ #到此未对label2做处理
+ } else {
+ #' use python generated graph
+ dir.create('results')
+ #' @param norm.list normalized data
+ res1 <- normalize_data(count.list)
+ #得到标准化表达矩阵以及高变基因集
+ norm.list <- res1[[1]]; hvg.features <- res1[[2]];
+
+ #' @param scale.list scaled data
+ scale.list <- scale_data(count.list,norm.list,hvg.features)
+ outputdir <- 'process_data'; dir.create(outputdir)
+ write.csv(hvg.features,file=paste0(outputdir,'/sel_features.csv'),quote=F,row.names=F)
+ N <