生信技能树单细胞数据挖掘笔记(1)

前言

最近有幸加入生信技能树的单细胞数据挖掘尝鲜群,借此机会给大家分享一下

生信技能树单细胞数据挖掘笔记(1)

生信技能树单细胞数据挖掘笔记(2)
生信技能树单细胞数据挖掘笔记(3)：降维与聚类
生信技能树单细胞数据挖掘笔记(4)：其他分析（周期判断、double诊断、细胞类型注释）
生信技能树单细胞数据挖掘笔记(5)：轨迹分析

课程大纲

入门基础知识

所需R包

getOption("BioC_mirror")
getOption("CRAN")
#CRAN基础包
options(CRAN="https://mirrors.ustc.edu.cn/CRAN/")## 设置镜像
cran_packages <- c('tidyverse',
                   'ggplot2'
                   ) 
for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)

#Bio分析包
Biocductor_packages <- c("Seurat",
                        "scran",
                        "scater",
                        "monocle",
                        "DropletUtils",
                        "SingleR"
                         )
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
#options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
# use BiocManager to install
for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#最后检查下成功与否
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}

### GEO 
# 
# GEO Platform (GPL)
# GEO Sample (GSM)
# GEO Series (GSE)
# GEO Dataset (GDS)

#一篇文章可以有一个或者多个GSE数据集，一个GSE里面可以有一个或者多个GSM样本。
#多个研究的GSM样本可以根据研究目的整合为一个GDS，不过GDS本身用的很少。
#而每个数据集都有着自己对应的芯片平台，就是GPL。(芯片名与基因名ID转换)
#https://blog.csdn.net/weixin_43569478/article/details/108079337
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=

BiocManager::install("GEOquery")
library(GEOquery)
gse1009 <- getGEO('GSE1009', destdir=".")
class(gse1009)
length(gse1009)
a <- gse1009[[1]]
class(gse1009[1])
a
b <- exprs(a)
c <- pData(a)
a$platform_id

下载、探索、数据整理

本次数据采用《Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients》中的单细胞数据进行分析（GSE号：84465）

### 1、下载、探索、整理数据----
## 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]
#行名为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),]

#原始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
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)

整理数据

可以发现两个数据不对应，a矩阵行名（sample）并非为GSM编号，而主要是由相应的plate_id与Well组合而成

## 1.2 整理数据 
# tumor and peripheral 分组信息
head(colnames(a))
[1] "X1001000173.G8" "X1001000173.D4" "X1001000173.B4" "X1001000173.A2" "X1001000173.E2" "X1001000173.F6"
head(b$plate_id)
[1] 1001000173 1001000173 1001000173 1001000173 1001000173 1001000173
head(b$Well)
#a矩阵行名（sample）并非为GSM编号，而主要是由相应的plate_id与Well组合而成

b.group <- b[,c("plate_id","Well","TISSUE","Patient_ID")]
b.group$sample <- paste0("X",b.group$plate_id,".",b.group$Well)
head(b.group)
    plate_id Well TISSUE Patient_ID         sample
1 1001000173   G8  Tumor      BT_S2 X1001000173.G8
2 1001000173   D4  Tumor      BT_S2 X1001000173.D4
3 1001000173   B4  Tumor      BT_S2 X1001000173.B4
4 1001000173   A2  Tumor      BT_S2 X1001000173.A2
5 1001000173   E2  Tumor      BT_S2 X1001000173.E2
6 1001000173   F6  Tumor      BT_S2 X1001000173.F6
identical(colnames(a),b.group$sample)

# 筛选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)

经过筛选之后得到所有的肿瘤细胞的表达矩阵

---

转载请注明：周小钊的博客>>>生信技能树单细胞数据挖掘笔记(1)