RNASEQ（二）: PCA主成分分析

对RNAsq的read count数据进行PCA分析

目的：PCA分析可以得到样本之间的相关性和离散程度。

内容：
1 . 基因表达量数据进行标准化，用tpm和fpkm两种方法进行相对定量，后续分析我们一般会用tpm。

2 . 使用标准化后的tpm数据做主成分分析（PCA）

数据：RNASEQ上游分析得到的read count矩阵。

工具：Rstudio。

步骤：

1. 下载人类基因组注释文件，提取每个基因的外显子长度，计算tpm,fpkm。

TPM=(Ni/Li)*1000000/sum(Ni/Li+……..+ Nm/Lm)

Ni：mapping到基因i上的read数； Li：基因i的外显子长度的总和

setwd("E:/")
#注释文件获取
txdb <- makeTxDbFromGFF("Homo_sapiens.GRCh38.98.chr.gtf", format="gtf")
exons_gene <- exonsBy(txdb,by="gene")
exons_gene_lens <- lapply(exons_gene, function(x){sum(width(reduce(x)))})
# lappy ()返回一个长度与输入列表对象相似的列表
#查看下
class(exons_gene_lens)
length(exons_gene_lens)
#转换成data frame
exons_gene_lens1 <- as.data.frame(exons_gene_lens)
class(exons_gene_lens1)
dim(exons_gene_lens1)
#转置
exons_gene_lens1 <- t(exons_gene_lens1)
dim(exons_gene_lens1)

##将文件输出
write.csv(exons_gene_lens1, file="gene_Length.csv")
##此时基因名字后面还存在小数点，所以需要转换一下。
#然后替换名字
#也可以不替换，直接跳过之一步骤。根据你自己的需求。
exons_gene_lens2 <- read.csv("gene_Length.csv", header=T)
colnames(exons_gene_lens2) <- c("Geneid","Length")
exons_gene_lens2[1:3,1:2]

write.csv(exons_gene_lens2, file="gene_Length.csv")
##载入reads counts的原始数据。
readscounts <- read.table("E:/原始数据/counts_matrix.txt", header=T)

##如果你自己的矩阵名字更刚才的gene_Length.csv一样，就不需要。
###将基因长度和原始矩阵合并
##取出两个样本之间都有的，也就是交集。
mycounts <- merge(exons_gene_lens2, readscounts, by="Geneid", all=F)
dim(mycounts)
mycounts[1:3,1:3]

##将第一列变成列名
##然后删除第一列。
rownames(mycounts)<-mycounts[,1]
mycounts<-mycounts[,-1]
mycounts[1:3,1:3]

##TPM计算
kb <- mycounts$Length / 1000
write.csv(kb, file="kb.csv")
##筛选出样本的表达数据
# 这一步骤是筛选出所有样本的表达量数据，第一列是length，所以不要。
countdata <- mycounts[,2:7]
##计算每个值的rpk
rpk <- countdata / kb

##计算TPM
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)
write.csv(tpm, file="mRNA_tpm.csv")
##计算FPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6) 
head(fpkm)
write.csv(fpkm, file="mRNA_fpkm.csv")

2. 主成分分析并作图。代码如下

tpm <- read.csv("E:/造血干/钟/mRNA_tpm.csv")
#使用prcomp函数进行PCA分析
# 数据转置，转换成行为样本，列为基因的矩阵
tpm <- t(tpm)
data.pca <- prcomp(tpm)
# 查看PCA的分析结果
summary(data.pca)
# 绘制主成分的碎石图
screeplot(data.pca, npcs = 10, type = "lines")

#查看主成分的结果
pca.scores <- as.data.frame(data.pca$x)
head(pca.scores)

####主成分分析可视化
library(ggpubr)
library(ggthemes)
#BiocManager::install("gmodels")
library(gmodels)
library(export)

setwd("E:/原始数据")
#导入计算完成的tpm数据
mytpm <- read.csv("mRNA_tpm.csv", header=TRUE)
##修改列名
rownames(mytpm) <-mytpm[,1]
mytpm <-mytpm[,-1]

#用gmodels里面自带的fast.prcomp对PC进行计算。
pca.info <- fast.prcomp(mytpm)
#查看PCA的结果。
head(pca.info)
summary(pca.info) #能看到每个PC具体的解释比例
head(pca.info$rotation) #特征向量，回归系数
head(pca.info$sdev) #特征值的开方
head(pca.info$x) #样本得分score

#将PCA的结果和数据的表型构建为一个数据矩阵。
##这个里面的Type=c(rep("control",50),rep("case",374))和列表里面样本的名字是对应的，根据自己数据的情况灵活选择。
pca.data <- data.frame(sample = rownames(pca.info$rotation), 
                       Type=c(rep("D4",2),rep("D8",2),rep("D12",2)), pca.info$rotation)

#绘图，不同的参数会有不同的结果。具体的参数以及含义自行百度。以下是两个例子。
ggscatter(pca.data,x="PC1", y="PC2", color="Type", ellipse=TRUE, ellipse.type="convex")
ggscatter(pca.data,x="PC1", y="PC2", color="Type", ellipse=TRUE, 
          size=4, palette=c("#1d419b","#CC0000"), 
          label="sample", repel=TRUE, main="PCA plot"+theme_base())


m=c("D4_1","D4_2","D8_1","D8_2","D12_1","D12_2")
##我最终画图用的。          
ggscatter(pca.data,x="PC1", y="PC2", color="Type", 
          size=4, ellipse.border.remove=TRUE,
          label=c(m), repel=TRUE, main="PCA plot"+theme_base())
#输出到文件
graph2ppt(file="PCA_2D.ppt", width=10, aspectr=1.0)

3. 结果

   
   image.png