很久没有跟新啦~~~忙着基金。。基金。。还是基金，然后就是文章。。文章。。还是文章。
终于有空整理之前的东西。
说绘图方法之前先要了解，转录组的PCA图的意义是什么？
为了检测样本之间的离散程度，也就是重复之间差的大不大。
1、绘图之前先解答样本重复的问题。
转录组测序一般情况下需要是三个重复。但是对于完全没有接触过测序的人来讲，就很疑惑：转录组测序为啥一定要生物学重复呢？我不要行不行？为啥大部分要3个重复？4个、5个、666个重复行不行呢？怎么样算重复？3只老鼠各测一次算3个重复，还是1只老鼠测三次算3个重复？一堆问题就很纠结，着实让人头大呀~~~
第一个问题：一定要生物学重复吗？
回答：最好要。
那什么情况下可以没有生物学重复呢？
1）科研经费太少，没钱测序。（这种情况就干脆不要测了，测1个就很鸡肋。）
2）实验证据绝对充分，然后想装点一下门面，看起来花哨一些。（那实验都做那么好了，那就多测几个嘛~要不就干脆不要测了。要不然本来能发nature，结果你“貂尾续狗”只能发个plosone，就很没必要。）
第二个问题：一定要3个重复吗？那我测2个或者4个行不行？
回答：重复的数量一定要≥3。
1）先回答设置重复的目的是什么？是为了：消除组内误差；增强结果的可靠性；检测离群样本。
1.1） 假如你给小鼠喂了一种药，不同小鼠对药的反应肯定不同，那么多个样本就可以消除小鼠之间本身的差异。
1.2）再假如，你给3只小鼠喂了药，但是其中一只就是天生免疫力极强，药物对它的影响极小。另外两只就比较相似，那后面分析的时候你就要把免疫力极强的那一只删掉，因为它的数据会对分析结果造成极大的偏差。
1.3）但是，如果你只有2只老鼠，其中一只天生免疫力极强，药物对它的影响极小。拿到测序数据一分析发现两只差别很大，那你选哪一只呢？有人说，那我肯定选免疫力正常的那一只呀。 哎呀这个问题，真的是。。。只有测序了之后你才能知道免疫力到底强不强，你给老鼠喂药之前你是不知道人家身体到底好不好的。这就是为啥不要选2个的原因。）
2）理论上来说重复越多越好，但是考虑实际情况，设置3个重复还是比较普适的方法。
具体原因参见以下文献：RNA-seq differential expression studies: more sequence or more replication?
3）动物或者植物之间样本的差异还是比较大的，所以可以多测一点，例如可以做5-10个重复之类的。土豪的话你可以测任何你觉得吉利的数字，譬如66、88、996甚至2333等等。（玩笑话哈
第三个问题：3只老鼠分别测序算重复，还是1只老鼠测3次算重复？
回答：3只老鼠各测一次。
搞清楚生物学重复和技术重复。（自行百度）
2，绘制PCA图
载入绘图的包
library(ggpubr)
library(ggthemes)
library(gmodels)
library(export)
设置运行路径并且导入你之前已经计算完的FPKM数据。
setwd("E:/8.资料/2.干细胞/2.差异表达分析")
#导入计算完成的FPKM数据
#myfpkm <- read.csv("RNAmatrix.txt", header=TRUE, row.names=1)#存在重复会报错。
myfpkm <- read.csv("RNAmatrix.txt", header=TRUE)
#设置行名
rownames(mycounts)<-mycounts[,1]
mycounts <- mycounts[,-1]
计算每个PCA的各项指数。
#用gmodels里面自带的fast.prcomp对PC进行计算。
pca.info <- fast.prcomp(myfpkm)
#查看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("control",50),rep("case",374)), pca.info$rotation)
利用ggscatter对PC进行绘图
#绘图，不同的参数会有不同的结果。具体的参数以及含义自行百度。以下是两个例子。
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())
##我最终画图用的。
ggscatter(pca.data,x="PC1", y="PC2", color="Type",
size=4, ellipse.border.remove=TRUE,
label="sample", repel=TRUE, main="PCA plot"+theme_base())
#输出到文件
graph2ppt(file="PCA_plot.ppt", width=10, aspectr=1.0)
或者可以试试3D绘制散点图
#3D：设置不同的参数可以得到不同的结果。
install.packages("scatterplot3d")
library(scatterplot3d)
scatterplot3d(pca.data$PC1, pca.data$PC2, pca.data$PC3, color=rep(c(1,2,3,4,5,6,7),each=2),pch=16)
scatterplot3d(pca.data[,3:5], color=rep(c(1,2,3,4,5,6,7),each=2),pch=16, cex.symbols=2.0)
##最终画图用的
scatterplot3d(pca.data$PC1, pca.data$PC2, pca.data$PC3,pch=16, color=c(rep("blue",50),rep("red",374)))
3D绘图的不好的地方就是，scatterplot3d里面没有参数让你展示每个点的名字，就很郁闷。
如果想实现就试试下面的方法，我也是google出来的。
#在点上面加上标签。按照pca.data[,1]的第一列的信息。
#这里貌似必须先显示图，然后将图赋值给s3d，再添加标签比较好。
s3d <- scatterplot3d(pca.data$PC1, pca.data$PC2, pca.data$PC3, pch=16, color=c(rep("blue",50),rep("red",374)))
text(s3d$xyz.convert(pca.data[2:5,3:5]), labels = pca.data[,1],
cex= 0.7, col = "black", adj=c(0,-1),font=2)
#输出结果
graph2ppt(file="PCA_plot_3D.ppt", width=10, aspectr=1.0)
网友评论