本笔记来源于B站@生信技能树-jimmy；学习视频链接: 「生信技能树」单细胞数据挖掘
- 下载、探索数据
## 1.1 下载、探索数据
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465
sessionInfo()
a <- read.table("../../rawdata/GSE84465_GBM_All_data.csv.gz")
a[1:4,1:4]
SingleCell_1.JPG
# 行名为symbol ID; 列名为sample，看上去像是两个元素的组合。
summary(a[,1:4])
boxplot(a[,1:4])
head(rownames(a))
tail(rownames(a),10)
# 可以看到原文的counts矩阵来源于htseq这个计数软件，所以有一些不是基因的行需要剔除：
# "no_feature" "ambiguous" "too_low_aQual" "not_aligned" "alignment_not_unique"
tail(a[,1:4],10)
a=a[1:(nrow(a)-5),] #去掉最后不相关的5行
#原始counts数据
# 3,589 cells of 4 human primary GBM samples, accession number GSE84465
# 2,343 cells from tumor cores and 1,246 cells from peripheral regions
# 下载并加载metadata
b <- read.table("../../rawdata/SraRunTable.txt",
sep = ",", header = T)
b[1:4,1:4]
table(b$Patient_ID) # 4 human primary GBM samples
table(b$TISSUE) # tumor cores and peripheral regions
table(b$TISSUE,b$Patient_ID)
SingleCell_2.JPG
- 整理数据
## 1.2 整理数据
# tumor and peripheral 分组信息
head(colnames(a))
head(b$plate_id)
head(b$Well)
# a的矩阵行名（sample）并非GSM编号，而主要是由相应的plate_id与Well组合而成
b.group <- b[,c("plate_id","Well","TISSUE","Patient_ID")]
paste0("X",b.group$plate_id[1],".",b.group$Well[1]) # 验证一下
b.group$sample <- paste0("X",b.group$plate_id,".",b.group$Well)
head(b.group)
identical(colnames(a),b.group$sample)
SingleCell_3.JPG
# 筛选tumor cell
index <- which(b.group$TISSUE=="Tumor")
length(index)
group <- b.group[index,] #筛选的是行
head(group)
a.filt <- a[,index] #筛选的是列。这里尤其需要注意！！！(原来还可以这样操作！)
dim(a.filt)
identical(colnames(a.filt),group$sample)
