rm(list=ls())
rm(list=ls())用于清除所有的變量，ls() 會(huì)返回當(dāng)前環(huán)境里所有的objects的名字迁客，list是rm() 里邊的options之一郭宝，意思是“想要?jiǎng)h除的變量的名字是。掷漱≌呈遥”。rm(list=ls()) 意思是我要?jiǎng)h除變量卜范，刪除變量的名字是ls()衔统。

讀數(shù)據(jù)的命令出錯(cuò)

BiocManager::install("data.table")
library(data.table)
用這個(gè)包讀數(shù)據(jù)很快，而且不會(huì)有很多格式的問題

• 讀數(shù)據(jù)的另一個(gè)命令
  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) 將數(shù)據(jù)儲(chǔ)存為txt格式后再讀

• 每次讀文件之前海雪，得告訴R你文件所在的位置锦爵，setwd("")
  setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop")

看數(shù)據(jù)的話，是fix(data)奥裸，或者h(yuǎn)ead(data[1:10,1:10])

fix(data)

data[1:10,1:10]

data<-fread("RNA-seq - 1.txt")
輸入數(shù)據(jù)：data<-fread("文件名.txt")
將xlsx數(shù)據(jù)另存為txt格式险掀。

火山圖

data$color <- ifelse(data$padj<0.05 & abs(data$log2FoldChange)>= 1,ifelse(data$log2FoldChange > 1,'red','blue'),'gray')
# abs：絕對值，對于有生物學(xué)重復(fù)的樣本湾宙，要求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

• 控制坐標(biāo)軸范圍：
  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圖

根據(jù)基因表達(dá)值進(jìn)行樣本間的PCA分析樟氢，確定樣本在PCA圖中的位置，R語言中能夠執(zhí)行PCA分析的方法有很多侠鳄，不過它們的算法都是統(tǒng)一的埠啃，隨便使用任何一個(gè)R包就可以

# 讀取數(shù)據(jù)
data <- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
# 提取第一列   (不用這一條也行）
symbol<-as.matrix(gene[,1])
# 選2-11列（自動(dòng)將第一列當(dāng)列名）
exp = as.matrix(data[,c(1:8)])
# 要用 log 轉(zhuǎn)化后的基因表達(dá)值，降低不同基因表達(dá)水平數(shù)量級相差過大的問題
exp<-log2(exp+1)
#我們使用 FactoMineR 包中的方法伟恶，實(shí)現(xiàn) PCA 分析和聚類添加
library(FactoMineR)
#樣本中基因表達(dá)值的 PCA 分析
exp.pca <- PCA(exp, ncp = 2, scale.unit = TRUE, graph = FALSE)
plot(exp.pca)

• 讀取數(shù)據(jù)遇到的一些問題：
  當(dāng)我查看數(shù)據(jù)時(shí)碴开，前面都是沒有內(nèi)容的，還是RNA-seq的數(shù)據(jù)

  
  
  
  

  代碼如下：

  
  
  
  
  
  解決如下：在運(yùn)行標(biāo)紅的代碼前博秫，只要截取前8列的表達(dá)譜數(shù)據(jù)即可潦牛，因?yàn)樽鯬CA圖需要的數(shù)據(jù)只有前8列。
  

上一步獲得了PCA分析結(jié)果挡育，并觀察到明顯的組間差異巴碗。接下來使用ggplot2包繪制一個(gè)好看的PCA圖。
需要事先在上述PCA分析結(jié)果中將關(guān)鍵信息提取出静盅，例如樣本點(diǎn)在PCA圖中的位置信息良价，以及PCA軸的貢獻(xiàn)度等寝殴。

#提取 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 前兩軸的貢獻(xiàn)度
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)  #作圖數(shù)據(jù)中包含了樣本坐標(biāo)和分組信息

#ggplot2 繪制二維散點(diǎn)圖
library(ggplot2)

p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 3) +  #根據(jù)樣本坐標(biāo)繪制二維散點(diǎn)圖
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')) +  #去除背景和網(wǎng)格線
labs(x =  paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '')  #將 PCA 軸貢獻(xiàn)度添加到坐標(biāo)軸標(biāo)題中
p

• 若出現(xiàn)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讀數(shù)據(jù)，data.table=F痊银，這樣也是數(shù)值型

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])  #計(jì)算主成分
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

• 一些報(bào)錯(cuò)的解決辦法
  data.pca <- prcomp(data)
  Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
  解決辦法：
  data<-as.numeric(as.matrix(data))
  這個(gè)錯(cuò)誤是在使用read.table()函數(shù)讀取文件遇到的抵蚊，經(jīng)過摸索，發(fā)現(xiàn)read.table()讀取的數(shù)據(jù)類型為data.frame, 通過將y的類型由data.frame轉(zhuǎn)換為numeric可解決報(bào)錯(cuò)溯革。
• 一些命令含義
  將基因表達(dá)值矩陣作個(gè)轉(zhuǎn)置贞绳，使行為樣本，列為基因

dim(df) # 查看行和列數(shù)
[1] 603  21
tmpp<-apply(tmpp,2,as.numeric)  # 2表示按列轉(zhuǎn)數(shù)值類型
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構(gòu)建標(biāo)簽

ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+ geom_point()

ggplot

去掉背景及網(wǎng)格線

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")