【生信技能树2020-10-31】单细胞数据挖掘学习笔记-1.1 下载、探索数据

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全家桶1年46，售后保障稳定

1.1 下载、探索数据

对照视频教材：

「生信技能树」单细胞数据挖掘_哔哩哔哩_bilibili

代码资源可在腾讯微云网盘上寻找：文件分享
 

本例从P4 2.1开始进行笔记记录及整理。本文档主要是操作中的细节补充，相应的代码及思路在原始的代码中已经注明，因此无需特殊补充。

文件下载 例如GSE84465，可于GEO网页上寻找

下载csv.gz文件，download的http网址即可

先基于setup0安装所有的包。

注意工作文档和下载下的文件最好在一个文件夹下，方便后续调取取用。

看1-4列及1-4行初步看一下基因表达矩阵。

列名此时是sample，是根据sample情况决定的。行名是symbol ID。

后续是对矩阵进行基本描述

### 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"   

Jetbrains全家桶1年46，售后保障稳定

# 需注意在GEO网页原文中寻找counts数的来源，用什么样的测序技术及计数软件和技术，本例中使用htseq来计数，会存在非基因的行（最后5行）。因此可以去掉最后5行。

# 可以看到原文的counts矩阵来源于htseq这个计数软件，所以有一些不是基因的行需要剔除：
#  "no_feature"           "ambiguous"            "too_low_aQual"        "not_aligned"          "alignment_not_unique"
tail(a[,1:4],10)

a=a[1:(nrow(a)-5),]

【关于基因注释的基本概念可参见网络材料教程，基本就是利用各类测序数据的基因序列应当相似的原理，注释到相应的基因，进而基于GO库来识别该基因的位置、分子功能与参与的通路】

【关于htseq，GFF/GTF的基本概念可参见：mRNA-seq学习（三）：htseq-count计数 – 简书、GTF/GFF文件的差异及其相互转换 – 简书】

【在CSDN中查询相应概念也是很好的方法。】

后续需明确细胞来源及类型，哪些是肿瘤或者正常组织细胞。该内容可在GEO数据库的样本下面SRA RUN Selector处获取。打开网页后即可看到patient_ID以及tissue便是组织细胞来源。此时plate.ID+Well便是每一个sample（列）样本的组合名称。
（具体可打开网址：https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA330719&o=acc_s%3Aa）

此时提供meta数据予以标注该类信息。将之可以下载下来。

检查原始的counts数据

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

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

整理数据，首先先看转录组矩阵的名称和meta数据的矩阵名称是否一致。

b[]可以将其中需要的信息提取出来。

paste0（）将这些名名称进行黏贴，检查是不是黏贴后与表达矩阵的格式一致。

将所有sample黏贴，然后制造一个新的sample列，对其进行检查。

identical()检查是否所有元素均为一致：a的每个colnames（每个样本的名字）与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)

sessionInfo()

which函数筛选metadata中的是tumor的行的信息

【此时用mode查看b.group可以发现b.group是list，class()查看可以发现是dataframe，因此which函数所返回的是一个索引值即为行的索引数】

length来确认此时的index数量是对的。

从b.group中来筛选相应行数的值，并予以head查看

因为前文中每个colnames的顺序与b.group的的sample顺序是一致的，因此index在表达矩阵的列顺序与colnames一致。

dim用以查看该矩阵的行列数量。dim可以查看该数据集的维度的问题。

最后进行identical的确认。