题目链接，生信菜鸟团：http://www.bio-info-trainee.com/3409.html
高级题来来回回已经做了三遍了，之前一直没有把笔记放上来，今天就统一上传一下。

1 安装几个R包

数据包： ALL, CLL, pasilla, airway 
软件包：limma，DESeq2，clusterProfiler 
工具包：reshape2
绘图包：ggplot2

安装这些包网上都有详细的说明，这里就不再赘述了，主要是遇到问题要根据提示来派出问题，下面粘贴一段经典的代码

source("https://bioconductor.org/biocLite.R")#镜像
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite(c("ALL","CLL", "pasilla", "airway")) #数据包
biocLite(c("limma","DESeq2", "clusterProfiler")) #软件包
install.packages("reshape2") #工具包
install.packages("ggplot2") #绘图包

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

rm(list = ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data("sCLLex")#取数据
e=exprs(sCLLex) #提取出来表达矩阵，赋给e
str(e) # 查看结构
head(e) # 查看前6行
dim(e) #查看表达矩阵大小

其中的信息

> sCLLex
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 22 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL
    (22 total)
  varLabels: SampleID Disease
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2

> dim(e)
[1] 12625    22

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

该步在2中已经完成

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

> ls("package:hgu95av2.db")
 [1] "hgu95av2"              "hgu95av2.db"          
 [3] "hgu95av2_dbconn"       "hgu95av2_dbfile"      
 [5] "hgu95av2_dbInfo"       "hgu95av2_dbschema"    
 [7] "hgu95av2ACCNUM"        "hgu95av2ALIAS2PROBE"  
 [9] "hgu95av2CHR"           "hgu95av2CHRLENGTHS"   
[11] "hgu95av2CHRLOC"        "hgu95av2CHRLOCEND"    
[13] "hgu95av2ENSEMBL"       "hgu95av2ENSEMBL2PROBE"
[15] "hgu95av2ENTREZID"      "hgu95av2ENZYME"       
[17] "hgu95av2ENZYME2PROBE"  "hgu95av2GENENAME"     
[19] "hgu95av2GO"            "hgu95av2GO2ALLPROBES" 
[21] "hgu95av2GO2PROBE"      "hgu95av2MAP"          
[23] "hgu95av2MAPCOUNTS"     "hgu95av2OMIM"         
[25] "hgu95av2ORGANISM"      "hgu95av2ORGPKG"       
[27] "hgu95av2PATH"          "hgu95av2PATH2PROBE"   
[29] "hgu95av2PFAM"          "hgu95av2PMID"         
[31] "hgu95av2PMID2PROBE"    "hgu95av2PROSITE"      
[33] "hgu95av2REFSEQ"        "hgu95av2SYMBOL"       
[35] "hgu95av2UNIGENE"       "hgu95av2UNIPROT"

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

> a = toTable(hgu95av2SYMBOL)
> a[which(a$symbol=="TP53"),]
      probe_id symbol
966    1939_at   TP53
997  1974_s_at   TP53
1420  31618_at   TP53

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

> length(a$symbol)
[1] 11460
> length(unique(a$symbol))
[1] 8585
> tail(sort(table(a$symbol)))# 基因最多对应8个探针

YME1L1  GAPDH INPP4A    MYB PTGER3  STAT1 
     7      8      8      8      8      8 
> table(sort(table(a$symbol)))# 最多有6555个基因对应一个探针

   1    2    3    4    5    6    7    8 
6555 1428  451  102   22   16    6    5

为什么有5个基因会分别有8个探针，而大部分6555个基因只对应一个探针？
A：不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量

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

> n_e <- e[!(rownames(e) %in% a$probe_id),]
> dim(n_e)
[1] 1165   22

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

> e1 = e[rownames(e)%in%a$probe_id,]
> dim(e1)
[1] 11460    22

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

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

> maxid = by(e1,a$symbol,function(x) rownames(x)[which.max(rowMeans(x))])
> uniid = as.character(maxid)
> uni_e = e1[rownames(e1) %in% uniid,]
> str(uni_e)  # 8585
 num [1:8585, 1:22] 5.74 2.29 3.31 7.54 5.08 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:8585] "1000_at" "1001_at" "1002_f_at" "1004_at" ...
  ..$ : chr [1:22] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" ...

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

rownames(uni_e) = ids[match(rownames(uni_e),ids$probe_id),2]
suppressPackageStartupMessages(library(reshape2))
# match exprSet
m_e = melt(uni_e)
colnames(m_e) = c('symbol','sample','value')
pd = pData(sCLLex) # 查看样本分组情况
Disease = pd[,2]
m_e$group = rep(Disease,each=nrow(uni_e))

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

suppressPackageStartupMessages(library(ggplot2))
ggplot(m_e,aes(x=sample,y=value,fill=group)) + geom_boxplot() 
ggplot(m_e,aes(value,fill=group)) + geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
ggplot(m_e,aes(value,col=group)) + geom_density()

1.png

2.png

3.png

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

# 平均值、中位数、最大值、最小值、标准差、变量、中位数绝对偏差
e_mean = tail(sort(apply(uni_e,1,mean)),50)
e_median = tail(sort(apply(uni_e,1,median)),50)
e_max = tail(sort(apply(uni_e,1,max)),50)
e_min = tail(sort(apply(uni_e,1,min)),50)
e_sd = tail(sort(apply(uni_e,1,sd)),50)
e_var = tail(sort(apply(uni_e,1,var)),50)
e_mad = tail(sort(apply(uni_e,1,mad)),50)

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

# 做热图之前需要将数据中心化、标准化
top50_gene = tail(sort(apply(uni_e,1,mad)),50)
top50_matrix = uni_e[top50_gene,]
top50_matrix2 = t(scale(t(top50_matrix))) # scale() 对数据进行标准化
# 标准化是原始分数减去平均数然后除以标准差，中心化是原始分数减去平均数。一般流程为先中心化再标准化

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

suppressPackageStartupMessages(library("UpSetR"))
all <- unique(c(names(e_mean),names(e_median),names(e_max),names(e_min),
                  names(e_sd),names(e_var),names(e_mad)))
e_all = data.frame(all,
                   e_mean=ifelse(all %in% names(e_mean),1,0),
                   e_median=ifelse(all %in% names(e_median),1,0),
                   e_max=ifelse(all %in% names(e_max),1,0),
                   e_min=ifelse(all %in% names(e_min),1,0),
                   e_sd=ifelse(all %in% names(e_sd),1,0),
                   e_var=ifelse(all %in% names(e_var),1,0),
                   e_mad=ifelse(all %in% names(e_mad),1,0)
)
upset(e_all,nsets = 7,sets.bar.color = "#BD1111")

4.png

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

> group_list = as.character(pd[,2])
> table(group_list)
group_list
progres.   stable 
      14        8

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

clust = t(e)
rownames(clust) = colnames(e)
clust_dist = dist(clust,method = "euclidean")
hc = hclust(clust_dist,"ward.D")
suppressPackageStartupMessages(library(factoextra))
fviz_dend(hc, k = 4 ,cex = 0.5,k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE, rect = TRUE)

5.png

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

pca_data <- prcomp(t(e),scale=TRUE)
pcx <- data.frame(pca_data$x)
pcr <- cbind(samples=rownames(pcx),group_list, pcx) 
ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)

6.png

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

gl = as.factor(group_list)
group1 = which(group_list == levels(gl)[1])
group2 = which(group_list == levels(gl)[2]) 
et1 = e1[,group1]
et2 = e1[,group2]
data_t = cbind(et1,et2)
pvals = apply(e1, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")  
data_mean_1 = rowMeans(et1) 
#progres是对照组
data_mean_2 = rowMeans(et2) 
#stable是使用药物处理后的——处理组
log2FC = data_mean_2-data_mean_1
DEG_t.test = cbind(data_mean_1, data_mean_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)
         data_mean_1 data_mean_2     log2FC        pvals     p.adj
36129_at    7.875615    8.791753  0.9161377 1.629755e-05 0.1867699
37676_at    6.622749    7.965007  1.3422581 4.058944e-05 0.2211373
33791_at    7.616197    5.786041 -1.8301554 6.965416e-05 0.2211373
39967_at    4.456446    2.152471 -2.3039752 8.993339e-05 0.2211373
34594_at    5.988866    7.058738  1.0698718 9.648226e-05 0.2211373
32198_at    4.157971    3.407405 -0.7505660 2.454557e-04 0.3192169

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

suppressPackageStartupMessages(library(limma))
design1=model.matrix(~factor(group_list)) 
colnames(design1)=levels(factor(group_list))
rownames(design1)=colnames(e)
fit1 = lmFit(e,design1)
fit1=eBayes(fit1)
options(digits = 3) 
mtx1 = topTable(fit1,coef=2,adjust='BH',n=Inf) 
# topTable 默认显示前10个基因的统计数据；使用选项n可以设置,n=Inf就是不设上限，全部输出
DEG_mtx1 = na.omit(mtx1) #去除缺失值
head(DEG_mtx1)

火山图(标准代码，照搬就完事了）

DEG=DEG_mtx1
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) ) 
DEG$result = 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), #round保留小数位数
                    '\nThe number of up gene is ',nrow(DEG[DEG$result =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$result =='DOWN',])
)
ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=result)) +
  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'))

7.png

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

DEG_t.test = DEG_t.test[rownames(DEG_mtx1),]
plot(DEG_t.test[,3],DEG_mtx1[,1])#可以看到两组呈现负相关
plot(DEG_t.test[,4],DEG_mtx1[,4])#可以看到两组数据差异不大
plot(-log10(DEG_t.test[,4]),-log10(DEG_mtx1[,4]))

8.png

9.png

10.png

R语言高级题作业笔记

R语言高级题作业笔记

1 安装几个R包

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

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

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

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

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

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

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

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语言高级题作业笔记

1 安装几个R包

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

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

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

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

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

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

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

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") 后显示的那些变量