rm(list=ls())
rm(list=ls())用于清除所有的变量，ls() 会返回当前环境里所有的objects的名字，list是rm() 里边的options之一，意思是“想要删除的变量的名字是。。”。rm(list=ls()) 意思是我要删除变量，删除变量的名字是ls()。
读数据的命令出错
BiocManager::install("data.table")
library(data.table)
用这个包读数据很快，而且不会有很多格式的问题
读数据的另一个命令
data <- read.table(choose.files(),header = TRUE,sep="\t")
data <- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)将数据储存为txt格式后再读每次读文件之前，得告诉R你文件所在的位置，setwd("")
setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop")
看数据的话，是fix(data)，或者head(data[1:10,1:10])
fix(data)
data[1:10,1:10]
data<-fread("RNA-seq - 1.txt")
输入数据：data<-fread("文件名.txt")
将xlsx数据另存为txt格式。
火山图
data$color <- ifelse(data$padj<0.05 & abs(data$log2FoldChange)>= 1,ifelse(data$log2FoldChange > 1,'red','blue'),'gray')
# abs：绝对值，对于有生物学重复的样本，要求padj﹤0.05
color <- c(red = "red",gray = "gray",blue = "blue")
p <- ggplot(data, aes(log2FoldChange, -log10(padj), col = color)) +
geom_point() +
theme_bw() +
scale_color_manual(values = color) +
labs(x="log2 (fold change)",y="-log10 (q-value)") +
geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
theme(legend.position = "none",
panel.grid=element_blank(),
axis.title = element_text(size = 16),
axis.text = element_text(size = 14))
p
控制坐标轴范围：
p+scale_x_continuous(limits = c(-3,3))+scale_y_continuous(limits = c(0,5)) xlim选最大值作为xlim的上下边界
x_lim <- max(logFC,-logFC)
y_lim<-max(-1*log10(padj))
- 确定x轴最大值：
colnames(data)
[1] "gene_id" "S10" "S2NAC"
[4] "S7" "S9" "M3"
[7] "M7" "M8" "M9"
[10] "baseMeanA" "baseMeanB" "baseMean"
[13] "log2FoldChange" "lfcSE" "stat"
[16] "pvalue" "padj" "type"
[19] "GO_BP" "GO_MF" "GO_CC"
[22] "KEGG_PATHWAY" "color"
> max(abs(data[,13]))
[1] 23.22173
- 绘制火山图
library(ggplot2)
library(RColorBrewer)
theme_set(theme_bw())
p <- ggplot(data,aes(log2FoldChange,-1*log10(padj),
color = sig))+geom_point(aes(color=sig))+
xlim(-23.22,22.23) + labs(x="log2 Fold Change",y="-log10(pvalue)")
p <- p + scale_color_manual(values =c("#0072B5","grey","#BC3C28"))+
geom_hline(yintercept=-log10(0.05),linetype=4)+
geom_vline(xintercept=c(0),linetype=4)
p <- p +theme(axis.line = element_line(size=0))+ylim(0,y_lim)
p <- p +theme(axis.text=element_text(size=20),axis.title=element_text(size=20))
p
PCA图
根据基因表达值进行样本间的PCA分析，确定样本在PCA图中的位置，R语言中能够执行PCA分析的方法有很多，不过它们的算法都是统一的，随便使用任何一个R包就可以
# 读取数据
data <- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
# 提取第一列 (不用这一条也行）
symbol<-as.matrix(gene[,1])
# 选2-11列（自动将第一列当列名）
exp = as.matrix(data[,c(1:8)])
# 要用 log 转化后的基因表达值，降低不同基因表达水平数量级相差过大的问题
exp<-log2(exp+1)
#我们使用 FactoMineR 包中的方法，实现 PCA 分析和聚类添加
library(FactoMineR)
#样本中基因表达值的 PCA 分析
exp.pca <- PCA(exp, ncp = 2, scale.unit = TRUE, graph = FALSE)
plot(exp.pca)
-
读取数据遇到的一些问题：
当我查看数据时，前面都是没有内容的，还是RNA-seq的数据
代码如下：
解决如下：在运行标红的代码前，只要截取前8列的表达谱数据即可，因为做PCA图需要的数据只有前8列。
上一步获得了PCA分析结果，并观察到明显的组间差异。接下来使用ggplot2包绘制一个好看的PCA图。
需要事先在上述PCA分析结果中将关键信息提取出，例如样本点在PCA图中的位置信息，以及PCA轴的贡献度等。
#提取 PCA 中的样本位置，即x轴（pca1轴）和y轴（pca2轴）
pca_sample <- data.frame(exp.pca$ind$coord[ ,1:2])
head(pca_sample)
Dim.1 Dim.2
Ciart 8.186812 -0.6273470
Fabp7 7.368796 0.5651271
Fmo2 4.965070 -1.0870700
Hif3a 7.467280 -0.7283460
Kirrel2 7.125757 -0.4987601
Plin4 6.743315 -0.8912442
#提取 PCA 前两轴的贡献度
pca_eig1 <- round(exp.pca$eig[1,2], 2)
pca_eig2 <- round(exp.pca$eig[2,2],2 )
读取并合并样本分组信息
group <- read.delim('pca.xlsx', row.names = 1, sep = '\t', check.names = FALSE)
group <- group[rownames(pca_sample), ]
pca_sample <- cbind(pca_sample, group)
pca_sample$samples <- rownames(pca_sample)
head(pca_sample) #作图数据中包含了样本坐标和分组信息
#ggplot2 绘制二维散点图
library(ggplot2)
p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 3) + #根据样本坐标绘制二维散点图
scale_color_manual(values = c('orange', 'purple')) + #自定义颜色
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) + #去除背景和网格线
labs(x = paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '') #将 PCA 轴贡献度添加到坐标轴标题中
p
- 若出现Error: unexpected symbol in "\u200b" ，将#注释去掉。
重新作图：
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1)
iris.pca<-prcomp(t(df[,1:8]))
summary(iris.pca)
plot(iris.pca$x,col=rep(c("red","blue"),c(4,4)),pch=19)
PCA
用fread读数据，data.table=F，这样也是数值型
library("data.table")
da<-fread("RNA-seq - 1.txt",data.table=F)
getwd()
"C:/Users/Administrator.DESKTOP-4UQ3Q0K/Desktop"
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
iris.pca<-prcomp(df[,1:8]) #计算主成分
iris_pcs <-data.frame(iris.pca$x)
head(iris_pcs)
PC1 PC2 PC3 PC4 PC5
Ciart -4934.580 5217.2257 1760.1106 -248.6604 262.260206
Fabp7 -4472.953 1544.6017 -349.0215 595.1240 -8.948593
Fmo2 -6893.075 676.7536 266.6638 -213.0876 132.676879
Hif3a -5814.847 3522.1495 1007.5449 -571.4270 489.318537
Kirrel2 -6008.083 2239.0950 625.8803 17.2822 -132.076437
Plin4 -6351.092 2396.1754 794.4072 -473.1440 365.604000
PC6 PC7 PC8
Ciart 161.56826 103.28029 -4.916623
Fabp7 -417.51164 -230.57635 -16.571842
Fmo2 19.73212 70.27110 -46.352493
Hif3a 169.11638 163.44406 -118.575331
Kirrel2 -34.05482 36.89955 61.624917
Plin4 172.99939 251.82729 61.521027
summary(iris.pca)
Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 1.640e+05 1.255e+03 639.23739 338.9 135 77.25
Proportion of Variance 9.999e-01 6.000e-05 0.00002 0.0 0 0.00
Cumulative Proportion 9.999e-01 1.000e+00 0.99999 1.0 1 1.00
PC7 PC8
Standard deviation 62.57 29.65
Proportion of Variance 0.00 0.00
Cumulative Proportion 1.00 1.00
# 做碎石图
library(ggplot2)
library(factoextra)
fviz_eig(iris.pca)
fviz_eig(iris.pca,addlabels = T)
碎石图
plot(exp.pca$x,cex = 2,main = "PCA analysis",
col = c(rep("red",3),rep("blue",3)),
pch = c(rep(16,3),rep(17,3)))
plot(exp.pca)
PCA analysis
- 一些报错的解决办法
data.pca <- prcomp(data)
Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
解决办法：
data<-as.numeric(as.matrix(data))
这个错误是在使用read.table()函数读取文件遇到的，经过摸索，发现read.table()读取的数据类型为data.frame, 通过将y的类型由data.frame转换为numeric可解决报错。 - 一些命令含义
将基因表达值矩阵作个转置，使行为样本，列为基因
dim(df) # 查看行和列数
[1] 603 21
tmpp<-apply(tmpp,2,as.numeric) # 2表示按列转数值类型
dim(tmpp)
[1] 8 603
按教程画一遍：
library(ggplot2)
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1)
head(df)
df_pca <- prcomp(df)
# Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
mode(df)
# "list"
df<-apply(df,2,as.numeric)
mode(df)
# "numeric"
df_pca<-prcomp(t(df[,1:8]))
df_pcs <-data.frame(df_pca$x, feature = rep(c("control","case"),c(4,4)))
# 用rep构建标签
ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+ geom_point()
ggplot
去掉背景及网格线
ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+
+ geom_point()+
+ theme_bw() +
+ theme(panel.border=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.line= element_line(colour = "black"))
image.png
添加PC1 PC2的百分比
percentage<-paste(colnames(df_pcs),"(", paste(as.character(percentage), "%", ")", sep=""))
ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+
+ geom_point(size=3)
+ xlab(percentage[1])
+ ylab(percentage[2])
image.png
使用ggplot2包绘制PCA图
library(ggplot2)
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1)
head(df)
df<-apply(df,2,as.numeric)
mode(df)
df_pca<-prcomp(t(df[,1:8]))
data.pca<-prcomp(t(df[,1:8]))
summary(data.pca)
pca.scores <- as.data.frame(data.pca$scores)
head(pca.scores)
library(ggplot2)
ggplot(pca.scores,aes(PC1,PC2,col=rownames(pca.scores)))
+ geom_point(size=3)
+ geom_text(aes(label=rownames(pca.scores)),vjust = "outward")
+ geom_hline(yintercept = 0,lty=2,col="red")
+ geom_vline(xintercept = 0,lty=2,col="blue",lwd=1)
+ theme_bw() + theme(legend.position = "none")
+ theme(plot.title = element_text(hjust = 0.5))
+ labs(x="PCA_1",y="PCA_2",title = "PCA analysis")
