R语言高级题作业笔记
2019-06-12 本文已影响74人
Michael_fzy
题目链接，生信菜鸟团：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
