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

Andyhigh

首先感谢下总版主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的基本函数，数据的读入输出，数据的处理（重复值的删除，缺失值的删除，数据的合并，排序，转置等）

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

1.

R语言中apply系列函数学习

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

这个对后续的分组很重要，此外还有很多，如以下命令中的相关函数，都是需要学习的
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上来。针对生信人员和研究者，应该都懂得相关处理，或者包，或者算法的原理，看看自己的数据是否符合要求！！！

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


最后上传几个图。。

qym059

嗯，不错！cytoscape教程如果有就更完美了

yjl

赞！R与Rstudio推荐默认安装，不然要自己设置.包的安装方式也可以参考群主的博客4种R包安装方式http://www.bio-info-trainee.com/1565.html

galen.cheng

赞一个,学习了

miniaowai

赞一个  学习了

lampirl

赞一个  学习了

bioant_2017

赞一个。