R语言小作业-中级

作者: DrKu | 来源:发表于2020-02-27 14:02 被阅读0次

作业 1

请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
ENSG00000000003
ENSG00000000005
ENSG00000000419
ENSG00000000457
ENSG00000000460
ENSG00000000938
提示：

library(org.Hs.eg.db)
 g2s=toTable(org.Hs.egSYMBOL)
 g2e=toTable(org.Hs.egENSEMBL)

g2se <-  merge(g2s,g2e,by.x="gene_id",by.y="gene_id")
index <- c("ENSG00000000003","ENSG00000000005","ENSG00000000419",
           "ENSG00000000457","ENSG00000000460","ENSG00000000938")

g2se[1:5,1:3]
# gene_id symbol      ensembl_id
#1       1   A1BG ENSG00000121410
#2      10   NAT2 ENSG00000156006
#3     100    ADA ENSG00000196839
#4    1000   CDH2 ENSG00000170558
#5   10000   AKT3 ENSG00000117020
g2se[g2se$ensembl_id %in% index,]
 gene_id   symbol      ensembl_id
# gene_id   symbol      ensembl_id
# 11530    2268      FGR ENSG00000000938
# 21666   55732 C1orf112 ENSG00000000460
# 22360   57147    SCYL3 ENSG00000000457
# 23819   64102     TNMD ENSG00000000005
# 25603    7105   TSPAN6 ENSG00000000003
# 29268    8813     DPM1 ENSG00000000419

作业 2

作业 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
提示：

library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
# 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

index2 <- c("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")
ids[ids$probe_id %in% index2,]

# 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

作业 3

找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图
提示：

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) 
group_list <- pData(sCLLex)
library(hgu95av2.db)

想想如何通过 ggpubr 进行美化。

ids3 <- toTable(hgu95av2SYMBOL)
ids3[ids3$symbol %in% "TP53",]
# probe_id symbol
# 966    1939_at   TP53
# 997  1974_s_at   TP53
# 1420  31618_at   TP53
TP53 <- c(ids3[ids3$symbol %in% "TP53",][,1])
exprSet_TP53 <- exprSet[rownames(exprSet) %in% TP53,]

library(reshape2)
exprSet_TP53_L <- melt(exprSet_TP53)
colnames(exprSet_TP53_L) <- c("probe","sample","value")
exprSet_TP53_L$group <- rep(group_list$Disease,each=nrow(exprSet_TP53))

library(ggplot2)
ggplot(data = exprSet_TP53_L)+geom_boxplot(mapping = aes(x=group,
                                                         y=value))

作业 4

找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets

自己上网操作即可

作业 5

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

自己上网操作即可

作业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，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。

library(GEOquery)
gset <- getGEO("GSE17215",destdir = ".",AnnotGPL = F,
               getGPL = F)
exprSet6 <- exprs(gset[[1]]) 
Pdata6 <- pData(gset[[1]]) 

library(hthgu133a.db)
ids6 <- toTable(hthgu133aSYMBOL)
index6 <- c("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")

probe_idex <- ids6[ids6$symbol %in% index6,][,1]
exprSet_filter <- exprSet6[rownames(exprSet6) %in% probe_idex,]
exprSet_filter= t(scale(t(exprSet_filter)))
pheatmap::pheatmap(exprSet_filter)

image.png

作业7

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

gset7 <- getGEO("GSE24673",destdir = ".",AnnotGPL = F,
               getGPL = F)

##没有分组情况下表达相关性初探
exprSet7 <- exprs(gset7[[1]]) 
Pdata7 <- pData(gset7[[1]]) 
cor_sample <- cor(exprSet7)
pheatmap::pheatmap(cor_sample)


#分组情况下表达相关性
group_list7 <- as.character(Pdata7[,8])
library(stringr)
group_list7=str_split(group_list7," -",simplify = T)[,2]
group_list7[10:11]=c("healthy","healthy")


cor_sample <- cor(exprSet7)
annotation_col = data.frame(
  group = group_list7
)
rownames(annotation_col) =colnames(exprSet7)

pheatmap::pheatmap(cor_sample,annotation_col = annotation_col)

image.png

作业8

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

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/"))
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F)
options()$repos
options()$BioC_mirror

#我本身有这个包，所有就直接加载了
library(hugene10sttranscriptcluster.db)

作业9

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

library(GEOquery)
gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,
               getGPL = F)
###注意变量名字不要和前面几道题的重复了
exprSet9 <- exprs(gset[[1]]) 
Pdata9 <- pData(gset[[1]]) 

library(hugene10sttranscriptcluster.db)
ids9 <- toTable(hugene10sttranscriptclusterSYMBOL)
tail(sort(table(ids9$symbol)))
#可以看到有些symbol对应了多个探针
# RPL41  UBTFL1  CDK11B  UBE2D3    IGKC LRRFIP1 
# 6       6       8       8      10      10
table(sort(table(ids9$symbol)))
#18072个symol和探针是一一对应的
# 1     2     3     4     5     6     8    10 
# 18072   599   132    16     5     6     2     2 
ids9=ids9[ids$symbol != "",]
#顺序一致了
# 1  7896759 LINC01128
# 2  7896761    SAMD11
# 3  7896779    KLHL17
# 4  7896798   PLEKHN1
# 5  7896817     ISG15
ids9=ids9[ids9$probe_id %in% rownames(exprSet9),]

dat9=exprSet9
dim(dat9)
#过滤前有33279个探针
# 33297     6 

dat9=dat9[ids9$probe_id,]
dim(dat9)
#过滤后19827个探针
# 19827     6
dat9[1:5,1:5]
# GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619
# 7896759    8.75126    8.61650    8.81149    8.32067    8.41445
# 7896761    8.39069    8.52617    8.43338    9.17284    9.10216
# 7896779    8.20228    8.30886    8.18518    8.13322    8.06453
# 7896798    8.41004    8.37679    8.27521    8.34524    8.35557
# 7896817    7.72204    7.74572    7.78022    7.72308    7.53797
ids9[1:5,1:2]



####按照平均值
ids9$mean <- apply(dat9,1,mean)
ids9$max <- apply(dat9,1,max)
ids9$sd <- apply(dat9,1,sd)
ids9_mean <- ids9
ids9_max <- ids9
ids_sd <- ids9

ids9_mean=ids9_mean[order(ids9$symbol,ids9$mean,decreasing = T),]
dim(ids9_mean)
#去重复之前
# 19827     5

ids9_mean=ids9_mean[!duplicated(ids9_mean$symbol),]
dat9=dat9[ids9_mean$probe_id,]
dim(dat9)
#去重复之后
# 18834     6
rownames(dat9) <- ids9_mean$symbol


###后面的方法和之前一样，注意重新导入dat9。按照最大值方法
dat9=exprSet9
dim(dat9)
dat9=dat9[ids9$probe_id,]
dim(dat9)

ids9_max=ids9_max[order(ids9_max$symbol,ids9_max$max,decreasing = T),]
dim(ids9_max)
ids9_max=ids9_max[!duplicated(ids9_max$symbol),]
dat9=dat9[ids9_max$probe_id,]
dim(dat9)
rownames(dat9) <- ids9_max$symbol


###按照方差
dat9=exprSet9
dim(dat9)
dat9=dat9[ids9$probe_id,]
dim(dat9)

ids_sd=ids_sd[order(ids_sd$symbol,ids_sd$sd,decreasing = T),]
dim(ids_sd)
ids_sd=ids_sd[!duplicated(ids_sd$symbol),]
dat9=dat9[ids_sd$probe_id,]
dim(dat9)
rownames(dat9) <- ids_sd$symbol

作业10

下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵
This entry was posted in 未分类 by ulwvfje. Bookmark the permalink.

#准备好三个文件：过滤去重后的表达矩阵（上一步的dat9）、design分组文件、contrast.matrix比较文件
exprSet9=dat9
Pdata9 <- pData(gset[[1]]) 
group_list9 <- str_split(Pdata9$title," ",simplify = T)[,4]
library(limma)
design <- model.matrix(~0+factor(group_list9))
colnames(design)=levels(factor(group_list9))
head(design)
# Control Vemurafenib
# 1       1           0
# 2       1           0
# 3       1           0
# 4       0           1
# 5       0           1
# 6       0           1
exprSet=dat9
rownames(design)=colnames(exprSet)
design
# Control Vemurafenib
# GSM1052615       1           0
# GSM1052616       1           0
# GSM1052617       1           0
# GSM1052618       0           1
# GSM1052619       0           1
# GSM1052620       0           1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`factor(group_list9)`
# [1] "contr.treatment"
contrast.matrix<-makeContrasts("Vemurafenib-Control",
                               levels = design)
contrast.matrix ##这个矩阵声明，我们要把 Tumor 组跟 Normal 进行差异分析比较
# Contrasts
# Levels        Vemurafenib-Control
# Control                      -1
# Vemurafenib                   1

  ##step1
  fit <- lmFit(exprSet,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)

  # logFC   AveExpr         t      P.Value    adj.P.Val        B
  # CD36   5.780170  7.370282  79.72556 1.209610e-16 2.278179e-12 26.75898
  # DUSP6 -4.212683  9.106625 -62.43810 1.804535e-15 1.323536e-11 25.01000
  # DCT    5.633027  8.763220  61.56547 2.108212e-15 1.323536e-11 24.89904
  # SPRY2 -3.801663  9.726468 -53.95479 9.056119e-15 4.264073e-11 23.80849
  # MOXD1  3.263063 10.171635  47.08154 4.074111e-14 1.305678e-10 22.59432
  # ETV4  -3.843247  9.667077 -46.99304 4.159535e-14 1.305678e-10 22.57698
exprSet["CD36",]
#验证，说明结果是正确的
# GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 
# 4.54610    4.40210    4.49239   10.25060   10.21480   10.31570

R语言小作业-中级

作业 1

作业 2

作业 3

作业 4

作业 5

作业6

作业7

作业8

作业9

作业10

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读