（一）RPKM标准化
这一部分要练习的是RPKM的标准化，因为在这篇文献上传的两个矩阵，一个是原始矩阵，一个是RPKM标准化之后的矩阵。为了尽可能的还原文章的数据，我决定还是用作者上传的原始矩阵来练习。
那么首先就是要先下载原始数据。然后解压，这个我就不具体说细节了。

http://www.metagenomics.wiki/pdf/definition/rpkm-calculation
这个网站给的公式一目了然，可以看一眼。

现在看这个公式，numReads我们矩阵里就有，totalNumReads也可以算出来，那这个geneLength咋算？
基因长度，有几种算法，参考（单细胞转录组学习笔记-12-RPKM概念及计算方法）：
-选择最长的转录本
-多个转录本的均值
-非冗余外显子长度之和
-非冗余CDS之和

教程里写了两种计算基因长度的代码，我就先用一种试试，有兴趣的同学可以去看（单细胞转录组学习笔记-12-RPKM概念及计算方法）

非冗余外显子长度之和方法计算基因长度

#先加载小鼠的TxDb包
>library("TxDb.Mmusculus.UCSC.mm10.knownGene")
>txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#这个包有好几个函数，我们想采用非冗余外显子加和的方法计算基因长度，因此选择exon和genes就好了
# 取出exon和gene，分别赋予一个变量
> exon_txdb=exons(txdb)
> genes_txdb=genes(txdb)
> genes_txdb#瞅一眼
GRanges object with 24402 ranges and 1 metadata column:
            seqnames              ranges strand |     gene_id
               <Rle>           <IRanges>  <Rle> | <character>
  100009600     chr9   21062393-21073096      - |   100009600
  100009609     chr7   84935565-84964115      - |   100009609
  100009614    chr10   77711457-77712009      + |   100009614
  100009664    chr11   45808087-45841171      + |   100009664
     100012     chr4 144157557-144162663      - |      100012
        ...      ...                 ...    ... .         ...
      99889     chr3   84496093-85887516      - |       99889
      99890     chr3 110246109-110250998      - |       99890
      99899     chr3 151730922-151749960      - |       99899
      99929     chr3   65528410-65555518      + |       99929
      99982     chr4 136550540-136602723      - |       99982
  -------
  seqinfo: 66 sequences (1 circular) from mm10 genome
> exon_txdb#再瞅一眼这个
GRanges object with 440347 ranges and 1 metadata column:
                 seqnames          ranges strand |   exon_id
                    <Rle>       <IRanges>  <Rle> | <integer>
       [1]           chr1 3073253-3074322      + |         1
       [2]           chr1 3102016-3102125      + |         2
       [3]           chr1 3252757-3253236      + |         3
       [4]           chr1 3466587-3466687      + |         4
       [5]           chr1 3513405-3513553      + |         5
       ...            ...             ...    ... .       ...
  [440343] chrUn_JH584304     55112-55701      - |    440343
  [440344] chrUn_JH584304     56986-57151      - |    440344
  [440345] chrUn_JH584304     58564-58835      - |    440345
  [440346] chrUn_JH584304     58564-59690      - |    440346
  [440347] chrUn_JH584304     59592-59667      - |    440347
  -------
  seqinfo: 66 sequences (1 circular) from mm10 genome

那么上面GRanges是啥玩意？答：存储基因组坐标和相关注释信息的"容器"。
可以看出这两个GRanges里的东西不一样，其中ranges那一列是染色体的坐标序列，也就是坐标，那么我们要把这两个GRanges里重合的部分拿出来。

找重合，用findOverlaps函数。

> overlapinfo = findOverlaps(exon_txdb,genes_txdb)
> overlapinfo
Hits object with 401998 hits and 0 metadata columns:
           queryHits subjectHits
           <integer>   <integer>
       [1]        18        6852
       [2]        19        6852
       [3]        20        6852
       [4]        21        6852
       [5]        22        6852
       ...       ...         ...
  [401994]    440343       18439
  [401995]    440344       18439
  [401996]    440345       18439
  [401997]    440346       18439
  [401998]    440347       18439
  -------
  queryLength: 440347 / subjectLength: 24402

exon_txdb写在前面，于是它就作为queryHits；那么它的第18个元素和genes_txdb的6852个元素存在交集，取出来看一下：

> genes_txdb[6852]
GRanges object with 1 range and 1 metadata column:
        seqnames          ranges strand |     gene_id
           <Rle>       <IRanges>  <Rle> | <character>
  18777     chr1 4807788-4886770      + |       18777
  -------
  seqinfo: 66 sequences (1 circular) from mm10 genome
> exon_txdb[18]
GRanges object with 1 range and 1 metadata column:
      seqnames          ranges strand |   exon_id
         <Rle>       <IRanges>  <Rle> | <integer>
  [1]     chr1 4807788-4807982      + |        18
  -------
  seqinfo: 66 sequences (1 circular) from mm10 genome
#看，genes_txdb里在1号染色体的4807788-4886770有个基因，而exon_txdb里这个元素的外显子在1号染色体4807788-4807982，正好被包含在这个基因里
#下面是提取gene和exon的信息。
>t1=exon_txdb[queryHits(overlapinfo)] #queryHits(x): Equivalent to as.data.frame(x)[[1]]
>t2=genes_txdb[subjectHits(overlapinfo)] #Equivalent to as.data.frame(x)[[2]]
#看一眼t1
> t1
    seqnames   start     end width strand exon_id
1       chr1 4807788 4807982   195      +      18
2       chr1 4807823 4807982   160      +      19
3       chr1 4807830 4807982   153      +      20
4       chr1 4807892 4807982    91      +      21
5       chr1 4807896 4807982    87      +      22
......
[ reached 'max' / getOption("max.print") -- omitted 401832 rows ]#全部是40多万行

看这个t1一共是6列，但是感觉缺了点啥。。。好像是基因ID，所以要给t1加上一列基因ID。要用上面的t2。

#mcols()这个函数在这里的意思是：提取t2这个dataframe里面，第一行包含的metadata的内容。
> t1$geneid=mcols(t2)[,1]
> t1
    seqnames   start     end width strand exon_id geneid
1       chr1 4807788 4807982   195      +      18  18777
2       chr1 4807823 4807982   160      +      19  18777
3       chr1 4807830 4807982   153      +      20  18777
4       chr1 4807892 4807982    91      +      21  18777
5       chr1 4807896 4807982    87      +      22  18777
......
#这回再看t1，多了一列geneid。

我这个t1里显示，18777对应着28个外显子，如何求这个18777的长度呢？可以看到，第3行和第4行有重叠的，所以用sum-overlap的方法。下面有点复杂，一步一步来：
（1）把所有对应到同一个基因的外显子都放到一块去：
利用split(x,f)，其中x是向量或数据框，f是分组的因子。
split(t1,as.factor(t1$geneid))
（2）对列表的每个元素取start到end的全部数值：
apply(t1,1,function(y){y[2]:y[3]})
（3）对列表去重、求长度（先去overlap，再求总长度）：
length(unique(unlist(tmp)))

那么把这3步放一起，就是一个循环：

>g_l = lapply(split(t1,t1$geneid),function(x){
  tmp=apply(x,1,function(y){
    y[2]:y[3]
  })
  length(unique(unlist(tmp)))
})
> g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))
> head(g_l)
    gene_id length
1 100009600   4819
2 100009609  10403
3 100009614    553
4 100009664   1643
5    100012   1865
6    100017   7010

现在把gene_id注释上具体的基因名字：

>library(org.Mm.eg.db)
>s2g=toTable(org.Mm.egSYMBOL)
>g_l=merge(g_l,s2g,by='gene_id')
#再瞅一眼
> head(g_l)
    gene_id length        symbol
1 100009600   4819         Zglp1
2 100009609  10403       Vmn2r65
3 100009614    553       Gm10024
4 100009664   1643 F630206G17Rik
5    100012   1865          Oog3
6    100017   7010       Ldlrap1

到此，基因长度我们就计算完成了～

RPKM标准化

#载入原始表达矩阵（这是作者上传的）
> rdata<- read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
> rdata[1:6,1:6]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2
0610005C13Rik              0              0              0              1              0              0
0610007P14Rik              0              0             18             11             17              0
0610009B22Rik              0              0              0              0              8              0
0610009L18Rik              0              0              0              0              0              0
0610009O20Rik              0              0              1              1             59             28
0610010B08Rik              0              0              0              0              0              0
#上面得到的g_l和原始表达矩阵的行名取交集
> ng=intersect(rownames(rdata),g_l$symbol)
> exprSet=rdata[ng,]
> lengths=g_l[match(ng,g_l$symbol),2]
> head(lengths)
[1] 3583  998  619 5343 2990 1445
> head(rownames(exprSet))
[1] "0610005C13Rik" "0610009B22Rik" "0610009L18Rik" "0610010F05Rik" "0610010K14Rik" "0610012G03Rik"
> exprSet[1:6,1:6]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2
0610005C13Rik              0              0              0              1              0              0
0610009B22Rik              0              0              0              0              8              0
0610009L18Rik              0              0              0              0              0              0
0610010F05Rik              0              0              0             11              0              0
0610010K14Rik              0              2              0              3              0              1
0610012G03Rik              0              0              0             15              0              9
#可以看出来，现在的这个矩阵的行名里，是原始矩阵的一部分，并不是全部的基因名了，因为我们取了交集了～
> dim(rdata)
[1] 24582   768
> dim(exprSet)
[1] 22449   768
#原始矩阵里有24582个基因，现在只有22449个基因

计算每一个文库的大小：

> total_count<- colSums(exprSet) #矩阵里每一列是一个细胞，所以是对每一列求它的cound的总和
> head(total_count) #可以看出每一个文库的大小都不一样
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2 
         96713          93427         162945         123748         266083         281996 

然后用循环来计算RPKM:

> rpkm <- t(do.call( rbind,
+                    lapply(1:length(total_count),#从第一个样本到最后一个
+                           function(i){
+                             10^9*exprSet[,i]/lengths/total_count[i] #exprSet实际上就是原始矩阵和g_l取交集后的那部分，i表示每一列的read数，lengths是基因长度，total_count是每一个样品总count数
+                           }) ))
> rpkm[1:6,1:6]
     [,1]     [,2] [,3]      [,4]     [,5]      [,6]
[1,]    0 0.000000    0  2.255355  0.00000  0.000000
[2,]    0 0.000000    0  0.000000 30.12606  0.000000
[3,]    0 0.000000    0  0.000000  0.00000  0.000000
[4,]    0 0.000000    0 16.636782  0.00000  0.000000
[5,]    0 7.159561    0  8.107965  0.00000  1.186003
[6,]    0 0.000000    0 83.885177  0.00000 22.086745

最后的最后，把计算好的这个RPKM矩阵加上行名和列名

> rownames(rpkm)=rownames(exprSet)
> colnames(rpkm)=colnames(exprSet)
> rpkm[1:6,1:6]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2
0610005C13Rik              0       0.000000              0       2.255355        0.00000       0.000000
0610009B22Rik              0       0.000000              0       0.000000       30.12606       0.000000
0610009L18Rik              0       0.000000              0       0.000000        0.00000       0.000000
0610010F05Rik              0       0.000000              0      16.636782        0.00000       0.000000
0610010K14Rik              0       7.159561              0       8.107965        0.00000       1.186003
0610012G03Rik              0       0.000000              0      83.885177        0.00000      22.086745

现在RPKM就计算好了，那怎么可以验证是不是正确的呢？
举个栗子：
上面exprSet矩阵里第一行第4列，0610005C13Rik基因在SS2_15_0048_A4里原始count数是1。A4样品的总count数是123748。基因长度是3583。
那么它的RPKM应该是：

10^9*exprSet[,i]/lengths/total_count[i] = 10^9×1/3583/123748 = 2.255355

和我们用循环算出来的一模一样。好了，我人生中第一个RPKM就算出来了。
计算好了，别忘了把这个新的RPKM矩阵保存下来，大功告成啦：

> write.table(rpkm,"allsample_rpkmNormalized.txt",sep="\t")

（二）重复文章里的图
在这篇文献里，有这么一张图（https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-07582-3/MediaObjects/41467_2018_7582_MOESM1_ESM.pdf
)：

这张图怎么看呢？先看横坐标，是RPKM的均值的log10计算，纵坐标的CV是变异系数。变异系数又称离散系数或相对偏差，这个相对偏差描述的是标准偏差与平均值之比，即：cv=sd/mean*100% 。CV的意义在哪里呢？Genes which are stably expressed across replicates/experiments as the CV is low (0.5)。图里的每个黑点代表每一个基因，红点代表spike-in。说的明白点，这张图就是用ERCC的数据做了一个技术误差检测，有点像我们熟悉的定量PCR里的标准曲线，测得基因在ERCC以下，说明我们测得基因sd值小于ERCC标准的，说明基因的技术误差也是在可接受范围之内。

要重复这张图，需要用到的数据是spike-in的数据，而我们在上面一部分RPKM计算的时候由于对基因进行annotation的时候没有加入spike-in的基因，所以最后实际上计算出来的RPKM矩阵是没有sike-in的。那么如何把spike-in的RPKM也计算出来并且和我们内源基因合并呢？

#先用作者上传的原始count矩阵提取ERCC的信息
> myrdata<- read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
> grep('ERCC',rownames(myrdata))
 [1] 24491 24492 24493 24494 24495 24496 24497 24498 24499 24500 24501 24502 24503 24504 24505 24506 24507 24508 24509
[20] 24510 24511 24512 24513 24514 24515 24516 24517 24518 24519 24520 24521 24522 24523 24524 24525 24526 24527 24528
[39] 24529 24530 24531 24532 24533 24534 24535 24536 24537 24538 24539 24540 24541 24542 24543 24544 24545 24546 24547
[58] 24548 24549 24550 24551 24552 24553 24554 24555 24556 24557 24558 24559 24560 24561 24562 24563 24564 24565 24566
[77] 24567 24568 24569 24570 24571 24572 24573 24574 24575 24576 24577 24578 24579 24580 24581 24582
#把含有ERCC的行单独存入一个新变量
> myrdata_ERCC<- subset(myrdata[24491:24582,])
#还记得计算RPKM需要什么吗？我们需要知道每一个样品的总read数
> total_count_ERCC<- colSums(myrdata_ERCC)
> head(total_count_ERCC)
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2 
         28010          34555          39081          12924          27404          21635 
#对于ERCC的基因长度，文章里没有提到具体的protocol，也没有给相关的信息，于是我为了练手，就暂时用ERCC92.gtf文件里的长度那一栏作为spike-in的基因长度，首先导入spike-in的gtf文件
> lengths_ERCC<- read.table("ERCC92.gtf",sep="\t")
#看一眼，我用的是第5列数据作为基因长度
> head(lengths_ERCC)
          V1   V2   V3 V4   V5 V6 V7 V8                                            V9
1 ERCC-00002 ERCC exon  1 1061  0  +  . gene_name ERCC-00002; transcript_id DQ459430;
2 ERCC-00003 ERCC exon  1 1023  0  +  . gene_name ERCC-00003; transcript_id DQ516784;
3 ERCC-00004 ERCC exon  1  523  0  +  . gene_name ERCC-00004; transcript_id DQ516752;
4 ERCC-00009 ERCC exon  1  984  0  +  . gene_name ERCC-00009; transcript_id DQ668364;
5 ERCC-00012 ERCC exon  1  994  0  +  . gene_name ERCC-00012; transcript_id DQ883670;
6 ERCC-00013 ERCC exon  1  808  0  +  . gene_name ERCC-00013; transcript_id EF011062;
> lengths_ERCC=lengths_ERCC[,5]
> head(lengths_ERCC)
[1] 1061 1023  523  984  994  808
> rpkm_ERCC <- t(do.call( rbind,
+                     lapply(1:length(total_count_ERCC),
+                     function(i){
+                     10^9*myrdata_ERCC[,i]/lengths_ERCC/total_count_ERCC[i] 
+                     }) ))

到这里，我们就把ERCC部分的RPKM计算出来了。
把这个矩阵加上行名和列名：

> colnames(rpkm_ERCC)=colnames(myrdata_ERCC)
> rownames(rpkm_ERCC)=rownames(myrdata_ERCC)
> rpkm_ERCC[1:6,1:6]
           SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2
ERCC-00002      121169.87      128522.45      103967.35      124559.12     133651.382     116707.947
ERCC-00003       13994.44       16266.02       14382.24       16791.15      21259.677      12470.290
ERCC-00004       72427.01       58321.41       76176.54       56663.07      68028.198      73971.916
ERCC-00009       30077.82       14852.02       36743.57       20130.19       8121.478       5965.567
ERCC-00012           0.00           0.00           0.00           0.00          0.000          0.000
ERCC-00013           0.00           0.00           0.00           0.00          0.000          0.000

下面就是把这个矩阵和我们之前得到的内源基因的RPKM矩阵合并在一起～

> sample_plus_ERCC_rpkm<- rbind(sample_rpkm,rpkm_ERCC)
> write.table(sample_plus_ERCC_rpkm,"sample_with_ERCC_rpkmNormalized.txt",sep="\t")

记得把合并之后的矩阵保存下来。

接下来就可以对新的矩阵进行变异系数相关性分析了：
因为cv=sd/mean*100%，所以我们需要先计算sd和mean的值。

#计算mean、sd值
>myanalysis<- read.table("sample_with_ERCC_rpkmNormalized.txt",header=T,sep="\t")
> mean_per_gene <- apply(myanalysis, 1, mean, na.rm = TRUE)#对表达矩阵每行求均值
> sd_per_gene <- apply(myanalysis, 1, sd, na.rm = TRUE)#对表达矩阵每行求标准差

构建数据框，计算cv值：

> cv_per_gene <- data.frame(mean = mean_per_gene,
+                           sd = sd_per_gene,
+                           cv = sd_per_gene/mean_per_gene)
#给个行名
> rownames(cv_per_gene) <- rownames(myanalysis)
> head(cv_per_gene)
                    mean         sd        cv
0610005C13Rik  0.1487278   3.393665 22.817956
0610009B22Rik 30.6373867 141.386313  4.614829
0610009L18Rik  4.5066664  23.589397  5.234334
0610010F05Rik  7.9219033  21.957539  2.771750
0610010K14Rik 11.1546269  30.684081  2.750794
0610012G03Rik 37.8874591  67.395019  1.778821

再在CV的数据框中添加两列：log10cv2和log10mean。原文中横坐标是0~4，所以再加个范围。

> cv_per_gene$log10cv2=log10(cv_per_gene$cv^2)
> cv_per_gene$log10mean=log10(cv_per_gene$mean)
> cv_per_gene=cv_per_gene[cv_per_gene$log10mean < 4, ]
> cv_per_gene=cv_per_gene[cv_per_gene$log10mean > 0, ]

下面就是画图了，一长串代码：

> library(ggpubr)
> ggscatter(cv_per_gene[-grep("ERCC",rownames(cv_per_gene)),], x = 'log10mean', y = 'log10cv2',
+           color = "black", shape = 16, size = 1, # Points color, shape and size
+           xlab = 'log10(mean RPKM)', ylab = "log10(CV^2)",
+           mean.point=T,
+           cor.coeff.args = list(method = "spearman"), 
+           label.x = 3,label.sep = "\n",
+           dot.size = 2,
+           ylim=c(-0.5, 2.5),
+           xlim=c(0,4),
+           # ggp参数的意思就是：增加一个ggplot图层。一个图层是样品的基因，另一个图层是spike-in类似于标准曲线的线。
+           ggp = ggscatter(cv_per_gene[grep("ERCC",rownames(cv_per_gene)),], x = 'log10mean', y = 'log10cv2',
+                           color = "red", shape = 16, size = 1, # Points color, shape and size
+                           xlab = 'log10(mean RPKM)', ylab = "log10(CV^2)",
+                           add = "loess", #添加拟合曲线
+                           mean.point=T,
+                           add.params = list(color = "red",fill = "lightgray"),
+                           cor.coeff.args = list(method = "pearson"), 
+                           label.x = 3,label.sep = "\n",
+                           dot.size = 2,
+           )
+ )

这张图跟原文不太一样，可能是因为作者用的是经过过滤后的rpkm矩阵进行画图的，所以他的图里红线并没有接触到x轴，而我的rpkm矩阵并没有经过过滤（我猜的。。。）我又用作者上传的经过过滤的rpkm矩阵进行了同样的分析操作，得到了下面这个图：

（三）看看两个板有没有批次效应
看批次效应，课程里用的是PCA的方法，两个板的PCA如果能分开，说明这两个板的批次效应很严重，如果吻合度很高，说明没有批次效应。
之前用到的几个函数：
计算距离的dist()函数，它是按行为操作对象；归一化的scale()函数虽然是对列进行操作；现在PCA也是对行/样本进行操作，而我们现在的rpkm矩阵的样品是列，所以需要先转置。

#前面这些和之前的代码是一样的，就是提取板的信息和分组信息。只不过我用了这个新的有ERCC的rpkm的矩阵
> myanalysis<- read.table("sample_with_ERCC_rpkmNormalized.txt",header=T,sep="\t")
> options(stringsAsFactors = F)
> hc=hclust(dist(t(myanalysis))) 
> clus = cutree(hc, 4)
> group_list= as.factor(clus) 
> table(group_list)
group_list
  1   2   3   4 
718  43   6   1 
> plate=do.call(rbind.data.frame,strsplit(colnames(myanalysis),"_"))[,3]
> n_g = apply(myanalysis,2,function(x) sum(x>1))
> meta=data.frame(g=group_list,plate=plate,n_g=n_g)
> meta$all='all'
> head(meta)
               g plate  n_g all
SS2_15_0048_A3 1  0048 2883 all
SS2_15_0048_A6 1  0048 2874 all
SS2_15_0048_A5 1  0048 3437 all
SS2_15_0048_A4 1  0048 4654 all
SS2_15_0048_A1 1  0048 4588 all
SS2_15_0048_A2 1  0048 5069 all
> plate=meta$plate
> table(plate) 
plate
0048 0049 
 384  384 
#这一步开始准备画PCA
> myanalysis_bk=myanalysis #备份矩阵
> dat=myanalysis_bk
> dat=t(dat) #PCA是对行进行操作，要求每一行是一个样本，所以要先转置一下
> dat=as.data.frame(dat)
> dat=cbind(dat,plate )
> dim(dat)
[1]   768 22542
> library("FactoMineR")
> library("factoextra")
> dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
>p=fviz_pca_ind(dat.pca ,repel =T,
   geom.ind = "point", # 只显示点，不显示文字
  col.ind = dat$plate,  # 按分组上色
  palette = c("#00AFBB", "#E7B800"),
  addEllipses = TRUE, # 加环
  legend.title = "Groups")+xlim(-50,50)+ylim(-50,50) #给横纵坐标给个limit

这个图两个板的点都重合在一起，没有区分开，说明两个板之间没有批次效应

（四）探索分组
利用三种方法：logCPM、RPKM、logRPKM。
对于单细胞测序，每个细胞都是一个样本，不像bulk-RNA-seq，有对照组和处理组。所以单细胞测序不能直接进行差异分析，需要先分组，看看哪些细胞离得更近，就划分为一组，最后在组之间比较。那么如何进行分组呢？

> cpmmatrix<- read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
#用作者上传的原始count矩阵进行CPM计算，再进行log2计算
> cpmmatrix=log2(edgeR::cpm(cpmmatrix)+1)
> hc.logCPM=hclust(dist(t(cpmmatrix))) 
#对logCPM分组可视化
> plot(hc.logCPM,labels = F)
> rpkmmatrix<-read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt",header=T,sep="\t")
#对作者上传的rpkm标准化后的矩阵进行分组可视化
> hc.RPKM=hclust(dist(t(rpkmmatrix))) 
> plot(hc.RPKM,labels = F)
#对rpkm矩阵进行log2计算后进行可视化
> hc.logRPKM=hclust(dist(t(log(rpkmmatrix+0.1)))) 
> plot(hc.logRPKM,labels = F)

logCPM矩阵分组

rpkm矩阵分组

log-rpkm矩阵分组

可以看出来，logCPM和log-rpkm矩阵的分组比较像，使用table()函数看一看：

> g1 = cutree(hc.logCPM, 4);table(g1)
g1
  1   2   3   4 
287 329 134  18 
> g2 = cutree(hc.RPKM, 4)  ;table(g2)
g2
  1   2   3   4 
112 606  15  35 
> g3 = cutree(hc.logRPKM, 4)  ;table(g3)
g3
  1   2   3   4 
177 210 363  18

接下来利用tSNE方法继续判断：
画tSNE，这里用的是Rtsne这个包。Rtsne函数是对行进行操作，因此我们原来的表达矩阵需要转置后运行。

> dat_matrix <- t(rpkmmatrix) # 矩阵转置
> dat_matrix =log2(dat_matrix+0.01) 
> set.seed(42) #因为tsne函数每次运行都会绘制新的坐标，因此需要设定随机种子，保证重复性
> library("Rtsne")
> tsne_out <- Rtsne(dat_matrix,pca=FALSE,
+                   perplexity=27,theta=0.5)
> plot(tsne_out$Y,col= g1,sub = 'hc.logCPM')
>plot(tsne_out$Y,col= g3,sub = 'hc.logRPKM')
plot(tsne_out$Y,col= g2,sub = 'hc.RPKM')

logCPM-tSNE

logRPKM-tSNE

RPKM-tSNE

结果可以看出来，logCPM的群分的是最开的，其次是logRPKM，最差的是rpkm。以上是利用log2RPKM的矩阵对三种分组进行作图。

（五）差异分析
在教程里（单细胞转录组学习笔记-13-差异分析及KEGG注释简介），大神用的是logCPM的矩阵，所以这里我也用同样的方法先试一下，主要是先熟悉一下单细胞测序的差异分析流程。

#清空之前的变量
>rm(list = ls()) 
>options(stringsAsFactors = F)
#读入作者上传的原始count数据
>a<-read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
#对原始矩阵进行logCPM的计算
>logcpm_data=log2(edgeR::cpm(a)+1)
#瞅一眼计算之后的矩阵
> logcpm_data[1:6,1:6]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1 SS2_15_0048_A2
0610005C13Rik              0              0       0.000000       3.022376       0.000000       0.000000
0610007P14Rik              0              0       6.458664       6.310622       5.843701       0.000000
0610009B22Rik              0              0       0.000000       0.000000       4.784226       0.000000
0610009L18Rik              0              0       0.000000       0.000000       0.000000       0.000000
0610009O20Rik              0              0       2.543677       3.022376       7.620886       6.506269
0610010B08Rik              0              0       0.000000       0.000000       0.000000       0.000000
#聚类，和之前的代码一样
> hc=hclust(dist(t(logcpm_data)))
> plot(hc,labels = FALSE) 
> clus = cutree(hc, 4)
> group_list= as.factor(clus)
#看看4组里每一组都有多少
> table(group_list)
group_list
  1   2   3   4 
287 329 134  18 
#把分组信息存入一个变量，一会儿能用到
> g=group_list

接下来将使用RNA-seq常用的limma包进行处理。之前的bulk-RNA-seq用的都是DEseq2进行差异分析的。根据教程里说的，单细胞测序的差异分析和常规的不一样～
下面要想办法得到差异基因：

# 刚才我们得到的logCPM矩阵要先对基因进行过滤，然后赋值给一个变量
>exprSet=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
# 然后因为是smart-seq2的数据，会有ERCC spike-in
> grep('ERCC',rownames(exprSet))
[1] 12139 12140 12141 12142 12143 12144 12145 12146 12147 12148 12149 12150 12151 12152 12153 12154 12155 12156 12157 12158 12159
[22] 12160 12161 12162 12163 12164 12165 12166 12167 12168 12169 12170 12171 12172 12173 12174 12175 12176 12177 12178 12179 12180
[43] 12181 12182 12183 12184 12185 12186 12187 12188 12189 12190 12191 12192 12193 12194 12195 12196 12197 12198
#去掉spike-in
> exprSet=exprSet[!grepl('ERCC',rownames(exprSet)),]

limma需要分组信息，这就是为什么上面我要给logCPM矩阵分组的原因。

# 定义分组信息
> group_list=ifelse(g==1,'me','other');table(group_list)
group_list
   me other 
  287   481 
#调用我们要用的包
>library(edgeR)
> library(limma)

接下来就是一长串代码，按照教程说的，我们只需要ctrl+C,ctrl+V就行了～

# 定义一个函数，输入是exprSet和group_list
>do_limma_RNAseq <- function(exprSet,group_list){
  suppressMessages(library(limma))
  design <- model.matrix(~0+factor(group_list))
  colnames(design)=levels(factor(group_list))
  rownames(design)=colnames(exprSet)
  design
  
  dge <- DGEList(counts=exprSet)
  dge <- calcNormFactors(dge)
  logCPM <- cpm(dge, log=TRUE, prior.count=3)
  
  v <- voom(dge,design,plot=TRUE, normalize="quantile")
  fit <- lmFit(v, design)
  
  group_list
  cont.matrix=makeContrasts(contrasts=c('me-other'),levels = design)
  fit2=contrasts.fit(fit,cont.matrix)
  fit2=eBayes(fit2)
  
  tempOutput = topTable(fit2, coef='me-other', n=Inf)
  DEG_limma_voom = na.omit(tempOutput)
  head(DEG_limma_voom) 
  return(DEG_limma_voom) #需要什么就返回什么
}
#接下来获取第一组差异基因
>group_list=ifelse(df$g==1,'me','other')
>deg1=do_limma_RNAseq(exprSet,group_list)
#第二组差异基因
> group_list=ifelse(g==2,'me','other');table(group_list)
group_list
   me other 
  329   439
> deg2=do_limma_RNAseq(exprSet,group_list)
#第三组差异基因
> group_list=ifelse(g==3,'me','other');table(group_list)
group_list
   me other 
  134   634
> deg3=do_limma_RNAseq(exprSet,group_list)
#第四组差异基因
> group_list=ifelse(g==4,'me','other');table(group_list)
group_list
   me other 
   18   750 
> deg4=do_limma_RNAseq(exprSet,group_list)

拿到了四组差异基因，之后就是画图了。为了对应原文，每组选top18差异基因。

#选top基因，当然要先排序，这里按照logFC排个序
>deg1=deg1[order(deg1$logFC,decreasing = T),]
>deg2=deg2[order(deg2$logFC,decreasing = T),]
>deg3=deg3[order(deg3$logFC,decreasing = T),]
>deg4=deg4[order(deg4$logFC,decreasing = T),]
#然后选前18个
>cg=c(head(rownames(deg1),18),
     head(rownames(deg2),18),
     head(rownames(deg3),18),
     head(rownames(deg4),18)
)

现在万事俱备，可以准备画图了

> library(pheatmap)
#矩阵赋值一个简单的名字
> mat=exprSet[cg,]
# order()的返回值是对应“排名”的元素所在向量中的位置
> mat=mat[,order(g)]
#加入分组信息
> ac=data.frame(group=g)
> rownames(ac)=colnames(exprSet)
#归一化
> n=t(scale(t(mat)))
> n[n>1]=1
> n[n< -2]= -2
#画图
> pheatmap(n,show_rownames = T,show_colnames = F, 
+          cluster_rows = F,cluster_cols = F,
+          annotation_col = ac)

下面做一个史上最丑的火山图：

par(mfrow=c(2,2))
with(deg1,plot( logFC,-log10( adj.P.Val)))
with(deg2,plot( logFC,-log10( adj.P.Val)))
with(deg3,plot( logFC,-log10( adj.P.Val)))
with(deg4,plot( logFC,-log10( adj.P.Val)))

分别对应1-4组的差异基因

这里教程是这样说的：

看到这里火山图的形状和我们平常见到的不太一样，这是因为我们得到差异基因的方法存在问题，这里的单细胞数据不单单是原来bulk转录组的 3v3样本这样，每个细胞都是一个样本，而我们只是又将相似的细胞聚在一起当成一个大组，后来我们也是根据大组进行的差异分析（以deg1为例，就相当于312个样本 vs 剩余的456个样本）。另外还是使用的limma包（原文用的ROTS包），于是产生的差异是可以理解的

总结：教程称这个分析过程不是真正的单细胞测序流程，算是入门了解。还需要继续学习单细胞测序的R包使用。不过对于小白的我来说已经有所收获了～好好学习，天天向上～