對(duì)RNAsq的read count數(shù)據(jù)進(jìn)行PCA分析
目的：PCA分析可以得到樣本之間的相關(guān)性和離散程度。
內(nèi)容：
1 . 基因表達(dá)量數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化翰铡，用tpm和fpkm兩種方法進(jìn)行相對(duì)定量赂蕴，后續(xù)分析我們一般會(huì)用tpm柳弄。
2 . 使用標(biāo)準(zhǔn)化后的tpm數(shù)據(jù)做主成分分析（PCA）
數(shù)據(jù)：RNASEQ上游分析得到的read count矩陣。
工具：Rstudio概说。
步驟：
- 下載人類(lèi)基因組注釋文件碧注，提取每個(gè)基因的外顯子長(zhǎng)度，計(jì)算tpm,fpkm糖赔。
TPM=(Ni/Li)*1000000/sum(Ni/Li+……..+ Nm/Lm)
Ni：mapping到基因i上的read數(shù)萍丐； Li：基因i的外顯子長(zhǎng)度的總和
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 ()返回一個(gè)長(zhǎng)度與輸入列表對(duì)象相似的列表
#查看下
class(exons_gene_lens)
length(exons_gene_lens)
#轉(zhuǎn)換成data frame
exons_gene_lens1 <- as.data.frame(exons_gene_lens)
class(exons_gene_lens1)
dim(exons_gene_lens1)
#轉(zhuǎn)置
exons_gene_lens1 <- t(exons_gene_lens1)
dim(exons_gene_lens1)
##將文件輸出
write.csv(exons_gene_lens1, file="gene_Length.csv")
##此時(shí)基因名字后面還存在小數(shù)點(diǎn)，所以需要轉(zhuǎn)換一下放典。
#然后替換名字
#也可以不替換逝变，直接跳過(guò)之一步驟。根據(jù)你自己的需求奋构。
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的原始數(shù)據(jù)壳影。
readscounts <- read.table("E:/原始數(shù)據(jù)/counts_matrix.txt", header=T)
##如果你自己的矩陣名字更剛才的gene_Length.csv一樣，就不需要弥臼。
###將基因長(zhǎng)度和原始矩陣合并
##取出兩個(gè)樣本之間都有的宴咧，也就是交集。
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計(jì)算
kb <- mycounts$Length / 1000
write.csv(kb, file="kb.csv")
##篩選出樣本的表達(dá)數(shù)據(jù)
# 這一步驟是篩選出所有樣本的表達(dá)量數(shù)據(jù)掺栅，第一列是length烙肺，所以不要。
countdata <- mycounts[,2:7]
##計(jì)算每個(gè)值的rpk
rpk <- countdata / kb
##計(jì)算TPM
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)
write.csv(tpm, file="mRNA_tpm.csv")
##計(jì)算FPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
head(fpkm)
write.csv(fpkm, file="mRNA_fpkm.csv")
- 主成分分析并作圖氧卧。代碼如下
tpm <- read.csv("E:/造血干/鐘/mRNA_tpm.csv")
#使用prcomp函數(shù)進(jìn)行PCA分析
# 數(shù)據(jù)轉(zhuǎn)置桃笙，轉(zhuǎn)換成行為樣本，列為基因的矩陣
tpm <- t(tpm)
data.pca <- prcomp(tpm)
# 查看PCA的分析結(jié)果
summary(data.pca)
# 繪制主成分的碎石圖
screeplot(data.pca, npcs = 10, type = "lines")
#查看主成分的結(jié)果
pca.scores <- as.data.frame(data.pca$x)
head(pca.scores)
####主成分分析可視化
library(ggpubr)
library(ggthemes)
#BiocManager::install("gmodels")
library(gmodels)
library(export)
setwd("E:/原始數(shù)據(jù)")
#導(dǎo)入計(jì)算完成的tpm數(shù)據(jù)
mytpm <- read.csv("mRNA_tpm.csv", header=TRUE)
##修改列名
rownames(mytpm) <-mytpm[,1]
mytpm <-mytpm[,-1]
#用gmodels里面自帶的fast.prcomp對(duì)PC進(jìn)行計(jì)算沙绝。
pca.info <- fast.prcomp(mytpm)
#查看PCA的結(jié)果搏明。
head(pca.info)
summary(pca.info) #能看到每個(gè)PC具體的解釋比例
head(pca.info$rotation) #特征向量，回歸系數(shù)
head(pca.info$sdev) #特征值的開(kāi)方
head(pca.info$x) #樣本得分score
#將PCA的結(jié)果和數(shù)據(jù)的表型構(gòu)建為一個(gè)數(shù)據(jù)矩陣闪檬。
##這個(gè)里面的Type=c(rep("control",50),rep("case",374))和列表里面樣本的名字是對(duì)應(yīng)的熏瞄，根據(jù)自己數(shù)據(jù)的情況靈活選擇。
pca.data <- data.frame(sample = rownames(pca.info$rotation),
Type=c(rep("D4",2),rep("D8",2),rep("D12",2)), pca.info$rotation)
#繪圖谬以，不同的參數(shù)會(huì)有不同的結(jié)果。具體的參數(shù)以及含義自行百度由桌。以下是兩個(gè)例子为黎。
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")
##我最終畫(huà)圖用的。
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)
-
結(jié)果
image.png
