新的作業(yè)
參考學(xué)習(xí)資料：
- https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103788
- RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的異同
- RNAseq數(shù)據(jù)，下載GEO中的FPKM文件后該怎么下游分析
按照教程拿到數(shù)據(jù)查看數(shù)據(jù)情況:
raw_counts
filtered_RPM
教程提到：
- RPM方法：
標(biāo)準(zhǔn)化了測(cè)序深度的影響，但沒有考慮轉(zhuǎn)錄本的長(zhǎng)度的影響。
- RPM適合于產(chǎn)生的read讀數(shù)不受基因長(zhǎng)度影響的測(cè)序方法，比如miRNA-seq測(cè)序，miRNA的長(zhǎng)度一般在20-24個(gè)堿基之間。
- RPKM/FPKM方法：
標(biāo)準(zhǔn)化了基因長(zhǎng)度的影響，
- RPKM/FPKM與RPM的區(qū)別：考慮了基因長(zhǎng)度對(duì)read讀數(shù)的影響。
- RPKM與FPKM的區(qū)別：RPKM值適用于單末端RNA-seq實(shí)驗(yàn)數(shù)據(jù)，F(xiàn)PKM適用于雙末端RNA-seq測(cè)序數(shù)據(jù)。
- RPKM/FPKM適用于基因長(zhǎng)度波動(dòng)較大的測(cè)序方法，如lncRNA-seq測(cè)序，lncRNA的長(zhǎng)度在200-100000堿基不等。
- TPM的計(jì)算方法也同RPKM/FPKM類似，首先使用式2計(jì)算每個(gè)基因的表達(dá)值，去除基因長(zhǎng)度的影響。隨后計(jì)算每個(gè)基因的表達(dá)量的百分比，最后再乘以
- TPM實(shí)際上改進(jìn)了RPKM/FPKM方法在跨樣品間定量的不準(zhǔn)確性。
- TPM的使用范圍與RPKM/FPKM相同。
原文作者是用的RPM標(biāo)化的數(shù)據(jù)并進(jìn)行了初步的過濾。
先用作者的數(shù)據(jù)進(jìn)行PCA和熱圖查看分組是否合理。
拿到數(shù)據(jù)后先看了下初步結(jié)果，發(fā)現(xiàn)組內(nèi)差異就比較明顯
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('GSE103788_filtered_RPM_genes.tsv.gz',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表達(dá)矩陣
b=read.table('GSE103788_raw_counts_genes.tsv.gz',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表達(dá)矩陣
rownames(a)=a[,1]
a <- a[,-1]
group_list=colnames(a)
exprSet <- a
boxplot(exprSet,outline=FALSE, notch=T, las=2)
boxplot
回去看看這12個(gè)樣本的簡(jiǎn)介：
RNA-seq from purified hepatocytes (wild type and peritumoral) and from livers exposed to toxic injury (24 and 48h after TAA).
意思是說作者將肝臟暴露在毒物TAA的作用下24，48小時(shí)后，分別純化收取野生型的肝臟細(xì)胞和癌旁的肝臟細(xì)胞然后進(jìn)行的RNA-seq。 由上圖可以看到腫瘤組3個(gè)樣本異質(zhì)性比較大，對(duì)照組也不太整齊，TAA暴露組24，48小時(shí)干預(yù)后的結(jié)果比對(duì)照有所增加，其中48小時(shí)干預(yù)組增加更明顯。
下面是化PCA圖和熱圖：
library("FactoMineR")
library("factoextra")
library(ggplot2)
library(pheatmap)
library(edgeR)
dat=log2(edgeR::cpm(exprSet)+1)
dat[1:4,1:4]
library(stringr)
group_list=str_split(colnames(exprSet),'_',simplify = T)[,3]
table(group_list)
save(exprSet,group_list,file = 'input.Rdata')
## 下面是畫PCA的必須操作，需要看說明書。
dat=t(dat)#畫PCA圖時(shí)要求是行名時(shí)樣本名，列名時(shí)探針名，因此此時(shí)需要轉(zhuǎn)換
dat=as.data.frame(dat)#將matrix轉(zhuǎn)換為data.frame
dat=cbind(dat,group_list) #cbind橫向追加，即將分組信息追加到最后一列
library("FactoMineR")#畫主成分分析圖需要加載這兩個(gè)包
library("factoextra")
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#現(xiàn)在dat最后一列是group_list，需要重新賦值給一個(gè)dat.pca,這個(gè)矩陣是不含有分組信息的
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dat$group_list, # color by groups
#palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
ggsave('all_samples_PCA_by_groups.png')
PCA
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'input.Rdata')
# 每次都要檢測(cè)數(shù)據(jù)
dat=log2(edgeR::cpm(exprSet)+1)
dat[1:4,1:4]
group_list=group_list
table(group_list)
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列?。┤∶恳恍械姆讲?，從小到大排序，取最大的1000個(gè)
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對(duì)那些提取出來的1000個(gè)基因所在的每一行取出，組合起來為一個(gè)新的表達(dá)矩陣
n=t(scale(t(dat[cg,]))) # 'scale'可以對(duì)log-ratio數(shù)值進(jìn)行歸一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #把a(bǔ)c的行名給到n的列名，即對(duì)每一個(gè)探針標(biāo)記上分組信息
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,filename = 'heatmap_top1000_sd.png')
pheatmap
雖然從PCA圖上看分組沒有太明顯的異常，但是從熱圖上還是可以看出是可以分開成2個(gè)矩陣進(jìn)行比較的結(jié)果。還是挺有意思的，一般來說干預(yù)樣品需要4-6個(gè)重復(fù)，這個(gè)實(shí)驗(yàn)只設(shè)計(jì)了2個(gè)重復(fù)，但是重復(fù)性還都挺好的，也是蠻厲害了。
