劉小澤寫于19.10.18
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索10X Genomics技術(shù)相關(guān)的分析
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二單元第10講：marker基因的表達(dá)量可視化
前言
這次的任務(wù)是模仿原文的：
腫瘤組織的差異分析(Fig. 4a蔽氨、c)
原文的圖
下載數(shù)據(jù)
之前的分析都是基于第一個(gè)病人的PBMC多柑，這次將基于這位病人的tumor：GSE117988_raw.expMatrix_Tumor.csv.gz
start_time <- Sys.time()
raw_dataTumor <- read.csv('./GSE117988_raw.expMatrix_Tumor.csv.gz', header = TRUE, row.names = 1)
end_time <- Sys.time()
end_time - start_time
# Time difference of 47.00589 secs
dim(raw_dataTumor) # 21861基因，7431細(xì)胞- already filtered
注意：之前PBMC包含了四個(gè)時(shí)間點(diǎn)的樣本,這里的Tumor包含了2個(gè)樣本（治療之前Pre和復(fù)發(fā)AR）
常規(guī)流程
step1: 歸一化
dataTumor <- log2(1 + sweep(raw_dataTumor, 2, median(colSums(raw_dataTumor))/colSums(raw_dataTumor), '*'))
> head(colnames(dataTumor))
[1] "AAACCTGAGGATGTAT.1" "AAACCTGCAGCGATCC.1" "AAACCTGGTACGAAAT.1" "AAACGGGAGCTGGAAC.1" "AAACGGGAGGAGTTGC.1"
[6] "AAACGGGAGTTTAGGA.1"
step2: 自定義劃分細(xì)胞類型
cellTypes <- sapply(colnames(dataTumor), function(x) ExtractField(x, 2, '[.]'))
cellTypes <-ifelse(cellTypes == '1', 'Tumor_Before', 'Tumor_AcquiredResistance')
> table(cellTypes)
cellTypes
Tumor_AcquiredResistance Tumor_Before
5188 2243
step3: 表達(dá)矩陣質(zhì)控
# 第一點(diǎn)：基因在多少細(xì)胞表達(dá)
> fivenum(apply(dataTumor,1,function(x) sum(x>0) ))
VP2 GPRIN2 EML3 ZNF140 RPLP1
1 8 103 566 7431
# 第二點(diǎn)：細(xì)胞中有多少表達(dá)的基因
> fivenum(apply(dataTumor,2,function(x) sum(x>0) ))
GGAACTTAGGAATCGC.1 TACGGTACAAGCCGCT.2 CCTACCAAGCGTGAGT.1 TGAGCCGAGACTAGAT.2 GATCGTAGTCATATGC.2
192.0 1059.0 1380.0 1971.5 5888.0
看到大部分基因在500多個(gè)細(xì)胞表達(dá)掏父，細(xì)胞平均能表達(dá)1000個(gè)基因以上
step4: 創(chuàng)建Seurat對(duì)象
tumor <- CreateSeuratObject(dataTumor,
min.cells = 1, min.features = 0, project = '10x_Tumor')
> tumor
An object of class seurat in project 10x_Tumor
21861 genes across 7431 samples.
step5: 添加metadata (nUMI 和 細(xì)胞類型)
tumor <- AddMetaData(object = tumor, metadata = apply(raw_dataTumor, 2, sum), col.name = 'nUMI_raw')
tumor <- AddMetaData(object = tumor, metadata = cellTypes, col.name = 'cellTypes')
step6: 聚類標(biāo)準(zhǔn)流程
start_time <- Sys.time()
tumor <- ScaleData(object = tumor, vars.to.regress = c('nUMI_raw'), model.use = 'linear', use.umi = FALSE)
tumor <- FindVariableGenes(object = tumor, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
tumor <- RunPCA(object = tumor, pc.genes = tumor@var.genes)
tumor <- RunTSNE(object = tumor, dims.use = 1:10, perplexity = 25)
end_time <- Sys.time()
end_time - start_time
# Time difference of 3.324982 mins
TSNEPlot(tumor, group.by = 'cellTypes', colors.use = c('#EF8A62', '#67A9CF'))
save(tumor,file = 'patient1.Tumor.V2.output.Rdata') # 3.6Gb 大小
分成兩群
接著進(jìn)行基因可視化
rm(list = ls())
options(warn=-1)
start_time <- Sys.time()
load('patient1.Tumor.V2.output.Rdata')
end_time <- Sys.time()
end_time - start_time
# Time difference of 19.83097 secs
取出log歸一化后的表達(dá)矩陣
count_matrix=tumor@data
> count_matrix[1:4,1:4]
AAACCTGAGGATGTAT.1 AAACCTGCAGCGATCC.1 AAACCTGGTACGAAAT.1 AAACGGGAGCTGGAAC.1
VP2 0.0000000 0 0 0
largeTAntigen 0.9670525 0 0 0
smallTAntigen 0.0000000 0 0 0
RP11-34P13.7 0.0000000 0 0 0
取出細(xì)胞分群信息
cluster=tumor@meta.data$cellTypes
> table(cluster)
cluster
Tumor_AcquiredResistance Tumor_Before
5188 2243
提取基因信息
文章主要探索了治療前和復(fù)發(fā)后的HLA-A和HLA-B的變化笋轨，于是我們先看看有沒有這兩個(gè)基因
allGenes = row.names(tumor@raw.data)
> allGenes[grep('HLA',allGenes)]
[1] "HHLA3" "HLA-F" "HLA-G" "HLA-A" "HLA-E" "HLA-C" "HLA-B" "HLA-DRA" "HLA-DRB5"
[10] "HLA-DRB1" "HLA-DQA1" "HLA-DQB1" "HLA-DQA2" "HLA-DQB2" "HLA-DOB" "HLA-DMB" "HLA-DMA" "HLA-DOA"
[19] "HLA-DPA1" "HLA-DPB1"
對(duì)HLA-A操作
FeaturePlot(object = tumor,
features.plot ='HLA-A',
cols.use = c("grey", "blue"),
reduction.use = "tsne")
HLA-A表達(dá)
看HLA-A表達(dá)量
> table(count_matrix['HLA-A',]>0, cluster)
cluster
Tumor_AcquiredResistance Tumor_Before
FALSE 2282 1057
TRUE 2906 1186
對(duì)HLA-B操作
FeaturePlot(object = tumor,
features.plot ='HLA-B',
cols.use = c("grey", "blue"),
reduction.use = "tsne")
可以看到治療前HLA-A基因有1186個(gè)表達(dá)，1057個(gè)不表達(dá)赊淑；復(fù)發(fā)后這個(gè)基因表達(dá)和不表達(dá)的數(shù)量也相近
HLA-B表達(dá)
> table(count_matrix['HLA-B',]>0, cluster)
cluster
Tumor_AcquiredResistance Tumor_Before
FALSE 4794 1258
TRUE 394 985
對(duì)HLA-B來講爵政，不管是治療前還是復(fù)發(fā)后陶缺，它的表達(dá)和不表達(dá)差異就很明顯钾挟。另外從治療前到復(fù)發(fā)饱岸，這個(gè)基因的表達(dá)數(shù)量的變化更顯著
小問題：HLA基因這么多掺出，我們?cè)趺凑业饺康木哂邢嗨颇Ｊ降幕颍?/h4>
所謂相似表達(dá)模式苫费，就是像HLA-B基因一樣汤锨，在一個(gè)群表達(dá)很多百框，另一個(gè)群表達(dá)很少闲礼，通過卡方檢驗(yàn)就能看出區(qū)別
> chisq.test(table(count_matrix['HLA-A',]>0, cluster))
Pearson's Chi-squared test with Yates' continuity correction
data: table(count_matrix["HLA-A", ] > 0, cluster)
X-squared = 6.1069, df = 1, p-value = 0.01347
> chisq.test(table(count_matrix['HLA-B',]>0, cluster))
Pearson's Chi-squared test with Yates' continuity correction
data: table(count_matrix["HLA-B", ] > 0, cluster)
X-squared = 1364.4, df = 1, p-value < 2.2e-16
看p值，HLA-B顯著性相比HLA-A就非常強(qiáng)柬泽，我們就是要挑出和HLA-B類似的慎菲，也就是極顯著的聂抢，設(shè)定p值的閾值為0.01
HLA_genes <- allGenes[grep('HLA',allGenes)]
# 將輸出結(jié)果保存在向量中
HLA_result <- c()
for (gene in HLA_genes) {
tmp <- chisq.test(table(count_matrix[gene,]>0, cluster))
if (tmp$p.value<0.01) {
HLA_result[gene] <- gene
}
}
> names(HLA_result)
[1] "HLA-F" "HLA-E" "HLA-C" "HLA-B" "HLA-DRA" "HLA-DRB5" "HLA-DRB1" "HLA-DQA1" "HLA-DQB1"
[10] "HLA-DQA2" "HLA-DMB" "HLA-DMA" "HLA-DOA" "HLA-DPA1" "HLA-DPB1"
> length(HLA_result)
[1] 15
