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

【生信技能树2020-10-31】单细胞数据挖掘学习笔记-1.1 下载、探索数据生信技能树-单细胞数据挖掘学习笔记2.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的确认。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请联系我们举报,一经查实,本站将立刻删除。

发布者:全栈程序员-站长,转载请注明出处:https://javaforall.net/210041.html原文链接:https://javaforall.net

(0)
全栈程序员-站长的头像全栈程序员-站长


相关推荐

  • 运行时错误10048,地址已在使用_无法指定被请求的地址

    运行时错误10048,地址已在使用_无法指定被请求的地址OSError:[WinError10049]是由于ip地址为空造成的设置端口验证的一行manager=QueueManager(address=(”,5000),authkey=b’abc’)中的地址添加127.0.0.1问题解决。

    2022年9月30日
    0
  • kafka主要用来做什么_kafka概念

    kafka主要用来做什么_kafka概念Kafka最初由LinkedIn公司开发的,并于2010年贡献给了Apache基金会,之后成为Apache顶级项目。目前Kafka已经定位为一个分布式流式处理平台,它以高吞吐、可持久化、可水平扩展、支持流数据处理等多种特性而被广泛使用。目前越来越多的开源分布式处理系统如Cloudera、Storm、Spark、Flink等都支持与Kafka集成。Kafka之所以受到越来越多的青睐,与它所“扮演”的三大角色是分不开的:消息系统:Kafka和传统的消息系统(也称作消息中

    2022年9月25日
    0
  • 【csma/ca协议和csma/cd协议的matlab仿真详解】

    【csma/ca协议和csma/cd协议的matlab仿真详解】首先你的熟悉csma/ca协议,csma/cd协议;csma/ca协议:点击打开链接csma/cd协议:点击打开链接这个课题有几个难点部分:1.就是需要考虑各种情况,并对每一种情况都必须做出相应的处理。2.怎么展示,怎么简单、直观而有效的展示你的程序正确性。鉴于本程序我采用动态图形形式展示csma/ca协议的运行过程。以下是我程序的运行结果的部分展示:1…

    2025年7月25日
    0
  • 多目标进化算法详述-MOEA/D与NSGA2优劣比较

    多目标进化算法详述-MOEA/D与NSGA2优劣比较多目标进化算法系列1.多目标进化算法(MOEA)概述2.多目标优化-测试问题及其Pareto前沿3.多目标进化算法详述-MOEA/D与NSGA2优劣比较4.多目标进化算法-约束问题的处理方法NSGA-II由KalyanmoyDeb等人于2002年在文章”AFastandElitistMultiobjectiveGeneticAlgorithm:…

    2022年5月19日
    160
  • 高级SQL查询-(聚合查询,分组查询,联合查询)[通俗易懂]

    高级SQL查询-(聚合查询,分组查询,联合查询)[通俗易懂]高级SQL查询-(聚合查询,分组查询,联合查询)

    2022年6月1日
    32
  • journalctl基本介绍

    journalctl基本介绍journalctl基础用法1、查看所有日志(默认显示本次启动的所有日志)[root@localhost~]#journalctl查看本次启动的所有日志也可以使用[root@localhost~]#journalctl-b2、查看内核日志[root@localhost~]#journalctl-k3、查看指定时间的日志通过–since和–until选项,可以过滤任意时间限制,显示指定条件之前、之后或之间的日志[root@localhost~]#jour

    2022年5月10日
    69

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号