程序员经验分享-hwhale

生信学习——R语言小作业-中级（附详细答案解读）

2023-10-27

题目目录


写在前面——自从上次做了初级的题目之后，就一直在看这个中级的题目。因为中间有事耽搁了许久，所以间隔了很多天才做完。虽然按照视频和百度磕磕绊绊的把这个题目写完了，但是脑子还是一团浆糊。知道代码是干嘛的，但是不知道为什么要这么做。革命尚未成功，同志仍需努力…

题目原文：http://www.bio-info-trainee.com/3750.html
视频教程：https://www.bilibili.com/video/BV1cs411j75B?p=13
优质答案：https://www.jianshu.com/p/e15ee2cd3174

注意：如果library(…)报错的话，是因为没有安装包，需要install.packegs(…)安装对应的包。

1. 请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)。

ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

# 清空Environment中的变量，也可以直接点小扫帚清空
rm(list = ls())

# 避免把字符串项当成因子
options(stringsAsFactors = FALSE)

a <- read.table("practice/e1.txt")

# 思路：先得到egSYMBOL和egENSEMBL数据框
library(org.Hs.eg.db)
g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)

# 保留a中的v1小数点前面的部分，并将其赋给a的ensembl_id
a$ensembl_id <- unlist(lapply(a$V1, function(x){
  strsplit(x, "[.]")[[1]][1]
}))

# 把a和g2e通过ensembl融合
tmp <- merge(a, g2e, by = "ensembl_id")
# 最后再根据gene_id进行融合
result <- merge(tmp, g2s, by = "gene_id")

在这里插入图片描述

在这里插入图片描述

2. 根据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

rm(list = ls())
options(stringsAsFactors = FALSE)
b <-  read.table("practice/e2.txt")

# 安装hgu133a.db包
# 若出错请查看 https://blog.csdn.net/narutodzx/article/details/119378949
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("hgu133a.db")

library(hgu133a.db)
ids <- toTable(hgu133aSYMBOL)
head(ids)

colnames(b) <- "probe_id"
result <- merge(ids, b, by="probe_id")

在这里插入图片描述

​

3. 找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图。想想如何通过 ggpubr 进行美化。

rm(list = ls())
options(stringsAsFactors = FALSE)

# BiocManager::install("CLL")
# BiocManager::install("hgu95av2.db")

# 隐藏包在加载的时候显示的信息
suppressPackageStartupMessages(library(CLL))
# 加载数据
data(sCLLex)
sCLLex

在这里插入图片描述

# 获得表达矩阵
# 用expr()提取assayData信息
exprSet <- exprs(sCLLex) 
# 获得临床信息
# 用pData()提取phenoData信息
pd <- pData(sCLLex)

library(hgu95av2.db)
ids <- toTable(hgu95av2SYMBOL)
head(ids)

# 在ids中搜索TP53

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

# 绘图
boxplot(exprSet['1939_at',] ~ pd$Disease) #signal
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

我看ggpubr中的ggboxplot函数都是对同一个数据框中的列进行绘图，不知道怎么对两个数据框中的列绘图，挖个坑，学会了再来填…

4. 找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况。

提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets

表达情况：http://www.cbioportal.org/results/plots?cancer_study_list=brca_tcga_pan_can_atlas_2018&Z_SCORE_THRESHOLD=2.0&RPPA_SCORE_THRESHOLD=2.0&profileFilter=mutations%2Cfusion%2Cgistic&case_set_id=brca_tcga_pan_can_atlas_2018_cnaseq&gene_list=BRCA1&geneset_list=%20&tab_index=tab_visualize&Action=Submit&plots_horz_selection=%7B%22dataType%22%3A%22clinical_attribute%22%2C%22selectedDataSourceOption%22%3A%22SUBTYPE%22%7D&plots_vert_selection=%7B%22selectedGeneOption%22%3A672%7D&plots_coloring_selection=%7B%7D
在这里插入图片描述

rm(list = ls())
options(stringsAsFactors = FALSE) 
d <- read.table("practice/plot.txt", sep = "\t", fill = T, header = T)

# 列重新命名
colnames(d) <- c("id", "subtype", "expression", "mut", "cna")

# 绘图
# install.packages("ggstatsplot")
library(ggstatsplot)
ggbetweenstats(data = d, x = subtype, y = expression)

在这里插入图片描述

​

5. 找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存。提示使用：http://www.oncolnc.org/

网站查询：http://www.oncolnc.org/kaplan/?lower=50&upper=50&cancer=BRCA&gene_id=7157&raw=TP53&species=mRNA
在这里插入图片描述

代码验证

rm(list = ls())
options(stringsAsFactors = FALSE) 
e <- read.csv("practice/BRCA_7157_50_50.csv", header = T)

library(ggplot2)
library(survival)
library(survminer) 

# table(e$Status)
e$Status <- ifelse(e$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=e)
sfit
summary(sfit)

# 绘制生存曲线
ggsurvplot(sfit, conf.int=F, pval=TRUE)
# ggsave('survival_TP53_in_BRCA_TCGA.png')

在这里插入图片描述

根据P值可以看出影响不大？
​

6. 下载数据集GSE17215的表达矩阵并且提取下面的基因画热图。

ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T

提示：根据基因名拿到探针ID，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。

rm(list = ls())
options(stringsAsFactors = FALSE)

# 下载GEO数据集遇到网络问题的解决方法：https://www.jianshu.com/p/2686e257b250#comments
# 加载数据
library(GEOquery)
# ?getGEO
# 获取GSE17215_series_matrix.txt.gz
gset <- getGEO("GSE17215", destdir = ".", getGPL = F)
# 保存数据，方便以后使用
save(gset, file = "GSE17215_eSet.Rdata")
# 加载文件，注意文件所在位置
load('GSE17215_eSet.Rdata')  

# 获得表达矩阵 方法同第3题
class(gset)
length(gset)
gset[[1]]
class(gset[[1]])

dat <- exprs(gset[[1]])
dim(dat)

在这里插入图片描述

在这里插入图片描述

获得表达矩阵还有一种方法，见：如何成功getGEO https://www.jianshu.com/p/73580b051fa9?utm_campaign=maleskine…&utm_content=note&utm_medium=seo_notes&utm_source=recommendation

# 处理数据 同第2题
library(hgu133a.db)
ids <- toTable(hgu133aSYMBOL)
# head(ids)
dat <- dat[ids$probe_id,]
head(dat)
dim(dat)

# 对dat中的每行数据取中位数，赋值给ids中新建的median中
ids$median <- apply(dat,1,median)
# 首先对symbol进行排序，然后按照median排序（降序）
ids <- ids[order(ids$symbol,ids$median,decreasing = T),]
# 去除重复的gene ，保留每个基因最大表达量结果
ids <- ids[!duplicated(ids$symbol),]
# 新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
dat <- dat[ids$probe_id,]
# 把ids的symbol这一列中的行名作为dat的行名
rownames(dat) <- ids$symbol
head(dat)
dim(dat)

数据处理前：在这里插入图片描述

数据处理后：在这里插入图片描述

# 绘制热图
ng="ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T"
# 将字符串按空格分隔开，得到一个list，然后取list中的[[1]]，得到分割后的字符串
ng=strsplit(ng,' ')[[1]]
ng

# 将ng中不存在的基因剔除
table(ng %in%  rownames(dat))
ng=ng[ng %in%  rownames(dat)]
table(ng %in%  rownames(dat))

# 根据ng的值在dat中进行匹配
dat=dat[ng,]

# 绘图
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')

在这里插入图片描述

在这里插入图片描述

​

7. 下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息。

rm(list = ls())
options(stringsAsFactors = FALSE)

# 下载并保存数据集
f='GSE24673_eSet.Rdata'
library(GEOquery)
# 设置网络，提高下载成功率
options( 'download.file.method.GEOquery' = 'libcurl' )
if(!file.exists(f)){
  gset <- getGEO('GSE24673', destdir=".",
                 AnnotGPL = F,    
                 getGPL = F)       
  save(gset,file=f)  
}
load('GSE24673_eSet.Rdata')

# 表达矩阵
dat=exprs(gset[[1]])
# dim(dat)
# head(dat)
# 提取phenoData信息
pd=pData(gset[[1]])

# 根据pd数据中的第8列(source_name_ch1)进行分组
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')

# 计算样本的相关性并且绘制热图
M=cor(dat)
pheatmap::pheatmap(M)

# 标记样本分组信息
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

在这里插入图片描述

在这里插入图片描述

​

8. 找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它。

如何查找对应的注释包：https://www.jianshu.com/p/73580b051fa9?utm_campaign=maleskine…&utm_content=note&utm_medium=seo_notes&utm_source=recommendation

# 安装注释包，固定格式
options()$repos
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# 引号中填写对应的注释包名，记得加“.db”
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F)
options()$repos
options()$BioC_mirror

​

9. 下载数据集GSE42872的表达矩阵，并且分别挑选出所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。

rm(list = ls())  
options(stringsAsFactors = F)

# 加载并保存数据
f='GSE42872_eSet.Rdata'
library(GEOquery)
options( 'download.file.method.GEOquery' = 'libcurl' )
if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     
                 getGPL = F)       
  save(gset,file=f)   
}
load('GSE42872_eSet.Rdata')

class(gset)
length(gset)
class(gset[[1]])

a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)

# 根据箱图可以大致看出数据的分布情况
boxplot(dat)
# 对数据的每行取mean,sd,mad，然后按降序排列，再取第一个最大值
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]

在这里插入图片描述

在这里插入图片描述

​

10. 下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵。

rm(list = ls()) 
options(stringsAsFactors = F)

# 同第9题
f='GSE42872_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   
}
load('GSE42872_eSet.Rdata') 

class(gset)
length(gset)
class(gset[[1]])

a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)

boxplot(dat)

# 设置分组信息
group_list=unlist(lapply(pd$title,function(x){
  strsplit(x,' ')[[1]][4]
}))
group_list

# 设计比较矩阵
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(dat)
design

# 自定义比较元素
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix 

# 参考：https://www.jianshu.com/p/d450f35db1cd
# 没看懂......
# step1
fit <- lmFit(dat,design)
# step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE

# step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

生信学习

数据分析

生信分析

R语言

生信学习——R语言小作业-中级（附详细答案解读） 的相关文章

随机推荐

  • string数组转int数组

    string数组类型转换为int数组 方法一 ConvertAll的用法 1 public static int StrToInt string str 2 3 return int Parse str 4 5 6 string arrs

  • 华为CE12808/S9700交换机istack/CSS堆叠主备倒换操作命令步骤

    一 华为CE12808交换机 istack堆叠状态 1 设备型号 交换机一 HUAWEI CE12808 交换机二 HUAWEI CE12808 2 istack堆叠主备倒换操作步骤 2 1 设备当前配置保存并进行备份 2 2 切换所用命令

  • Flink java模拟生成自定义流式数据

    思路如下 定义一个POJO类 注意flink里使用的类必须有一个无参的构造方法 自定义DataSource实现SourceFunction接口 使用ctx collect 传入想要发送的数据就可以了 首先定义一个POJO类 class My

  • 渗透测试工程师的一些面试题3（同样适合刚入门的小白看哦~~~！）

    SQL注入防护 1 使用安全的API 2 对输入的特殊字符进行Escape转义处理 3 使用白名单来规范化输入验证方法 4 对客户端输入进行控制 不允许输入SQL注入相关的特殊字符 5 服务器端在提交数据库进行SQL查询之前 对特殊字符进行

  • TypeError: 'method' object is not subscriptable

    TypeError method object is not subscriptable 此错误一般是函数没加括号导致 下图写错了应该是 case id self sheet cell row 1 value

  • C++ - 重载函数与模板函数(function template)

    参考 CSDN C C 中函数重载的理解 Essential C 一 重载函数 1 1 重载函数的意义 重载函数通常用来在同一个作用域内用同一个函数名命名一组功能相似的函数 这样做减少了函数名的数量 避免了名字空间的污染 对于程序的可读性有

  • 什么是DDoS攻击

    DDoS攻击的定义 DDoS Distributed Denial of Service 分布式拒绝服务 攻击 是指攻击者控制僵尸网络中的大量僵尸主机向攻击目标发送大流量数据 耗尽攻击目标的系统资源 导致其无法正常地响应服务请求 如下图所示

  • 记一次Linux服务器上查杀木马经历

    开篇前言 Linux服务器一直给我们的印象是安全 稳定 可靠 性能卓越 由于一来Linux本身的安全机制 Linux上的病毒 木马较少 二则由于宣称Linux是最安全的操作系统 导致很多人对Linux的安全性有个误解 以为它永远不会感染病毒

  • P1229 遍历问题洛谷

    洛谷P1229 遍历问题 题目描述 我们都很熟悉二叉树的前序 中序 后序遍历 在数据结构中常提出这样的问题 已知一棵二叉树的前序和中序遍历 求它的后序遍历 相应的 已知一棵二叉树的后序遍历和中序遍历序列你也能求出它的前序遍历 然而给定一棵二

  • 计算机网络总结（2）

    四 分类的IP地址 IPv4的地址长度为32bit 标准分类的IP地址是由网络号和主机号组成 用点分十进制表示 IP地址的指派范围 一般不使用的特殊IP地址 五 数据报封装和分用 1 IP数据报格式 IP数据报的格式能够说明IP协议都具有什

  • 启动终端快捷键

    Ctrl Alt t

  • 【SPDK】六、vhost子系统

    vhost子系统在SPDK中属于应用层或叫协议层 为虚拟机提供vhost blk vhost scsi和vhost nvme三种虚拟设备 这里我们以vhost blk为分析对象 来讨论vhost子系统基本原理 vhost子系统初始化 vho

  • chemdraw怎么连接两个结构_科研小站

    化学结构绘图必备软件 ChemDraw了解一下 软件简介 ChemDraw是由CambridgeSoft公司制作的ChemOffice系列软件 是目前国内外最流行 最有应用价值的化学绘图软件 可以快速 精确地绘制化学结构 是各期刊指定的化学

  • C语言---形参所导致的段错误

    前言 今天刷B站 无意之间看到一个宣称90 人都会错的嵌入式面试题 感兴趣就看了一下 卡了十多分钟才想明白 只是一个小知识点 但还是分享一下 题目 include

  • 项目：UDP聊天室

    UDP UDP User Datagram Protocol 是一种无连接 不可靠 面向数据报的传输协议 与TCP相比 UDP更加轻量级 不提供像TCP那样的可靠性和流控制机制 但具备较低的通信延迟和较少的开销 UDP具有以下几个特点 1

  • Mosaic 【HDU - 4819】【二维线段树】

    题目链接 这道题难就只是难在题目难读 题意读懂后就是一道普通的二维线段树更新查找问题 题意 给你一个N N的矩阵 并且已经建立了初始值 然后给你个点以及L 很多人不解其义 其实就是给你个点 然后查的是以 x y 为基础的点 在以左上角 x

  • Python 图像处理入门（识别物体轮廓）

    经过灰度化 高斯滤波 二值化后得到大致轮廓 根据轮廓大小筛选部分干扰元素 尽可能只保留主要物体的轮廓 1 灰度化 图像灰度化是指将图像从彩色转换为灰度图像的过程 这通常是通过将每个像素的 RGB 值转换为单个灰度值来完成的 灰度值通常是使用

  • 剑指offer——连续子数组的最大和

    题目链接 https www nowcoder com practice 459bd355da1549fa8a49e350bf3df484 tpId 13 tqId 11183 rp 1 ru 2Factivity 2Foj qru 2Ft

  • -------Python中ConfigArgParse模块介绍---------

    来源 https pypi org project ConfigArgParse import configargparse p configargparse ArgParser default config files etc app c

  • 生信学习——R语言小作业-中级（附详细答案解读）

    题目目录 1 请根据R包org Hs eg db找到下面ensembl 基因ID 对应的基因名 symbol 2 根据R包hgu133a db找到下面探针对应的基因名 symbol 3 找到R包CLL内置的数据集的表达矩阵里面的TP53基因

热门标签