搜索
查看: 2761|回复: 6

【R处理芯片数据】--总结

[复制链接]

13

主题

33

帖子

235

积分

版主

Rank: 7Rank: 7Rank: 7

积分
235
发表于 2017-1-16 15:24:03 | 显示全部楼层 |阅读模式
本帖最后由 Andyhigh 于 2017-1-16 15:24 编辑



首先感谢下总版主Jimmy,是他的“R处理芯片数据”教程带我入了生信分析的门,

现对自己的学习总结如下(几个关键点)(自己学习的过程中也总结了不少的帖子。。。)

希望论坛越来越棒,群主的理念能得到推广!!!!!!


1. R,Rstudio的安装;

可能大家这个安装设置比较简单,但如果设置好了,会对后面的下载速度很有影响,是个关键的地方。

注意以下几个地方的协调就好。。。
R中 Rprofile.site 文件中的镜像,R中的镜像设置,Rstudio中的镜像设置,需要相互匹配,才不会报错,且下载速度很快。。。

具体可见链接:http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=159530

2. 包的安装
对于初学者,这个可能比较头痛,出现问题可以考虑以下几个方面:

2.1 代码写错,如bioconductor中的source(“https:xxxx")   需要注意字母和符号(括号、”.“,”“)的中英文,大小写等,最保险就是去bioconductor的官网去复制代码;

2.2 文件太大,网速不给力,这个认定标准为下载的大小是否与实际大小相符。。。这个没办法,只有换一个网速好的地方或者网速好的时间段。。。

2.3 如果还安装不了,可考虑关闭R,Rstudio重新安装;在Rstudio中不能安装,在R中尝试安装等。。。

2.4 借用总版主的一句话,没有包是安装不了的

这方面自己写了个简单的链接,大家可以看下http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=160148


3. R的基本函数,数据的读入输出,数据的处理(重复值的删除,缺失值的删除,数据的合并,排序,转置等)

这一块只有多试了,自己提出些基本的





这个对后续的分组很重要,此外还有很多,如以下命令中的相关函数,都是需要学习的
group_list <- unlist(lapply(metadata2$title, function(x){strsplit(as.character(x)," ")[[1]][1]}))

这里的函数也会涉及数据的合并merge(),  转置等 t()

2. 数据的读入,输出

For example
TT<-read.table("2D3D.csv",header=TRUE, sep=",")
write.table(D,file="3D2.csv",sep=",")

3.数据的处理
#####################################
###数据读入后,可能需要去除重复的数据,方法如下:
a<-read.table("xxx.csv", header=TRUE, sep=",")
index<-duplicated(a$Symbol)
index
a2<-a[!index,]  
###删除重复值

###数据读入后,根据具体需要,确定是否删除NA值,方法如下:
a3<-a2[na.omit(a2$Symbol),]
write.table(a3,file="single cell ding.csv",sep = ",",quote = F)
a4<-na.omit(a2)
write.table(a4,file="a4.csv",sep = ",",quote = F)
###删除缺失值
####################################################


4. R处理芯片数据教程的应用
学习的目的主要是应用,于是我们需要利用自己学到的知识,分析与自己相关的课题。。。

由于自己这边涉及到造血干细胞的工作,今天上午在GEO 数据库中检索了下,进行了一个两个分组的学习,其中是针对 “5行代码的教程“,根据自己的需要做了一些修改。基本上囊括了目前市面上芯片数据分析的大部分图。。。

具体代码如下:

[AppleScript] 纯文本查看 复制代码
rm(list=ls())
setwd("C:/users/lixia/Desktop/learn/GSE35008")

source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("GEOquery")
biocLite("hugene10sttranscriptcluster.db") 
###出现包安装不上,可用goole搜索这个包的名字,复制bioconductor中的代码,可能原因是输入中的代码不对。

library(humanid)
library(GEOquery)
library(hugene10sttranscriptcluster.db)
library(limma)
eSet <- GEOquery::getGEO('GSE35008', destdir='./',getGPL = F)

exprSet <- get_symbol_exprSet(eSet[[1]],'hugene10sttranscriptcluster.db') 
## GPLhugene10sttranscriptcluster =annotation
metadata <- pData(eSet[[1]])

选取前10组进行分析
lixia<-exprSet[,1:10]
metadata2 <- metadata[1:10,]

group_list <- unlist(lapply(metadata2$title, function(x){
  strsplit(as.character(x)," ")[[1]][1]
})) 


DEG <- do_DEG_2groups(prefix='GSE35008',
                      exprSet= lixia,
                      group_list=group_list,
                      method='limma')

write.table(DEG,file="long short.xls", sep="\t",quote=F) ###永久存入文件到指定目录;

DEG2 <- DEG[order(DEG$adj.P.Val),]  ###排序
write.table(DEG2,file="long short2.xls", sep="\t",quote=F)
sum(DEG$P.Value < 0.05, na.rm=TRUE) ###统计符合条件的gene个数

lx<-summary(DEG)  ###统计中位数
lx


DEG$symbol <- rownames(DEG)  #作用! 无此命立,则出现“no lines available in input”错误

format_DEG(DEG,prefix='GSE35008',GOstats = T)

QCexpressionMatrix(lixia,group_list,prefix='GSE35008')

createGSEAinput(prefix='GSE35008',lixia,group_list)


5. 数据或问题的可视化

可视化,自己认为就是将数据,用图表展示,

自己也写了些帖子:

R画聚类图与PCA图 http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=160301

R画网络图  http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=160421

R画热图 http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=159345


此外,附上一个ID转换的代码,感觉还是比较有用的:
代码为
[AppleScript] 纯文本查看 复制代码
分析数据时,ID转换在所难免,本章节将介绍如何将NCBI 中的Entrez Gene ID转换为Symbol 和 Esembel ID,利用R中的BiomaRt包,具体代码如下:
 
source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)
biocLite("biomaRt")
library(biomaRt)
ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")
entrzID=c("672","1")
getBM(attributes=c("entrezgene","hgnc_symbol","ensembl_gene_id"), 
      filters = "entrezgene", values =entrzID, mart=ensembl )
genesymbol<-getBM(attributes=c("entrezgene","hgnc_symbol",
                               "ensembl_gene_id"), filters = "entrezgene", 
                  values =entrzID, mart=ensembl )
write.table(genesymbol, file="symbol.xls", sep="\t",quote=F)

tmp<-read.table("geneid.txt")  ###需要转换的geneid 集
tmp2<-getBM(attributes=c("entrezgene","hgnc_symbol",
                               "ensembl_gene_id"), filters = "entrezgene", 
                  values =tmp, mart=ensembl )
write.table(tmp2, file="2d3d.xls", sep="\t",quote=F)


http://scu.zju.edu.cn/redir.php?catalog_id=58400&object_id=162298


6. 数据的解读与意义

个人认为,高通量测序的数据,只是为研究者们提供一个思路,各种模型、算法提出的预测都需要验证!!!


针对高通量测序的具体数据,自己认为不能太看重一两个基因的表达,尤其是当下的单细胞数据分析结果,二是应该将重点关注在一个gene set或者moudle上来。针对生信人员和研究者,应该都懂得相关处理,或者包,或者算法的原理,看看自己的数据是否符合要求!!!

数据具体的解读,分析方法和意义 还需要多读文章了,这个没办法!!!


最后上传几个图。。


















本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:fastq数据perl脚本过滤
下一篇:生信编程直播第六题:下载最新版的KEGG信息,并且解析好
回复

使用道具 举报

0

主题

1

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2017-1-16 23:15:47 | 显示全部楼层
嗯,不错!cytoscape教程如果有就更完美了
回复 支持 反对

使用道具 举报

5

主题

18

帖子

142

积分

注册会员

Rank: 2

积分
142
发表于 2017-1-17 11:14:36 | 显示全部楼层
赞!R与Rstudio推荐默认安装,不然要自己设置.包的安装方式也可以参考群主的博客4种R包安装方式http://www.bio-info-trainee.com/1565.html
回复 支持 反对

使用道具 举报

3

主题

43

帖子

202

积分

中级会员

Rank: 3Rank: 3

积分
202
发表于 2017-3-23 18:21:53 | 显示全部楼层
赞一个,学习了
回复 支持 反对

使用道具 举报

2

主题

17

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-4 20:56:48 | 显示全部楼层
赞一个  学习了
回复 支持 反对

使用道具 举报

1

主题

3

帖子

42

积分

新手上路

Rank: 1

积分
42
发表于 2017-7-5 23:25:31 | 显示全部楼层

赞一个  学习了
回复 支持 反对

使用道具 举报

0

主题

3

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2017-12-7 19:14:37 | 显示全部楼层
赞一个。
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-12-18 18:48 , Processed in 0.136154 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.