前言

从两个月前老板说要搞一个简单的生信分析到现在，经过一个转录组分析的训练后，感觉可以试试看了。这个项目的要求是从一个单细胞数据里比较不同细胞在不同通路上转录水平的差异。同样的，对代码只是一知半解，留下了一些坑，下学期闲下来得开始系统地学学R和Linux了。唔，文献阅读已经鸽了两周了，咕咕咕。还需要交代一下大背景，我直接用我电脑上的Rstudio来完成的，我电脑8G内存，处理这个项目略显吃力。

致谢&参考资料

这次的参考主要来自@刘小泽和@生信start_site的相关教程，深表谢意。

CH0数据准备

回过头来看，这一部分居然是我碰到的最大困难。其实文献我早就挑好了，Single-cell transcriptomic profiling of the aging mouse brain。这篇文章提供了一个可视化的网站，看了几个我关心的通路基因，感觉没啥差异，算了，还是老老实实自己弄吧。
我先尝试像上回做转录组一样去GEO下了SRR的数据，那数据太大了，而且不知道为啥我的aspera在高速下载完一个后又会切换回http，下回还是用循环语句下载吧hhh。
看到其他教程都是下载csv格式的原始矩阵作为起点，但是到我这变成了变成了这个：

不管先下载下来解压看看，大概是这样，不知道怎么往下处理了。

于是我开始找有没有带csv格式的数据，找了一圈也没合适的。就只好回到这个不知道是啥的数据，感谢Jimmy老师的文章救了我，原来这个似乎是一个更高级的矩阵。只需要在R中输入最简单的读取表格的代码就可以了。

##我先用年轻小鼠的8个数据开始分析，用一个sapply函数批量读取表格，建立一个list
memory.limit(8000000) #先修改一下内存限制，我这直接改到了8G，否则内存不够直接报错
path = ".\\GSE129788_RAW" 
fileNames <- dir(path)
filePath <- sapply(fileNames, function(x){paste(path,x,sep='\\')})
data <- lapply(filePath, function(x){read.table(x, header=T)})

这样就把数据导入到了R

CH1seurat标准流程

Seurat大概是用的最广泛的10XscRNA-seq处理的包了，原文的数据也是基于这个包处理的。对于count matrix表达矩阵，可以直接转变成Seurat对象。批量处理算法还要学习一下，我就只好一个一个转了，还好就8个。

library(Seurat)
count1 <- CreateSeuratObject(counts = data[[1]])
count2 <- CreateSeuratObject(counts = data[[2]])
count3 <- CreateSeuratObject(counts = data[[3]])
count4 <- CreateSeuratObject(counts = data[[4]])
count5 <- CreateSeuratObject(counts = data[[5]])
count6 <- CreateSeuratObject(counts = data[[6]])
count7 <- CreateSeuratObject(counts = data[[8]])
count7 <- CreateSeuratObject(counts = data[[7]])
count8 <- CreateSeuratObject(counts = data[[8]])

##做一个合并
ids <- c("y1", "y2", "y3", "y4", "y5", "y6", "y7", "y8") #名命一下
sce.big <- merge(count1,
y = c(count2, count3, count4, count5, count6, count7, count8),
add.cell.ids = ids,
project = "brain")

##存一下
save(sce.big,file = 'sce.big_merge_young.brain.Rdata')

2021.6.8更新：发现了一个优雅的函数可以合并掉数据框里的数据

data_merge <- do.call(cbind, data)

这样就可以将一个数据框转成一个计数矩阵

下面就开始Seurat标准流程

sce.big <- NormalizeData(sce.big)
sce.big <- FindVariableFeatures(sce.big, selection.method = "vst", nfeatures = 2000)
all.genes <- row.names(sce.big)
sce.big <- ScaleData(sce.big, features = all.genes)
sce.big <- RunPCA(sce.big, features = VariableFeatures(object = sce.big))
sce.big <- FindNeighbors(sce.big, dims = 1:10)
sce.big_cluster <- FindClusters(sce.big, resolution = 0.6) #注意这里的resolution是分辨率，越大最后得到的群越多，可以根据需求调整
sce.big_cluster <- RunTSNE(object = sce.big_cluster, dims = 1:10)

下面就可以画一个聚类结果的图了

DimPlot(sce.big_cluster, label = T)

CH2数据注释

下面我们就来看一些marker gene的表达情况

markergens <- c("Gfap","Gja1","Syt1","Tubb3","Csf1r","Cx3cr1", "Tmem119")
markergens
FeaturePlot(object = sce.big_cluster, features = markergens, cols = c("grey", "blue"))

我的后续分析只需要三种细胞，所以就注释了三种细胞

new.cluster.ids <- c("astrocyte", "1", "microglia", "3", "neuron","5", "6", "neuron", "8", "9", "10", "astrocyte", "12", "13", "14", "15", "16")
names(new.cluster.ids) <- levels(sce.big_cluster)
sce.big_cluster_new <- RenameIdents(sce.big_cluster, new.cluster.ids)
DimPlot(sce.big_cluster_new, label = T)

说一句，如果好几个subgroup都是同一种细胞的话，直接把对应的编号改成一样的名字，它们会自己合并（这个问题困扰了我好久，师姐居然告诉我要调前面那个resolution才能解决）

把这三种细胞挑出来，存一下开始下一步吧

sce.big_for_DEG = SubsetData(sce.big_cluster_new, ident.use = c("neuron", "astrocyte", "microglia"))
save(sce.big_for_DEG, file = "young_mouse_analysis.Rdata")

CH3monocle差异分析

做一些准备工作

count_matrix = sce.big_for_DEG@assays[["RNA"]]@data
cluster <- sce.big_for_DEG@active.ident
gene_annotation <- as.data.frame(rownames(count_matrix))

开始monocle流程

expr_matrix <- as.matrix(count_matrix)

sample_ann <- data.frame(cells = names(as.data.frame(count_matrix)), cellType = cluster)
rownames(sample_ann) <- names(as.data.frame(count_matrix))

gene_ann <- as.data.frame(rownames(count_matrix))
rownames(gene_ann) <- rownames(count_matrix)
colnames(gene_ann) <- “generic.skeleton()”
colnames(gene_ann) <- "genes"

pd <- new("AnnotatedDataFramw", data = sample_ann)
pd <- new("AnnotatedDataFrame",
data=sample_ann)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
sc_cds <- newCellDataSet(
expr_matrix,
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=1)
cds = detectGenes(cds, min_expr = 0.1)
cds = sc_cds
cds <- detectGenes(cds, min_expr = 0.1)

继续对cds对象进行处理

##过滤
expressed_gene <- row.names(subset(cds@featureData@data, num_cells_expressed >=5))
length(expressed_gene)
cds <- cds[expressed_gene]

##筛选
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
disp_table <- dispersionTable(cds)
unsup_clusteringene <- subset(disp_table, mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clusteringene$gene_id)
plot_ordering_genes(cds) #画个图看看

## 降维&聚类
cds <- reduceDimension(cds, max_components = 2, num_dim = 6, reduction_method = 'tSNE', verbose = T)
cds <- clusterCells(cds, num_clusters = 4)
plot_cell_clusters(cds, 1,2, color = "cellType") #再画个图看聚类效果

差异分析

diff_test_res <- differentialGeneTest(cds, fullModelFormulaStr = "~cellType")
sig_gene <- subset(diff_test_res, qval<0.1)

CH4作图

这个项目的目的是比较不同细胞在一个通路中的转录组差异，于是我们需要挑选这个通路上的所有基因画一个热图，具体通路我这保密了。

获得gene list

我们这里用KEGG的数据库导出一个通路上的所有基因，我们需要先获取通路的编号，首先访问KEGG pathway

library(KEGGREST)
gs <- keggGet('mmu*')
gs[[1]]$GENE
genes <- unlist(lapply(gs[[1]]$GENE, function(x) strsplit(x, ";")))
genelist <- genes[1:length(genes)%%3 == 2]

下面开始作图

##比较一下通路基因和差异分析结果
genelist %in% rownames(sig_gene)
#去掉不要的

##画图
library(pheatmap)
dat=count_matrix[htmapGenes,]
pheatmap(dat,show_colnames =F,
         show_rownames = T,
         cluster_cols = F, 
         cluster_rows = T)

最后大概长这样

尾声

后续希望利用这篇文章看看aging的影响，以及R真的要开始系统学了，希望早日闲下来吧

2021.6.8更新：通路富集

时隔一年，这个课题的实验已经做的差不多了，临近毕业终于有时间来学习了，刚好分析这部分的内容还有一些工作要做，就接着做了一个通路富集，看看有没有运气碰到我们关心的通路。这篇文章参考了https://www.jianshu.com/p/e47913f5fe0e中的相关方法。

首先需要确定我们之前找到的差异基因分属于哪一群，这里也参考了一个算法：

#首先包一个函数
function(count, clustering, DE_genes){
    cluster_nb <- unique(clustering)
    mean_per_cluster <- vector()
    DE_genes <- DE_genes[order(rownames(DE_genes)),]
    count <- count[order(rownames(count)),]
    count_de_genes <- count[rownames(count) %in% DE_genes$genes,]
    print(dim(count_de_genes))
    for (clusters in cluster_nb) {
        # print(head(count_de_genes[,
        #       colnames(count_de_genes) %in% names(clustering[clustering==clusters])
        #   ]))
        mean <- rowMeans(
            as.matrix(count_de_genes[,
                                     colnames(count_de_genes) %in% names(clustering[clustering==clusters])
            ])
        )
        names(mean) <- clusters
        mean_per_cluster <- cbind(
            mean_per_cluster,
            mean
        )
    }
    colnames(mean_per_cluster) <- cluster_nb
    up_reg_cluster <- colnames(mean_per_cluster)[apply(mean_per_cluster,1,which.max)]
    de_genes_table <- data.frame(
        DE_genes,
        mean_per_cluster,
        cluster=up_reg_cluster
    )
    
    return(de_genes_table)
}

#接下来用这个函数进行分配，这里需要我们前面准备的三个东西：差异基因、分群信息、计数矩阵
de_cluster <- get_up_reg_clusters(count_matrix,cluster,sig_gene)

然后就是GO分析的部分了：

#先进行基因ID的转换
library(clusterProfiler)
entrez <- bitr(sig_gene$genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
de_gene_clusters <- de_cluster[de_cluster$genes %in% entrez$SYMBOL, c("genes", "cluster")]
de_gene_clusters <- data.frame(ENTREZID=entrez$ENTREZID[entrez$SYMBOL %in% de_gene_clusters$genes],cluster=de_gene_clusters$cluster )
list_de_gene_clusters <- split(de_gene_clusters$ENTREZID, 
                               de_gene_clusters$cluster)

#GO分析
formula_res <- compareCluster(
  ENTREZID~cluster, 
  data=de_gene_clusters, 
  fun="enrichGO", 
  OrgDb="org.Mm.eg.db",
  ont          = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.01,
  qvalueCutoff  = 0.05
)

# 可视化
pdf('DEG_GO_each_cluster.pdf',width = 11,height = 6)
dotplot(formula_res, showCategory=5)
dev.off()

最后得到的图

通路富集