R语言小作业
一、根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
提示：
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
-
Ensembl基因组数据库项目是欧洲生物信息学研究所和Wellcome Trust Sanger研究所之间的一个联合科学项目，该项目于1999年启动，以应对即将完成的人类基因组计划。[2] Ensembl旨在为遗传学家，分子生物学家和其他研究我们自己的物种和其他脊椎动物和模式生物的基因组的研究人员提供集中资源。[3] Ensembl是用于检索基因组信息的几种众所周知的基因组浏览器之一。
-
Ensembl ID主要由五部分组成：
1、ENS：这个开头表示这个是Ensembl id
2、物种：也是几个英文字母，MUS代表小鼠，如果是人则为空
3、类型：E代表exon，FM代表Ensembl protein family，G代表gene，GT代表gene tree，P代表protein，R代表regulatory feature，T代表transcript
一系列数字
5、版本号
其实只要注意前面3个部分就行了
第一步：删除已存在变量和使用命令（stringsAsFactors = FALSE）以防止出错（R often uses a concept of factors to re-encode strings. This can be too early and too aggressive. Sometimes a string is just a string.To avoid problems delay re-encoding of strings by using stringsAsFactors = FALSE when creating data.frames.）
> rm(list=ls())
> options(stringsAsFactors = F)
第二步：导入数据：
e1<-read.table("clipboard",header=T,sep=',') #读取剪切板的内容即其他地方复制后，直接使用该命令调取复制的内容。
或者直接新建.txt文档，将内容复制进去：
image.png
在将数据读进去
> a=read.table('e1.txt')第三步：加载包
> library(org.Hs.eg.db)
载入需要的程辑包：AnnotationDbi
载入需要的程辑包：stats4
载入需要的程辑包：BiocGenerics
载入需要的程辑包：parallel
载入程辑包：‘BiocGenerics’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport,
clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply,
parSapply, parSapplyLB
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, append, as.data.frame, basename, cbind, colMeans,
colnames, colSums, dirname, do.call, duplicated, eval, evalq, Filter,
Find, get, grep, grepl, intersect, is.unsorted, lapply, lengths, Map,
mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
Position, rank, rbind, Reduce, rowMeans, rownames, rowSums, sapply,
setdiff, sort, table, tapply, union, unique, unsplit, which, which.max,
which.min
载入需要的程辑包：Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with 'browseVignettes()'.
To cite Bioconductor, see 'citation("Biobase")', and for packages
'citation("pkgname")'.
载入需要的程辑包：IRanges
载入需要的程辑包：S4Vectors
载入程辑包：‘S4Vectors’
The following object is masked from ‘package:base’:
expand.grid
载入程辑包：‘IRanges’
The following object is masked from ‘package:grDevices’:
windows
了解一下这个包的作用> ?org.Hs.eg.db
image.png
根据末尾提示：
image.png
我们查看一下这个包内有什么内容：
image.png
image.png
并不知道是什么鬼，但是可以找到和题目提示一样的字样。按提示操作，并查看g2s和g2e发现：
> g2s=toTable(org.Hs.egSYMBOL)
> g2e=toTable(org.Hs.egENSEMBL)
image.png
image.png
image.png
发现我们已有的信息ensembl_id,并且得知symbol（对象）这一列表示的是基因名，由此确定答题方向，通过ensembl_id确定gene_id,再通过gene_id确定基因名。
image.png
我们在g2e和我们已知的数据a的ensembl_id不一样，区别在于最后的版本号，我们已有数据有版本号，而得到的g2e没有版本号，所以先将其版本号去掉。
> lapply(a$V1,function(x){strsplit(x,'[.]')})
[[1]]
[[1]][[1]]
[1] "ENSG00000000003" "13"
[[2]]
[[2]][[1]]
[1] "ENSG00000000005" "5"
[[3]]
[[3]][[1]]
[1] "ENSG00000000419" "11"
[[4]]
[[4]][[1]]
[1] "ENSG00000000457" "12"
[[5]]
[[5]][[1]]
[1] "ENSG00000000460" "15"
[[6]]
[[6]][[1]]
[1] "ENSG00000000938" "11"
> lapply(a$V1,function(x){strsplit(x,'[.]')[[1]][1]
+ })
[[1]]
[1] "ENSG00000000003"
[[2]]
[1] "ENSG00000000005"
[[3]]
[1] "ENSG00000000419"
[[4]]
[1] "ENSG00000000457"
[[5]]
[1] "ENSG00000000460"
[[6]]
[1] "ENSG00000000938"
> tmp=merge(a,g2e,by='ensembl_id')
> tmp
ensembl_id V1 gene_id
1 ENSG00000000003 ENSG00000000003.13 7105
2 ENSG00000000005 ENSG00000000005.5 64102
3 ENSG00000000419 ENSG00000000419.11 8813
4 ENSG00000000457 ENSG00000000457.12 57147
5 ENSG00000000460 ENSG00000000460.15 55732
6 ENSG00000000938 ENSG00000000938.11 2268
> daan=merge(tmp,g2s,by='gene_id')
> daan
gene_id ensembl_id V1 symbol
1 2268 ENSG00000000938 ENSG00000000938.11 FGR
2 55732 ENSG00000000460 ENSG00000000460.15 C1orf112
3 57147 ENSG00000000457 ENSG00000000457.12 SCYL3
4 64102 ENSG00000000005 ENSG00000000005.5 TNMD
5 7105 ENSG00000000003 ENSG00000000003.13 TSPAN6
6 8813 ENSG00000000419 ENSG00000000419.11 DPM1
- 在R中可以使用merge()函数去合并数据框，其强大之处在于在两个不同的数据框中标识共同的列或行。
merge函数的用法：
merge(x, y, …)
# S3 method for default
merge(x, y, …)
# S3 method for data.frame
merge(x, y, by = intersect(names(x), names(y)),
by.x = by, by.y = by, all = FALSE, all.x = all, all.y = all,
sort = TRUE, suffixes = c(".x",".y"), no.dups = TRUE,
incomparables = NULL, …)
x,y:用于合并的两个数据框
by,by.x,by.y：指定依据哪些行合并数据框,默认值为相同列名的列.
all,all.x,all.y：指定x和y的行是否应该全在输出文件.
sort： by指定的列是否要排序.
suffixes： 指定除by外相同列名的后缀.
incomparables： 指定by中哪些单元不进行合并.
二、根据R包hgu133a.db找到下面探针对应的基因名(symbol)
1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_at
- 第一步骤和第一问题一样需要新建一个包含题目探针id的txt文件,在此不重复了，直接读取文件了。
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('e2.txt') #读取需要得到基因名的探针id列表
library(hgu133a.db)
#加载包（需要通过学习包或者前辈告诉你这个包含有包含探针id和基因名相符合的列表）
ids=toTable(hgu133aSYMBOL) #取出包含探针id和基因名相符合的列表
colnames(a)='probe_id' #将a的列名改为探针的id
daan=merge(a,ids,by='probe_id') # 将两个表整合在一起，即得出答案
答案为：
> daan
probe_id symbol
1 1053_at RFC2
2 117_at HSPA6
3 121_at PAX8
4 1255_g_at GUCA1A
5 1316_at THRA
6 1320_at PTPN21
7 1405_i_at CCL5
8 1431_at CYP2E1
9 1438_at EPHB3
10 1487_at ESRRA
11 1494_f_at CYP2A6
12 1598_g_at GAS6
13 160020_at MMP14
14 1729_at TRADD
15 177_at PLD1
在最后合并两个表格除了使用merge函数，还可以使用match函数
> match(a$probe_id,ids$probe_id)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
#match函数将a中的probe_id这一列与ids中probe_id这一列匹配的做成索引（即第一个到第十五个都是匹配的）。
> tmp=ids[match(a$probe_id,ids$probe_id),]
# 再使用索引取出匹配的两列的内容（即取出ids的一到十五行）
> tmp==daan # 与答案一致
probe_id symbol
1 TRUE TRUE
2 TRUE TRUE
3 TRUE TRUE
4 TRUE TRUE
5 TRUE TRUE
6 TRUE TRUE
7 TRUE TRUE
8 TRUE TRUE
9 TRUE TRUE
10 TRUE TRUE
11 TRUE TRUE
12 TRUE TRUE
13 TRUE TRUE
14 TRUE TRUE
15 TRUE TRUE
三、找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图
> rm(list = ls())
> options(stringsAsFactors = F)
> options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/") # 修改镜像
> BiocManager::install('CLL') # 下载bioconductor的cll包
> suppressPackageStartupMessages(library(CLL)) #目前不知道什么意思？
> data('sCLLex') # 加载sCLLex这个数据集
> exprSet=exprs(sCLLex) # 通过exprs得到表达矩阵并赋值给exprSet
> pd=pData(sCLLex) # 通过pData获得分组信息并赋值给pd
> options(BioC_mirror = "http://mirrors.ustc.edu.cn/bioc/")
> BiocManager::install('hgu95av2.db') # 下载bioconductor的hgu95av2.db包
> library(hgu95av2.db)
> ids=toTable(hgu95av2SYMBOL)
中间的失误：
> ids=toTable(hgu95av2.db) # 这不是一个对象，其中包含很多东西
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘toTable’ for signature ‘"ChipDb"’
> ids=toTable(hgu95av2SYMBOL) # 这才是正确的。
image.png
通过ids数据找到tp53对应的探针名，便可画图了
> boxplot(exprSet['1939_at',] ~ pd$Disease)
image.png
> boxplot(exprSet['1974_s_at',] ~ pd$Disease)
image.png
> boxplot(exprSet['31618_at',] ~ pd$Disease)
image.png
四、找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets
打开http://www.cbioportal.org/，操作如下：
image.png
image.png
image.png
image.png
下载数据到本地，放到R能够识别的地方（我的路径：此电脑-文档）
> install.packages("ggstatsplot",repos="http://mirror.bjtu.edu.cn/ ")
Installing package into ‘C:/Users/lijian/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
> install.packages("ggstatsplot")
Installing package into ‘C:/Users/lijian/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
also installing the dependencies ‘ps’省略大量依赖包, ‘WRS2’
There is a binary version available but the source version is later:
binary source needs_compilation
zip 2.0.0 2.0.1 TRUE
Binaries will be installed
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/ps_1.3.0.zip'
Content type 'application/zip' length 307358 bytes (300 KB)
downloaded 300 KB
省略大量依赖包下载链接
package ‘ggstatsplot’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\lijian\AppData\Local\Temp\RtmpYblEDs\downloaded_packages
a=read.table('plot.txt',sep='\t',fill = T)
> rm(list = ls())
> options(stringsAsFactors = F)
> a=read.table('plot.txt',sep = '\t',fill = T,header = T)
> colnames(a)=c('id','subtype','expression','mut')
> dat=a
> library(ggstatsplot)
> ggbetweenstats(data = dat, x = subtype, y = expression)
Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
Note: Shapiro-Wilk Normality Test for expression : p-value = < 0.001
Note: Bartlett's test for homogeneity of variances for factor subtype : p-value = < 0.001
得到另一种形式的图片，但是与网页制作的图片是一致的。
image.png
五、找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存。
提示使用：http://www.oncolnc.org/
打开提示网址：
image.png
image.png
image.png
把下载得到的数据放到Rstudio能够识别的地方。
执行一下这些代码：
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
dat=a
library(ggplot2)
library(survival)
library(survminer)
table(dat$Status)
dat$Status=ifelse(dat$Status=='Dead',1,0)
sfit <- survfit(Surv(Days,Status)~Group, data=dat)
sfit
summary(sfit)
ggsurvplot(sfit,conf.int = F,pval = TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.PNG')
画出和网页一致的图（图片还需进一步查资料了解）
image.png
b=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
colnames(b)=c('Patient','subtype','expression','mut')
b$Patient=substring(b$Patient,1,12)
网友评论