单细胞转录组学习笔记-5-熟悉文献作者提供的两个表达矩阵

转载来自：刘小泽链接：https://www.jianshu.com/p/e659a2f164f7

刘小泽写于19.6.17-第二单元第三讲：熟悉文献作者提供的两个表达矩阵

笔记目的：根据生信技能树的单细胞转录组课程探索smart-seq2技术相关的分析技术
链接在：https://www.bilibili.com/video/BV1dt411Y7nn?p=6

过滤后的操作

上次得到的dat表达矩阵过滤掉低表达基因后，剩下12198个基因

看看其中的spike-in情况

grep('ERCC',rownames(dat))
[1] 12139 12140 12141 12142 12143 12144 12145 12146 12147 12148 12149 12150
[13] 12151 12152 12153 12154 12155 12156 12157 12158 12159 12160 12161 12162
[25] 12163 12164 12165 12166 12167 12168 12169 12170 12171 12172 12173 12174
[37] 12175 12176 12177 12178 12179 12180 12181 12182 12183 12184 12185 12186
[49] 12187 12188 12189 12190 12191 12192 12193 12194 12195 12196 12197 12198

一、去除细胞文库大小差异 (归一化 cpm法)

每个细胞测得数据大小不同，这样是没办法看高表达还是低表达的，必须先保证基数一样才能比较，cpm（counts per million）这个算法就是做这个事情的。

cpm是归一化的一种方法，代表每百万碱基中每个转录本的count值

注意：这个算法只是校正文库差异，而没有校正基因长度差异。要注意我们分析的目的就是：比较一个基因在不同细胞的表达量差异，而不是考虑一个样本中不同两个基因的差异，因为"没有两片相同的树叶”这个差异是正常的。但是同一个基因由于某种条件发生了改变，背后的生物学意义是更值得探索的。

用起来很简单，有现成的函数cpm() ，然后我们再用log将数据降个维度，但保持原有数据形状不变：log2(edgeR::cpm(dat)+1)

意思就是：cpm需要除以测序总reads数，而这个值作为分母会导致结果千差万别，有的特别大有的很小。为了后面可视化不受极值的影响，用log转换一下可以将数值变小，并且原来大的数值最后还是大，并不改变这个现实

那么具体这个函数做了什么事，才是真正需要了解的：

# 先看看前4行4列的数据
>   dat[1:4,1:4] 
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
0610007P14Rik              0              0             18             11
0610009B22Rik              0              0              0              0
0610009L18Rik              0              0              0              0
0610009O20Rik              0              0              1              1
# 比如先计算一下第三个样本的总测序量
> sum(dat[,3])
[1] 206831 #结果是0.2M
# 那么对于第三个样本SS2_15_0048_A5的第一个基因0610007P14Rik(结果是18)
# 计算它的cpm值：count值*1000000/总测序reads
> 18*1000000/206831
[1] 87.02757
# 和标准公式比较看看，结果完全相同
> edgeR::cpm(dat[,3])[1]
[1] 87.02757
# 因此最后就是
dat=log2(edgeR::cpm(dat)+1)

二、归一化后聚类

主要步骤为：
1、用dist函数计算细胞直接的距离，用hclus进行层次聚类，用cutree提取分群的数目（比如分成4个cluster）。注释到每个细胞所在的cluster数。
2、用strsplit提取每个细胞的板号，用apply计算每个细胞的基因表达总数。

第一步：理解dist函数

首先理解，它是计算距离用的，正如函数名称所描述的一样：distance

# 先构建一个测试矩阵
x=1:5
y=2*x
z=52:56
tmp=data.frame(x,y,z)
> tmp
  x  y  z
1 1  2 52
2 2  4 53
3 3  6 54
4 4  8 55
5 5 10 56
# 可以看到，x和y是有一定相关性的，而z和它们很难扯上关系
# 然后尝试计算x、y、z之间的距离，来验证我们的猜想
>   dist(tmp)
         1        2        3        4
2 2.449490                           
3 4.898979 2.449490                  
4 7.348469 4.898979 2.449490         
5 9.797959 7.348469 4.898979 2.449490
# 好像得到的不是我们想要的。我们想要的是x、y、z距离结果，而计算给出的是以"行"为单位的结果
# 因此，猜测dist应该是以行为输入。因此修改一下tmp，让x、y、z为行，其实也就是转置一下，转置函数用t()
>   dist(t(tmp))
           x          y
y   7.416198           
z 114.039467 107.377838

>   cor(tmp)
           x           y           z
x  1.00000000  1.00000000 -0.09844392
y  1.00000000  1.00000000 -0.09844392
z -0.09844392 -0.09844392  1.00000000

>  dist(t(scale(tmp)))
      x        y
y 0.000000         
z 4.446571 4.446571

# 可以看到dist函数计算样本直接距离和cor函数计算样本直接相关性，是完全不同的概念。虽然我都没有调它们两个函数的默认的参数。
  # 
  # 总结：
  # 
  # - dist函数计算行与行（样本）之间的距离
  # - cor函数计算列与列（样本）之间的相关性
  # - scale函数默认对每一列（样本）内部归一化
  # - 计算dist之前，最好是对每一个样本（列）进行scale一下
  #

同样的，我们这里的dat数据，是要计算细胞间的距离，也就是列与列之间的距离，使用dist(t(dat)) 计算。数据中有768个细胞，也就是要计算768个细胞核768个细胞之间的距离，计算量还是很大的。

关于dist计算距离的方法：主要有6种：”欧式euclidean”, “切比雪夫距离maximum”, “绝对值距离manhattan”, “Lance距离canberra”, “定型变量距离binary” or “明可夫斯基距离minkowski（使用时要指定p值）”。

默认使用第一种欧氏距离，它计算的是：几何空间中两点之间的距离。思想类似于勾股定理求第三条斜边的长度=》平方和再开方。

第二步：理解hclust函数

它是进行层次聚类（系谱聚类）的方法

关于hclust聚类的方法：”离差平方和法ward”, “最短距离法single”, “最长距离法complete”,”类平均法average”, “相似法mcquitty”, “中间距离法median” or “重心法centroid”。默认使用complete算法

hc=hclust(dist(t(dat))) 
# 如果要进行可视化
plot(hc,labels = FALSE) #labels这个选项的意思是不显示各个样本名称，因为样本太多，会让图看起来很乱

image

另外hclust函数还有一个亲戚：cutree，顾名思义，就是对聚类树进行修剪。我们知道聚类结果是分群的，cutree就是指定输出哪些群(结果是从大群到小群排列)

# 例如要看看分的4大群
clus = cutree(hc, 4)
group_list= as.factor(clus) #得到的这个因子型变量group_list中样本顺序和输入的顺序一致，并且属于第几类都有记录
>   table(group_list) 
group_list
  1   2   3   4 
312 300 121  35

提取批次信息

在上一步操作结果中，可以看到，样本名都是有规律的，例如：

> head(colnames(dat))
[1] "SS2_15_0048_A3" "SS2_15_0048_A6" "SS2_15_0048_A5" "SS2_15_0048_A4"
[5] "SS2_15_0048_A1" "SS2_15_0048_A2"

其中SS2_15都是一样的，Pxx也不需要管，重要的是中间的0048、0049，表示两个384孔板编号

那么如何提取？

使用strsplit函数，strsplit(x, split, fixed = FALSE) ，需要注意两点：

字符串切分后，返回的是一个列表，如果要再还原成字符串，需要用unlist()

默认情况下它是使用正则表达式的，如果不想用，可以指定fixed = TRUE

> unlist(strsplit("a.b.c", "."))
[1] "" "" "" "" ""
> unlist(strsplit("a.b.c", ".", fixed = TRUE))
[1] "a" "b" "c"

# 方法一：纯base包(思路就是：将拆分得到的list变成数据框)
options(stringsAsFactors = F)
plate=do.call(rbind.data.frame,strsplit(colnames(dat),"_"))[,3] 
# 方法二：stringr包
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3]

统计每个样本的基因表达数目

主要使用cutree剪下来的层次聚类信息、细胞板批次信息、每个样本的基因表达信息

前两个已经具备，下面进行第三个：每个样本的基因表达信息

# 还记得之前对基因进行过滤时，我们是对行进行操作
apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50))
# 这里检测每个样本中有多少基因是表达的，count值以1为标准，rpkm值可以用0为标准
n_g = apply(a,2,function(x) sum(x>1))
# 对于单细胞转录组，一般会有超过半数的基因不会表达(这个在下面构建完数据框还可以再看一下)

最后新建细胞的属性信息

可以构建数据框了：

meta=data.frame(g=group_list,plate=plate,n_g=n_g)
# 然后再添加一列，目前用不到，后续会介绍
meta$all='all'