生信人的20个R语言习题及其答案-土豆学习笔记

作者: 土豆学生信 | 来源:发表于2019-01-28 18:26 被阅读92次

生信人的20个R语言习题 http://www.bio-info-trainee.com/3409.html
生信人的20个R语言习题答案 http://www.bio-info-trainee.com/3415.html
参考答案：
https://cldiao.github.io/2018/09/27/R%E8%AF%AD%E8%A8%8020%E9%A2%98/
http://rvdsd.top/2018/08/30/R/20%E4%B8%AAR%E8%AF%AD%E8%A8%80%E4%B9%A0%E9%A2%98/
https://www.jianshu.com/p/788010093c90

1. 安装一些R包：

数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
不同领域的R包使用频率不一样，在生物信息学领域，尤其需要掌握bioconductor系列包。

if(F){
source("http://bioconductor.org/biocLite.R")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")#修改镜像，安装会加速
BiocManager::install("clusterProfiler")
BiocManager::install("ComplexHeatmap")
BiocManager::install("maftools")
BiocManager::install("ggplot2")
BiocManager::install("jmzeng1314/biotrainee")
}
#或者如下：
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
BiocManager::install(c('ALL','CLL','pasilla','clusterProfiler')) 
BiocManager::install(c('airway','DESeq2','edgeR','limma'))
install.packages("reshape2", "ggplot2")

2.了解ExpressionSet对象

比如CLL包里面就有data(sCLLex) ，找到它包含的元素，提取其表达矩阵(使用exprs函数)，查看其大小

参考：http://www.bio-info-trainee.com/bioconductor_China/software/limma.html

参考：https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)
##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
dim(exprSet)
# [1] 12625    22
exprSet[1:5,1:5]
#            CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904
# 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170
# 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113
# 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243
# 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932

3.了解 str,head,help函数，作用于第二步提取到的表达矩阵

str(exprSet)
# str: Compactly display the internal structure of an R object, a diagnostic function and an alternative to summary (and to some extent, dput).
head(exprSet)

4. 安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后显示的那些变量

hgu95av2.db是一个注释包，它为hgu95av2平台的芯片提供注释，这个包中有很多注释文件，如下所示：

BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ls("package:hgu95av2.db")
#[1] "hgu95av2"              "hgu95av2.db"           "hgu95av2_dbconn"       "hgu95av2_dbfile"       "hgu95av2_dbInfo"       "hgu95av2_dbschema"    
 [7] "hgu95av2ACCNUM"        "hgu95av2ALIAS2PROBE"   "hgu95av2CHR"           "hgu95av2CHRLENGTHS"    "hgu95av2CHRLOC"        "hgu95av2CHRLOCEND"    
[13] "hgu95av2ENSEMBL"       "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"      "hgu95av2ENZYME"        "hgu95av2ENZYME2PROBE"  "hgu95av2GENENAME"     
[19] "hgu95av2GO"            "hgu95av2GO2ALLPROBES"  "hgu95av2GO2PROBE"      "hgu95av2MAP"           "hgu95av2MAPCOUNTS"     "hgu95av2OMIM"         
[25] "hgu95av2ORGANISM"      "hgu95av2ORGPKG"        "hgu95av2PATH"          "hgu95av2PATH2PROBE"    "hgu95av2PFAM"          "hgu95av2PMID"         
[31] "hgu95av2PMID2PROBE"    "hgu95av2PROSITE"       "hgu95av2REFSEQ"        "hgu95av2SYMBOL"        "hgu95av2UNIGENE"       "hgu95av2UNIPROT"

5. 理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID

?hgu95av2SYMBOL
?toTable
summary(hgu95av2SYMBOL)
#SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
|
| Lkeyname: probe_id (Ltablename: probes)
|    Lkeys: "1000_at", "1001_at", ... (total=12625/mapped=11460)
|
| Rkeyname: symbol (Rtablename: gene_info)
|    Rkeys: "A1BG", "A2M", ... (total=61050/mapped=8585)
|
| direction: L --> R

ids <- toTable(hgu95av2SYMBOL)
View(ids)
library(dplyr)
# 方法1：
filter(ids, symbol=="TP53") #用dplyr包的筛选功能，找到 TP53 基因对应的探针ID
#     probe_id symbol
#1   1939_at   TP53
#2 1974_s_at   TP53
#3  31618_at   TP53

#方法2：
ids[grep("^TP53$", ids$symbol),]
#       probe_id symbol
# 966    1939_at   TP53
# 997  1974_s_at   TP53
# 1420  31618_at   TP53
#方法1，2虽然结果相同，但是定义的对象是不同的

hug95av2SYMBOL是一个R对象，它提供的是芯片生产厂家与基因缩写之间的映射信息。这个映射的信息主要依据Entrez Gene数据库。现在我们通过mappedkeys()这个函数，得到映射到基因上的探针信息。

6.理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？

length(unique(ids$symbol))
# [1] 8585
tail(sort(table(ids$symbol)))
#YME1L1  GAPDH INPP4A    MYB PTGER3  STAT1 
#     7      8      8      8      8      8 
table(sort(table(ids$symbol)))
#  1    2    3    4    5    6    7    8 
# 6555 1428  451  102   22   16    6    5

image.png

不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量。

7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

提示：有1165个探针是没有对应基因名字的。

%in% 逻辑判断

用法 a %in% table
a值是否包含于table中，为真输出TURE，否者输出FALSE

table(rownames(exprSet)) %in% ids$probe_id
# %in% is a more intuitive interface as a binary operator, which returns a logical vector indicating if there is a match or not for its left operand.
n_exprSet <- exprSet[!(rownames(exprSet) %in% ids$probe_id),]
dim(n_exprSet)
# [1] 1165   22
View(n_exprSet)
# These probes are not in the package.

8.过滤表达矩阵，删除那1165个没有对应基因名字的探针。

方法1：%in% 逻辑判断

exprSet <- exprSet[rownames(exprSet) %in% ids$probe_id, ]
dim(exprSet)
[1] 11460    22
View(exprSet)
# These probes are in the package.

方法2 mappedkeys() 映射关系

 length(hgu95av2SYMBOL)
[1] 12625
probe_map <- hgu95av2SYMBOL
length(probe_map)
[1] 12625
#全部的探针数目
# [1] 12625
probe_info <- mappedkeys(probe_map)
length(probe_info)
[1] 11460
#探针与基因产生映射的数目
gene_info <- as.list(probe_map[probe_info])
# 转化为数据表
length(gene_info)
[1] 11460
gene_symbol <- toTable(probe_map[probe_info])
# 从hgu95av2SYMBOL文件中，取出有映射关系的探针，并生成数据框给gene_symbol
head(gene_symbol)
   probe_id  symbol
1   1000_at   MAPK3
2   1001_at    TIE1
3 1002_f_at CYP2C19
4 1003_s_at   CXCR5
5   1004_at   CXCR5
6   1005_at   DUSP1

mappedkeys用法示例，帮助理解。

library(hgu95av2.db)
  x <- hgu95av2GO
  x
  length(x)
  count.mappedkeys(x)
  x[1:3]
  links(x[1:3])

  ## Keep only the mapped keys
  keys(x) <- mappedkeys(x)
  length(x)
  count.mappedkeys(x)
  x # now it is a submap

  ## The above subsetting can also be achieved with
  x <- hgu95av2GO[mappedkeys(hgu95av2GO)]

  ## mappedkeys() and count.mappedkeys() also work with an environment
  ## or a list
  z <- list(k1=NA, k2=letters[1:4], k3="x")
  mappedkeys(z)
  count.mappedkeys(z)

  ## retrieve the set of primary keys for the ChipDb object named 'hgu95av2.db'
  keys <- keys(hgu95av2.db)
  head(keys)

9.整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。

A. 提示，理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
B. 然后根据得到探针去过滤原始表达矩阵

ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
View(head(exprSet))

10.把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。

rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
#  CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# MAPK3    5.743132  6.219412  5.523328  5.340477  5.229904
# TIE1     2.285143  2.291229  2.287986  2.295313  2.662170
# CYP2C19  3.309294  3.318466  3.354423  3.327130  3.365113
# CXCR5    1.085264  1.117288  1.084010  1.103217  1.074243
# CXCR5    7.544884  7.671801  7.474025  7.152482  6.902932

library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
#  probe    sample    value    group
#1   MAPK3 CLL11.CEL 5.743132 progres.
#2    TIE1 CLL11.CEL 2.285143 progres.
#3 CYP2C19 CLL11.CEL 3.309294 progres.
#4   CXCR5 CLL11.CEL 1.085264 progres.
#5   CXCR5 CLL11.CEL 7.544884 progres.
#6   DUSP1 CLL11.CEL 5.083793 progres.
View(head(exprSet))

11. 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ，然后画所有样本的这些图

参考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html

理解ggplot2的绘图语法，数据和图形元素的映射关系

### ggplot2
library(ggplot2)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)

image.png

12.理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。

注意：这个题目出的并不合规，请仔细看。

g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)

 [1] "DUSP5"   "IGFBP4"  "H1FX"    "ENPP2"   "FLNA"    "CLEC2B"  "TSPYL2"  "ZNF266"  "S100A9"  "NR4A2"   "TGFBI"   "ARF6"    "APBB2"   "VCAN"    "RBM38"  
[16] "CAPG"    "PLXNC1"  "RGS2"    "RNASE6"  "VAMP5"   "CYBB"    "GNLY"    "CCL3"    "OAS1"    "ENPP2"   "TRIB2"   "ZNF804A" "H1FX"    "IGH"     "JUND"   
[31] "SLC25A1" "PCDH9"   "VIPR1"   "COBLL1"  "GUSBP11" "S100A8"  "HBB"     "FOS"     "LHFPL2"  "FCN1"    "ZAP70"   "IGLC1"   "LGALS1"  "HBB"     "FOS"    
[46] "SLAMF1"  "TCF7"    "DMD"     "IGF2BP3" "FAM30A"

13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。

## heatmap
library(pheatmap)
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

image.png

14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。

## UpSetR
# https://cran.r-project.org/web/packages/UpSetR/README.html
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0),
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)

image.png

15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。

pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
# [1] "progres." "stable"   "progres." "progres." "progres." "progres." "stable"   "stable"   "progres." "stable"   "progres." "stable"   "progres." "stable"  
# [15] "stable"   "progres." "progres." "progres." "progres." "progres." "progres." "stable"  
dim(exprSet)
# [1] 8585   22
exprSet[1:5,1:5]
#      CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
MAPK3    5.743132  6.219412  5.523328  5.340477  5.229904
TIE1     2.285143  2.291229  2.287986  2.295313  2.662170
CYP2C19  3.309294  3.318466  3.354423  3.327130  3.365113
CXCR5    7.544884  7.671801  7.474025  7.152482  6.902932
DUSP1    5.083793  7.610593  7.631311  6.518594  5.059087

16.对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)

## hclust
colnames(exprSet)=paste(group_list,1:22,sep='')
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
                cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)

image.png

17.对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。

# install.packages("ggfortify")
library(ggfortify)
exprSet <- exprs(sCLLex)
df <- as.data.frame(t(exprSet))
df$group <- group_list
# autoplot uses ggplot2 to draw a particular plot for an object of a particular class in a single command.
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group')

image.png

18.根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格

## t.test
dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
#           avg_1    avg_2     log2FC        pvals     p.adj
36129_at 7.875615 8.791753  0.9161377 1.629755e-05 0.2057566
37676_at 6.622749 7.965007  1.3422581 4.058944e-05 0.2436177
33791_at 7.616197 5.786041 -1.8301554 6.965416e-05 0.2436177
39967_at 4.456446 2.152471 -2.3039752 8.993339e-05 0.2436177
34594_at 5.988866 7.058738  1.0698718 9.648226e-05 0.2436177
32198_at 4.157971 3.407405 -0.7505660 2.454557e-04 0.3516678

19.使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图)。

# DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix
##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
##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)

## volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)

g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)

image.png

20.对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。

### different P values
head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1])
plot(DEG_t.test[,4],nrDEG[,4])
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))

image.png

rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
exprSet['GAPDH',]
exprSet['ACTB',]
exprSet['DLEU1',]
library(ggplot2)
library(ggpubr)
my_comparisons <- list(
  c("stable", "progres.")
)
dat=data.frame(group=group_list,
               sampleID= names(exprSet['DLEU1',]),
               values= as.numeric(exprSet['DLEU1',]))
ggboxplot(
  dat, x = "group", y = "values",
  color = "group",
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")

image.png

## heatmap
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

image.png

生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！
B站链接：https://m.bilibili.com/space/338686099
YouTube链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程师入门最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
学徒培养：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw
生信技能树 - 简书 https://www.jianshu.com/u/d645f768d2d5

生信人的20个R语言习题及其答案-土豆学习笔记

1. 安装一些R包：

2.了解ExpressionSet对象

3.了解 str,head,help函数，作用于第二步提取到的表达矩阵

4. 安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后显示的那些变量

5. 理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID

6.理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？

7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

%in% 逻辑判断

8.过滤表达矩阵，删除那1165个没有对应基因名字的探针。

方法1：%in% 逻辑判断

方法2 mappedkeys() 映射关系

9.整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。

10.把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。

11. 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ，然后画所有样本的这些图

12.理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。

13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。

14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。

15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。

16.对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)

17.对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。

18.根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格

19.使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图)。

20.对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

R语言做生信

R语言在生物信息学的应用

作图

数据-R语言-图表-决策-Linux-Python

R生信

R语言与统计分析

Cook R

生物信息学与算法

生信人的20个R语言习题及其答案-土豆学习笔记

1. 安装一些R包：

2.了解ExpressionSet对象

3.了解 str,head,help函数，作用于 第二步提取到的表达矩阵

4. 安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 显示的那些变量

5. 理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID

6.理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？

7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

%in% 逻辑判断

8.过滤表达矩阵，删除那1165个没有对应基因名字的探针。

方法1：%in% 逻辑判断

方法2 mappedkeys() 映射关系

9.整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。

10.把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。

11. 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ， 然后画所有样本的 这些图

12.理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。

13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。

14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。

15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。

16.对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)

17.对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。

18.根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格

19.使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图)。

20.对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

3.了解 str,head,help函数，作用于第二步提取到的表达矩阵

4. 安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后显示的那些变量

11. 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ，然后画所有样本的这些图