marker基因的表达量可视化

课程笔记

粉丝：有单细胞线上课程吗？

小编：什么？ 我们的单细胞转录组分析线上课程已经上线好久了，你们竟然都不知道吗，每篇推文后面的课程推荐没人看的吗，小编已哭晕在厕所

好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接） 

这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获！

作者 | 单细胞天地小编 刘小泽

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二单元第10讲：marker基因的表达量可视化

前言

这次的任务是模仿原文的：

肿瘤组织的差异分析(Fig. 4a、c)

下载数据

之前的分析都是基于第一个病人的PBMC，这次将基于这位病人的tumor：GSE117988_raw.expMatrix_Tumor.csv.gz

1start_time <- Sys.time()
2raw_dataTumor <- read.csv('./GSE117988_raw.expMatrix_Tumor.csv.gz', header = TRUE, row.names = 1)
3end_time <- Sys.time()
4end_time - start_time
5# Time difference of 47.00589 secs
6dim(raw_dataTumor) # 21861基因，7431细胞- already filtered

注意：之前PBMC包含了四个时间点的样本,这里的Tumor包含了2个样本（治疗之前Pre和复发AR）

常规流程

step1: 归一化

1dataTumor <- log2(1 + sweep(raw_dataTumor, 2, median(colSums(raw_dataTumor))/colSums(raw_dataTumor), '*'))
2> head(colnames(dataTumor))
3[1] "AAACCTGAGGATGTAT.1" "AAACCTGCAGCGATCC.1" "AAACCTGGTACGAAAT.1" "AAACGGGAGCTGGAAC.1" "AAACGGGAGGAGTTGC.1"
4[6] "AAACGGGAGTTTAGGA.1"

step2: 自定义划分细胞类型

1cellTypes <- sapply(colnames(dataTumor), function(x) ExtractField(x, 2, '[.]'))
2cellTypes <-ifelse(cellTypes == '1', 'Tumor_Before', 'Tumor_AcquiredResistance')
3> table(cellTypes)
4cellTypes
5Tumor_AcquiredResistance             Tumor_Before 
6                    5188                     2243 

step3: 表达矩阵质控

1# 第一点：基因在多少细胞表达 
2> fivenum(apply(dataTumor,1,function(x) sum(x>0) ))
3   VP2 GPRIN2   EML3 ZNF140  RPLP1 
4     1      8    103    566   7431
5# 第二点：细胞中有多少表达的基因
6> fivenum(apply(dataTumor,2,function(x) sum(x>0) ))
7GGAACTTAGGAATCGC.1 TACGGTACAAGCCGCT.2 CCTACCAAGCGTGAGT.1 TGAGCCGAGACTAGAT.2 GATCGTAGTCATATGC.2 
8             192.0             1059.0             1380.0             1971.5             5888.0 

看到大部分基因在500多个细胞表达，细胞平均能表达1000个基因以上

step4: 创建Seurat对象

1tumor <- CreateSeuratObject(dataTumor, 
2                           min.cells = 1, min.features = 0, project = '10x_Tumor')
3> tumor 
4An object of class seurat in project 10x_Tumor 
5 21861 genes across 7431 samples.

step5: 添加metadata (nUMI 和 细胞类型)

1tumor <- AddMetaData(object = tumor, metadata = apply(raw_dataTumor, 2, sum), col.name = 'nUMI_raw')
2tumor <- AddMetaData(object = tumor, metadata = cellTypes, col.name = 'cellTypes')

step6:  聚类标准流程

 1start_time <- Sys.time()
 2tumor <- ScaleData(object = tumor, vars.to.regress = c('nUMI_raw'), model.use = 'linear', use.umi = FALSE)
 3tumor <- FindVariableGenes(object = tumor, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
 4tumor <- RunPCA(object = tumor, pc.genes = [e-mail])
 5tumor <- RunTSNE(object = tumor, dims.use = 1:10, perplexity = 25)
 6end_time <- Sys.time()
 7end_time - start_time
 8# Time difference of 3.324982 mins
 9TSNEPlot(tumor, group.by = 'cellTypes', colors.use = c('#EF8A62', '#67A9CF'))
10save(tumor,file = 'patient1.Tumor.V2.output.Rdata') # 3.6Gb 大小

接着进行基因可视化

1rm(list = ls()) 
2options(warn=-1) 
3start_time <- Sys.time()
4load('patient1.Tumor.V2.output.Rdata')
5end_time <- Sys.time()
6end_time - start_time
7# Time difference of 19.83097 secs

取出log归一化后的表达矩阵

1count_matrix=tumor@data
2> count_matrix[1:4,1:4]
3              AAACCTGAGGATGTAT.1 AAACCTGCAGCGATCC.1 AAACCTGGTACGAAAT.1 AAACGGGAGCTGGAAC.1
4VP2                    0.0000000                  0                  0                  0
5largeTAntigen          0.9670525                  0                  0                  0
6smallTAntigen          0.0000000                  0                  0                  0
7RP11-34P13.7           0.0000000                  0                  0                  0

取出细胞分群信息

1cluster=[e-mail]$cellTypes
2> table(cluster)
3cluster
4Tumor_AcquiredResistance             Tumor_Before 
5                    5188                     2243

提取基因信息

文章主要探索了治疗前和复发后的HLA-A和HLA-B的变化，于是我们先看看有没有这两个基因

1allGenes = row.names([e-mail])
2> allGenes[grep('HLA',allGenes)]
3 [1] "HHLA3"    "HLA-F"    "HLA-G"    "HLA-A"    "HLA-E"    "HLA-C"    "HLA-B"    "HLA-DRA"  "HLA-DRB5"
4[10] "HLA-DRB1" "HLA-DQA1" "HLA-DQB1" "HLA-DQA2" "HLA-DQB2" "HLA-DOB"  "HLA-DMB"  "HLA-DMA"  "HLA-DOA" 
5[19] "HLA-DPA1" "HLA-DPB1"

对HLA-A操作

1FeaturePlot(object = tumor, 
2            features.plot ='HLA-A', 
3            cols.use = c("grey", "blue"), 
4            reduction.use = "tsne")

看HLA-A表达量

1> table(count_matrix['HLA-A',]>0, cluster)
2       cluster
3        Tumor_AcquiredResistance Tumor_Before
4  FALSE                     2282         1057
5  TRUE                      2906         1186

对HLA-B操作

1FeaturePlot(object = tumor, 
2            features.plot ='HLA-B', 
3            cols.use = c("grey", "blue"), 
4            reduction.use = "tsne")

可以看到治疗前HLA-A基因有1186个表达，1057个不表达；复发后这个基因表达和不表达的数量也相近

1> table(count_matrix['HLA-B',]>0, cluster)
2       cluster
3        Tumor_AcquiredResistance Tumor_Before
4  FALSE                     4794         1258
5  TRUE                       394          985

对HLA-B来讲，不管是治疗前还是复发后，它的表达和不表达差异就很明显。另外从治疗前到复发，这个基因的表达数量的变化更显著

小问题：HLA基因这么多，我们怎么找到全部的具有相似模式的基因？

所谓相似表达模式，就是像HLA-B基因一样，在一个群表达很多，另一个群表达很少，通过卡方检验（http://rpubs.com/chixinzero/490992）就能看出区别

 1>  chisq.test(table(count_matrix['HLA-A',]>0, cluster))
 2
 3    Pearson's Chi-squared test with Yates' continuity correction
 4
 5data:  table(count_matrix["HLA-A", ] > 0, cluster)
 6X-squared = 6.1069, df = 1, p-value = 0.01347
 7
 8>  chisq.test(table(count_matrix['HLA-B',]>0, cluster))
 9
10    Pearson's Chi-squared test with Yates' continuity correction
11
12data:  table(count_matrix["HLA-B", ] > 0, cluster)
13X-squared = 1364.4, df = 1, p-value < 2.2e-16

看p值，HLA-B显著性相比HLA-A就非常强，我们就是要挑出和HLA-B类似的，也就是极显著的，设定p值的阈值为0.01

 1HLA_genes <- allGenes[grep('HLA',allGenes)]
 2# 将输出结果保存在向量中
 3HLA_result <- c()
 4for (gene in HLA_genes) {
 5  tmp <- chisq.test(table(count_matrix[gene,]>0, cluster))
 6  if (tmp$p.value<0.01) {
 7    HLA_result[gene] <- gene
 8  }
 9}
10> names(HLA_result)
11 [1] "HLA-F"    "HLA-E"    "HLA-C"    "HLA-B"    "HLA-DRA"  "HLA-DRB5" "HLA-DRB1" "HLA-DQA1" "HLA-DQB1"
12[10] "HLA-DQA2" "HLA-DMB"  "HLA-DMA"  "HLA-DOA"  "HLA-DPA1" "HLA-DPB1"
13> length(HLA_result)
14[1] 15

“

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

”

---