新的作业
参考学习资料：
- https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103788
- RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同
- RNAseq数据，下载GEO中的FPKM文件后该怎么下游分析
按照教程拿到数据查看数据情况:
raw_counts
filtered_RPM
教程提到：
- RPM方法：
标准化了测序深度的影响，但没有考虑转录本的长度的影响。
- RPM适合于产生的read读数不受基因长度影响的测序方法，比如miRNA-seq测序，miRNA的长度一般在20-24个碱基之间。
- RPKM/FPKM方法：
标准化了基因长度的影响，
- RPKM/FPKM与RPM的区别：考虑了基因长度对read读数的影响。
- RPKM与FPKM的区别：RPKM值适用于单末端RNA-seq实验数据，FPKM适用于双末端RNA-seq测序数据。
- RPKM/FPKM适用于基因长度波动较大的测序方法，如lncRNA-seq测序，lncRNA的长度在200-100000碱基不等。
- TPM的计算方法也同RPKM/FPKM类似，首先使用式2计算每个基因的表达值，去除基因长度的影响。随后计算每个基因的表达量的百分比，最后再乘以
- TPM实际上改进了RPKM/FPKM方法在跨样品间定量的不准确性。
- TPM的使用范围与RPKM/FPKM相同。
原文作者是用的RPM标化的数据并进行了初步的过滤。
先用作者的数据进行PCA和热图查看分组是否合理。
拿到数据后先看了下初步结果，发现组内差异就比较明显
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('GSE103788_filtered_RPM_genes.tsv.gz',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表达矩阵
b=read.table('GSE103788_raw_counts_genes.tsv.gz',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表达矩阵
rownames(a)=a[,1]
a <- a[,-1]
group_list=colnames(a)
exprSet <- a
boxplot(exprSet,outline=FALSE, notch=T, las=2)
boxplot
回去看看这12个样本的简介：
RNA-seq from purified hepatocytes (wild type and peritumoral) and from livers exposed to toxic injury (24 and 48h after TAA).
意思是说作者将肝脏暴露在毒物TAA的作用下24，48小时后，分别纯化收取野生型的肝脏细胞和癌旁的肝脏细胞然后进行的RNA-seq。 由上图可以看到肿瘤组3个样本异质性比较大，对照组也不太整齐，TAA暴露组24，48小时干预后的结果比对照有所增加，其中48小时干预组增加更明显。
下面是化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图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
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')
# 每次都要检测数据
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个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
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) #把ac的行名给到n的列名，即对每一个探针标记上分组信息
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,filename = 'heatmap_top1000_sd.png')
pheatmap
虽然从PCA图上看分组没有太明显的异常，但是从热图上还是可以看出是可以分开成2个矩阵进行比较的结果。还是挺有意思的，一般来说干预样品需要4-6个重复，这个实验只设计了2个重复，但是重复性还都挺好的，也是蛮厉害了。
