1.1 下载、探索数据
对照视频教材：
代码资源可在腾讯微云网盘上寻找：文件分享
本例从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" # 需注意在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的确认。
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/210041.html原文链接：https://javaforall.net
