高级题目
1.安装一些R包:
if (!require("ALL")) {BiocManager::install("ALL")}
if (!require("CLL")) {BiocManager::install("CLL")}
if (!require("pasilla")) {BiocManager::install("pasilla")}
if (!require("airway")) {BiocManager::install("airway")}
if (!require("limma")) {BiocManager::install("limma")}
if (!require("DESeq2")) {BiocManager::install("DESeq2")}
if (!require("clusterProfiler")) {BiocManager::install("clusterProfiler")}
if (!require("reshape2")) {BiocManager::install("reshape2")}
以及对这些包的一些介绍:
1、ALL: DFCI丽兹实验室的T细胞和B细胞急性淋巴细胞白血病数据(包括2004年4月的版本)
2.CLL: CLL软件包包含慢性淋巴细胞性白血病(CLL)基因表达数据。CLL数据包含24个样本,这些样本在疾病进展方面被分类为进行性或稳定。
3、pasilla:该软件包提供了从RNA-seq数据中选择的基因计算出的每个外显子和每个基因的读数计数,这些数据在Brooks AN,Yang L,Duff MO,Hansen的文章“果蝇和哺乳动物之间的RNA调节图的保守性”中发表。
4、airway:该程序包提供了一个RangedSummarizedExperiment对象,该对象是基因的读取计数,用于在用地塞米松治疗的四种人呼吸道平滑肌细胞系上进行RNA-Seq实验。小插图包装中提供了有关基因模型和读数计数程序的详细信息。
5、limma: 数据分析,线性模型和微阵列数据的差异表达。
6、DESeq2:估计高通量测序分析中计数数据的方差均值依赖性,并基于使用负二项式分布的模型测试差异表达。
7、clusterProfiler:该软件包实现了分析和可视化基因和基因簇的功能概况(GO和KEGG)的方法。
8、reshape2:使用两个函数灵活地重组和聚合数据:melt和'dcast'(或'acast')。
2.了解ExpressionSet对象,比如CLL包里面就有data(sCLLex) ,找到它包含的元素,提取其表达矩阵(使用exprs函数),查看其大小
1.什么是ExpressionSet对象?
在Biobase基础包中,ExpressionSet是非常重要的类,因为 Bioconductor设计之初是为了对基因芯片数据进行分 析,而ExpressionSet正是Bioconductor为基因表达数据格式所定制的标准。它是所有涉及基因表达量相关数据在 Bioconductor中进行操作的基础数据类型,比如affyPLM, affy, oligo, limma, arrayMagic等等
ExpressionSet的组成:assayData--头文件--experimentData
suppressPackageStartupMessages(library(CLL))#加载包
data(sCLLex)
exprSet=exprs(sCLLex) ##提取表达矩阵,它列是各个样本,行是每个探针ID
samples=sampleNames(sCLLex) ##提取名字
pdata=pData(sCLLex) ##提取临床信息
group_list=as.character(pdata[,2]) ##强制转换类型,转换为字符型
dim(exprSet) ##查看维度.
[1] 12625 22 #行列
exprSet[1:5,1:5] #查看
###再做一个QC检测
par(cex = 0.7) #控制图片的文字和点的大小.
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main="expression value",las=2)做一个箱形图,可以查看一些特异的点
3.了解 str,head,help函数,作用于 第二步提取到的表达矩阵
> help('str') #str()函数来查看这个数据框当中有哪些变量,以及变量的类型
> help('head') #查看前面部分
> help('help') #帮助函数
str(exprSet)
head(exprSet)
4.安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 显示的那些变量
Affymetrix人类基因组U95设置注释数据(芯片hgu95av2),使用来自公共存储库的数据组装而成
BiocManager::install("hgu95av2.db")
library('hgu95av2.db')
ls("package:hgu95av2.db")
str("package:hgu95av2.db")
5.理解 head(toTable(hgu95av2SYMBOL)) 的用法,找到 TP53 基因对应的探针ID
ls("package:hgu95av2.db") #查看hgu95av2有什么变量
head(toTable(hgu95av2SYMBOL)) #查看hgu95av2SYMBOL的前面部分的数据
ids <- toTable(hgu95av2SYMBOL) #换为表格
TP53 <- ids[ids$symbol=='TP53',] #选择需要的内容
6.理解探针与基因的对应关系,总共多少个基因,基因最多对应多少个探针,是哪些基因,是不是因为这些基因很长,所以在其上面设计多个探针呢?
①.探针和基因的对应关系:
在用基因芯片中,是通过一系列已知的核酸探针,然后和带有荧光标记的核酸序列互补配对,通过测定荧光强度最强的位置,获取一组互补的探针序列.后期的分析,通过探针组的结果,去分析对应的基因.一一对应.
②.总共多少个基因,基因最多对应多少个探针,是哪些基因
attach(ids) #这样后面访问的时候,就不用加$了
length(unique(symbol)) #查看长度,也即是个数
[1] 8584 #一共有8584个基因symbol
table(symbol)[table(symbol)==max(table(symbol))] 查找众数,也就是出现最多的元素。table可以用来对数据统计出现的次数。
7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵,找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。
noin_exprSet <- exprSet[!row.names(exprSet)%in%ids$probe_id,] #首先是提取exprSet的行名字,也就是探针,然后用%in%进行一个逻辑判断,是否相同.probe_id是在hgu95av2SYMBOL的行名.然后一比对
dim(noin_exprSet)
[1] 1166 22
8.过滤表达矩阵,删除那1165个没有对应基因名字的探针。
exprSet <- exprSet[rownames(exprSet) %in% ids$probe_id, ] ##原理和前面一题一样。!是取反向操作,这题不用
dim(exprSet)
9.整合表达矩阵,多个探针对应一个基因的情况下,只保留在所有样本里面平均表达量最大的那个探针。
提示,理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
然后根据得到探针去过滤原始表达矩阵
思路:先把是多个探针对应一个基因的那些探针找出来——然后再比较表达量——再过滤。X
参考答案,应该用by函数,进行分类,把多个探针对应一个基因的这样的情况分类出来——再比较表达量
a <- by(exprSet,ids$symbol, ##exprSet是需要处理的数据,ids$symbol是分类的标注。分类出每一个小的data
function(x) ##对每一个小的data进行函数处理。x就应该是代表着每一个分类出来的data
row.names(x)[which.max(rowMeans(x))]) ##先去平均值,再找最大的值,过滤开来。
head(a) ##看出a还是有探针和基因对应的,但是我们只要探针就可以了。
INDICES
AADAC AAK1 AAMP AANAT AARS1 AASDHPPT
"36512_at" "39463_at" "38434_at" "36332_at" "36185_at" "35761_at"
b = as.character(a) ##转换数据类型
[1] "36512_at" "39463_at" "38434_at" "36332_at" "36185_at" "35761_at" ##这样就是有探针了,后面可以用来对exprSet进行筛选。
exprSet <- exprSet[row.names(exprSet)%in%b,] 从前面选出来的b,对exprSet进行筛选。
dim(exprSet)
[1] 8584 22
10.把过滤后的表达矩阵更改行名为基因的symbol,因为这个时候探针和基因是一对一关系了。
row.names(exprSet)=ids[rownames(exprSet)%in%ids$probe_id,2] ##逻辑判断一下,然后更换行名。注意这个“2”。应该是要去第二列的意思。
11、对第10步得到的表达矩阵进行探索,先画第一个样本的所有基因的表达量的boxplot,hist,density , 然后画所有样本的 这些图
参考:http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
理解ggplot2的绘图语法,数据和图形元素的映射关系
思路:先画第一个样本的图:看着参考,要画的图有:boxplot、vioplot、boxplot、histogram
density、gpairs、cluster、PCA、heatmap、DEG && volcano plot————再看看能不能做个循环,把所有的样本都画出来。
##箱形图
p=ggplot(long_exprSet,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
image.png
大概是这样的图。x轴是22个样本,y轴是表达量。不同颜色是每样本人的疾病进展情况(progres和stable)
可以得出的结论是:
①、看出每个样本的异常值的情况。有几个样本的异常值特别多,比如16号样本??有什么用?
②、看出平均值,每个样本的表达量的平均值的??
###vioplot
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
image.png
什么是vioplot?
vioplot叫做小提琴图。小提琴图 (Violin Plot) 用于显示数据分布及其概率密度。
得出什么结论?
①、????
##直方图
p=ggplot(long_exprSet,aes(value,fill=group))+geom_histogram()+facet_wrap(~sample,nrow = 3)##做个直方图,其中facet_wrap是一个分面函数。
print(p)
image.png
得出的结论:
①、大部分的样本,的表达量分布在3-4左右。
##密度图
p=ggplot(long_exprSet,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
image.png
可以看出,表达量value值,在大概3-4的时候,fre频率是最大的。
gpairs????看不懂怎么看这图呀,相关性分析?
12.理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标,最后按照mad值排序,取top 50 mad值的基因,得到列表。
注意:这个题目出的并不合规,请仔细看。
思路:可以很简单的计算出一个基因的这些统计指数,如何计算多个基因的呢??
——所以需要用到apply函数才可以
mean <- apply(exprSet, 1,mean)
mean <- as.data.frame(mean)
median <- apply(exprSet, 1,median)
median <- as.data.frame(median)
max <- apply(exprSet, 1,max)
max <- as.data.frame(max)
min <- as.data.frame(apply(exprSet, 1,min))
var <- as.data.frame(apply(exprSet, 1,var))
mad <- as.data.frame(apply(exprSet, 1,mad))
merge(mean,median,max,min,sd,var,mad,by=rownames)??想合并起来,看起来会比较好看,不过好像不行,因为dim()看出来,只有一列。还是算了,不合并了。
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
names(g_mad)
13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集,并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。
思路:先根据g_mad来在exprest里面取出部分数据——然后再对这份数据进行做个热图。
TOP50 <- exprSet[g_mad,]
library(pheatmap)
pheatmap(TOP50)
image.png
14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表,使用UpSetR包来看他们之间的overlap情况。
UpSetR包:对集合,进行可视化的一个包。
思路:怎么去做一个,符合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), ##ifelse的一个循环。
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)
dim(pdata)
Disease <- pdata$Disease
print(Disease)
16.对所有样本的表达矩阵进行聚类并且绘图,然后添加样本的临床表型数据信息(更改样本名)
思路:如何聚类?——绘图?——更改样本名。
怎么看聚类图。从大的往里面切,看看自己分几类
那分类的依据是什么?不知道用的是什么算法,但分类的依据应该是,基因表达量。
colnames(exprSet)=paste(group_list,1:22,sep='') ##把CLL11.CEL——换为progres.这样的病例信息
exprSet <- t(exprSet) ##需要吧行和列转换一下。因为接下来要多个dist,求距离的,是对dist的是计算每一行的之间的距离,然后是以每一列的作为变量作为标注
hc <- hclust(dist(exprSet)) ##聚类
plot(as.dendrogram(hc), horiz = TRUE) ##horiz = TRUE是吧竖向的换为横向的。
image.png
怎么看这个图:从左往右看,看看那些样本是一大类的。——这样又能得出什么有用的信息呢?
17.对所有样本的表达矩阵进行PCA分析并且绘图,同样要添加表型信息。
思路:什么是PCA分析。ggfortify这个包怎么用。
PCA分析:主成分分析,数据的降维?看疾病组和对照组能否分开
install.packages("ggfortify")
library(ggfortify)
exprSet <- exprs(sCLLex)
df <- as.data.frame(t(exprSet))
df$group <- group_list
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group') ???
image.png
第一成分和第二成分加起来才30%左右,证明这个是不太好的?
另外一个就是progres组和stable组是没有明显差异的,因为他们有交集,距离比较近。是否有差异还得看p值
18.根据表达矩阵及样本分组信息进行批量T检验,得到检验结果表格
如何进行T检验?
还是不懂??双样本的检验?
19.使用limma包对表达矩阵及样本分组信息进行差异分析,得到差异分析表格,重点看logFC和P值,画个火山图(就是logFC和-log10(P值)的散点图)。
??如何用这个包?
对哪些数据进行分析?
suppressMessages(library(limma)) #使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组。这时,如果逐一使用两两比较求差异基因则略显复杂。其实开发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)
fit
##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) 去掉NA值
head(nrDEG )
#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)