单细胞测序下游分析 P 2.1 下载 探索 整理

「生信技能树」单细胞数据挖掘_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1pa4y1s76J?p=4这部分主要是下载GEO数据后，进行数据特征的探索以及整理。

原始代码可参见up主之前的上传github地址。

GEO数据下载地址：

GEO Accession viewerNCBI's Gene Expression Omnibus (GEO) is a public archive and resource for gene expression data.https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465GEO数据下载的meta信息地址:

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA330719https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA330719里面SRA里面有每个sample的信息以及分组、组织来源等信息。

点击select里面download里面total的metadata下载即可。

一、样本表达矩阵及样本信息下载及查看

### 1、下载、探索、整理数据----
## 1.1 下载、探索数据
#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465
sessionInfo()
setwd("")
getwd()
# 读取文件耗时比较长，请耐心等待
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),]

之前都是对a这个数据集进行内容的基本探索，发现a最后5行不是基因ID，因此删除。

> tail(rownames(a),10)
 [1] "ZYG11B"               "ZYX"                  "ZZEF1"               
 [4] "ZZZ3"                 "tAKR"                 "no_feature"          
 [7] "ambiguous"            "too_low_aQual"        "not_aligned"         
[10] "alignment_not_unique"

 这里最后5行不是基因名称。但每个芯片可能并不相同。

---

b文件此时为metadata即样本及文件注释信息。

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

查看不同样本及细胞的信息。

---

二、将metadata中的信息与表达矩阵的id信息相对应起来。

## 1.2 整理数据 
# tumor and peripheral 分组信息
head(colnames(a))
head(b$plate_id)
head(b$Well)
#a矩阵行名（sample）并非为GSM编号，而主要是由相应的plate_id与Well组合而成

a矩阵行名（sample）并非为GSM编号，而主要是由相应的plate_id与Well组合而成

因此可以从b中构建与a相同的样本名称信息的内容，与a进行匹配，进而对a的样本信息进行标注。

按照a的样本名称，使用paste()语句进行b中样本信息的黏贴。

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)

identical()函数来验证样本的列名与sample名一致。

三、筛选肿瘤细胞

# 筛选tumor cell
index <- which(b.group$TISSUE=="Tumor")
length(index)
group <- b.group[index,] #筛选的是行
head(group)

筛选meta信息中肿瘤对应的行的信息。which函数筛选肿瘤所在行，并提取为group。

a.filt <- a[,index] #筛选的是列
dim(a.filt)
identical(colnames(a.filt),group$sample)

sessionInfo()

筛选肿瘤表达矩阵中的列的相关信息。

必须要求meta信息以及肿瘤列的顺序一致才能使用同样的index。

使用identical来验证是否两者顺序一致。