转录组数据绘制PCA图以及生物学重复的相关问题

很久没有跟新啦~~~忙着基金。。基金。。还是基金，然后就是文章。。文章。。还是文章。
终于有空整理之前的东西。

说绘图方法之前先要了解，转录组的PCA图的意义是什么？
为了检测样本之间的离散程度，也就是重复之间差的大不大。

1、绘图之前先解答样本重复的问题。
转录组测序一般情况下需要是三个重复。但是对于完全没有接触过测序的人来讲，就很疑惑：转录组测序为啥一定要生物学重复呢？我不要行不行？为啥大部分要3个重复？4个、5个、666个重复行不行呢？怎么样算重复？3只老鼠各测一次算3个重复，还是1只老鼠测三次算3个重复？一堆问题就很纠结，着实让人头大呀~~~

第一个问题：一定要生物学重复吗？
回答：最好要。
那什么情况下可以没有生物学重复呢？
1）科研经费太少，没钱测序。（这种情况就干脆不要测了，测1个就很鸡肋。）
2）实验证据绝对充分，然后想装点一下门面，看起来花哨一些。（那实验都做那么好了，那就多测几个嘛~要不就干脆不要测了。要不然本来能发nature，结果你“貂尾续狗”只能发个plosone，就很没必要。）

第二个问题：一定要3个重复吗？那我测2个或者4个行不行？
回答：重复的数量一定要≥3。
1）先回答设置重复的目的是什么？是为了：消除组内误差；增强结果的可靠性；检测离群样本。
11） 假如你给小鼠喂了一种药，不同小鼠对药的反应肯定不同，那么多个样本就可以消除小鼠之间本身的差异。
12）再假如，你给3只小鼠喂了药，但是其中一只就是天生免疫力极强，药物对它的影响极小。另外两只就比较相似，那后面分析的时候你就要把免疫力极强的那一只删掉，因为它的数据会对分析结果造成极大的偏差。
13）但是，如果你只有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图
载入绘图的包

设置运行路径并且导入你之前已经计算完的FPKM数据。

计算每个PCA的各项指数。

利用ggscatter对PC进行绘图

或者可以试试3D绘制散点图

3D绘图的不好的地方就是，scatterplot3d里面没有参数让你展示每个点的名字，就很郁闷。
如果想实现就试试下面的方法，我也是google出来的。

RNA-seq下游分析之 PCA图_欧阳火火的博客-CSDN博客

rm(list=ls())

mydata <- readtable("C:/Users/gao/Desktop/allidtxt",header =TRUE,quote='\t',skip = 1)

smpleNames <- c('mesc_1','mesc_1','mesc_2','mesc_2','mesc_3','mesc_3','mesc_4','mesc_4','mesc_5','mesc_5','mesc_6','mesc_6','mesc_7','mesc_7','mesc_8','mesc_8')

names(mydata)[7:22] <- smpleNames

head(mydata)

countmatrix <- asmatrix(mydata[7:22])

rownames(countmatrix) <- mydata$Geneid

head(countmatrix)

save(countmatrix,file="exprRdata")

table2 <- dataframe(name = c('mesc_1','mesc_1','mesc_2','mesc_2','mesc_3','mesc_3','mesc_4','mesc_4','mesc_5','mesc_5','mesc_6','mesc_6','mesc_7','mesc_7','mesc_8','mesc_8'),condition = c('A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B'))

dds2 <- DESeqDataSetFromMatrix(countmatrix,colData=table2,design =~condition) 

dds2 <- dds[rowSums(counts(dds2)) > 1,]

head(dds2)

rld1<- rlog(dds2)

plotPCA(rld1, intgroup=c( "name","condition"))

library(ggplot2)

data1 <- plotPCA(rld1, intgroup=c("condition","name"), returnData=TRUE) 

percentVar1 <- round(100 attr(data1, "percentVar"))

p1<- ggplot(data1, aes(PC1, PC2, color=name)) +

  geom_point(size=3) +

  xlab(paste0("PC1: ",percentVar[1],"% variance")) +

  ylab(paste0("PC2: ",percentVar[2],"% variance"))

p1

主成分分析（Principal components analysis，PCA）是一种统计分析、简化数据集的方法。 它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些 不相关变量称为主成分（Principal Components）。 具体地，主成分可以看做一个线性方程，其包含一系列线性系数来指示投影方向（如图）。PCA对原始数据的正则化或预处理敏感（相对缩放）。PCA是最简单的以特征量分析多元统计分布的方法。通常情况下，这种运算可以被看作是揭露数据的内部结构，从而更好的解释数据的变量的方法。

主坐标分析（Principal Coordinates Analysis，PCoA） ，即经典多维标度（Classical multidimensional scaling），用于研究数据间的相似性。PCoA与PCA都是降低数据维度的方法， 但是差异在在于PCA是基于原始矩阵，而PCoA是基于通过原始矩阵计算出的距离矩阵。 因此，PCA是尽力保留数据中的变异让点的位置不改动，而PCoA是尽力保证原本的距离关系不发生改变，也就是使得原始数据间点的距离与投影中即结果中各点之间的距离尽可能相关（如图）。

R中有很多包都提供了PCA和PCoA,比如常用的ade4包。本文将基于该包进行PCA和PCoA的分析，数据是自带的deug，该数据提供了104个学生9门课程的成绩（见截图）和综合评定。综合评定有以下几个等级：A+,A,B,B-,C-,D。
让我们通过PCA和PCoA来看一看这样的综合评定是否合理，是否确实依据这9门课把这104个学生合理分配到不同组（每个等级一个组）。

（1）PCA分析及作图
前文已经介绍了PCA是基于原始数据，所以直接进行PCA分析即可。由于前面已经介绍过散点图的绘制方法，这里不再细讲，PCA分析完毕后我们直接作图展示结果。

整体看起来还不错，就是B-和C-的学生似乎难以区分。

（2）PCoA分析及作图

有时候PCA和PCoA的结果差不多，有时候某种方法能够把样本有效分开而另一种可能效果不佳，这些都要看样本数据的特性。
因为没有现成可供分享的微生物组数据，所以用了这个成绩的数据集。通常来说在微生物组的研究中，我们会根据物种丰度的文件对数据进行PCA或者PCoA分析，也是我们所说的beta-diveristy分析，根据PCA或者PCoA的结果看疾病组和对照组能否分开，以了解微生物组的总体变化情况。

首先参考 这篇文章 理解PCA的原理。
使用ggord包进行PCA分析和作图
R语言用列数据画PCA带椭圆图
接着参考这个 教程 熟悉DESeq2的用法。
下面直接上我的例子，用DESeq2作PCA分析：

PCA其他有关信息可以参考 这篇文章

再举个形象的栗子，假如你是一本养花工具宣传册的摄影师，你正在拍摄一个水壶。水壶是三维的，但是照片是二维的，为了更全面的把水壶展示给客户，你需要从不同角度拍几张。下图是你从四个方向拍的照片：

第一张图里水壶的背面可以看到，但是看不到前面。
第二张图是拍前面，可以看到壶嘴，这张图可以提供了第一张图缺失的信息，但是壶把看不到了。
第三张俯视图既可以看到壶嘴，也可以看到壶把，但是无法看出壶的高度。
第四张图是你打算放进目录的，水壶的高度，顶部，壶嘴和壶把都清晰可见。
PCA的设计理念与此类似，它可以将高维数据集映射到低维空间的同时，尽可能的保留更多变量。

使用 R 语言能做出像 SIMCA-P 一样的 PCA 图吗？
答案是肯定的，使用 R 语言不仅能做出像 SIMCA-P 一样的 PCA 图，还能做出比 SIMCA-P 更好看的图，而且好看的上限仅取决于个人审美风格。

主成分分析图 = 散点图 + 置信椭圆 ，散点的横纵坐标对应 PCA 的第一主成分、第二主成分。

接下来想给散点加上分类颜色：

颜色是加上了，但是椭圆咋变成了 3 个？

原来是 stat_ellipse 函数默认对每个类别的数据计算自己的置信区间。如何对多类样本只计算一个置信区间呢？查看 stat_ellipse 的帮助文档：

原来是 stat_ellipse 函数默认会继承 ggplot 中的 aes 设置，如果希望 stat_ellipse 使用自己的 aes 设置，需要将参数 inheritaes 设置为 FALSE。

接下来对样式进行微调：为不同类别样本自定义着色，添加 x 轴、y 轴标题，添加 title：

将作图结果和 SIMCA-P 对比，散点、椭圆基本完全一致，只是比它更顺眼一些罢了~

欢迎留言讨论，如果本文有帮助到你，点个赞就更好啦！

[1] Master Machine Learning With scikit-learn

[1] R 数据可视化：水平渐变色柱状图
[2] R 数据可视化：双坐标系柱线图
[3] R 数据可视化：BoxPlot
[4] R 数据可视化：环形柱状图

1、打开Origin软件和做好的图。
2、在左边的工具栏找到图形按钮，鼠标单击长按，在d出的菜单中选择Circletool，圆形图案按钮。
3、按住Ctrl键，鼠标光标在图上单击，不要松开，然后拖动鼠标。就会看到圆形图案的轮廓了。
4、调整到合适的大小，松开鼠标左键，就得到了圆形图案了。

---