（一）RPKM標準化
這一部分要練習的是RPKM的標準化，因為在這篇文獻上傳的兩個矩陣谓苟，一個是原始矩陣彬祖，一個是RPKM標準化之后的矩陣。為了盡可能的還原文章的數(shù)據(jù)纽竣，我決定還是用作者上傳的原始矩陣來練習墓贿。
那么首先就是要先下載原始數(shù)據(jù)茧泪。然后解壓，這個我就不具體說細節(jié)了聋袋。

http://www.metagenomics.wiki/pdf/definition/rpkm-calculation
這個網(wǎng)站給的公式一目了然队伟，可以看一眼。

現(xiàn)在看這個公式舱馅，numReads我們矩陣里就有缰泡，totalNumReads也可以算出來，那這個geneLength咋算代嗤？
基因長度棘钞，有幾種算法，參考（單細胞轉(zhuǎn)錄組學習筆記-12-RPKM概念及計算方法）：
-選擇最長的轉(zhuǎn)錄本
-多個轉(zhuǎn)錄本的均值
-非冗余外顯子長度之和
-非冗余CDS之和

教程里寫了兩種計算基因長度的代碼干毅，我就先用一種試試宜猜，有興趣的同學可以去看（單細胞轉(zhuǎn)錄組學習筆記-12-RPKM概念及計算方法）

非冗余外顯子長度之和方法計算基因長度

#先加載小鼠的TxDb包
>library("TxDb.Mmusculus.UCSC.mm10.knownGene")
>txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#這個包有好幾個函數(shù)，我們想采用非冗余外顯子加和的方法計算基因長度硝逢，因此選擇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是啥玩意？答：存儲基因組坐標和相關(guān)注釋信息的"容器"渠鸽。
可以看出這兩個GRanges里的東西不一樣叫乌，其中ranges那一列是染色體的坐標序列，也就是坐標徽缚，那么我們要把這兩個GRanges里重合的部分拿出來憨奸。

找重合，用findOverlaps函數(shù)凿试。

> 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()這個函數(shù)在這里的意思是：提取t2這個dataframe里面礼仗，第一行包含的metadata的內(nèi)容。
> 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，多了一列g(shù)eneid元践。

我這個t1里顯示韭脊，18777對應著28個外顯子，如何求這個18777的長度呢单旁？可以看到沪羔，第3行和第4行有重疊的，所以用sum-overlap的方法象浑。下面有點復雜蔫饰，一步一步來：
（1）把所有對應到同一個基因的外顯子都放到一塊去：
利用split(x,f)，其中x是向量或數(shù)據(jù)框愉豺，f是分組的因子篓吁。
split(t1,as.factor(t1$geneid))
（2）對列表的每個元素取start到end的全部數(shù)值：
apply(t1,1,function(y){y[2]:y[3]})
（3）對列表去重、求長度（先去overlap蚪拦，再求總長度）：
length(unique(unlist(tmp)))

那么把這3步放一起杖剪，就是一個循環(huán)：

>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

現(xiàn)在把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
#可以看出來驰贷，現(xiàn)在的這個矩陣的行名里盛嘿，是原始矩陣的一部分，并不是全部的基因名了括袒，因為我們?nèi)×私患恕?> dim(rdata)
[1] 24582   768
> dim(exprSet)
[1] 22449   768
#原始矩陣里有24582個基因次兆，現(xiàn)在只有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 

然后用循環(huán)來計算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數(shù)，lengths是基因長度城须，total_count是每一個樣品總count數(shù)
+                           }) ))
> 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

現(xiàn)在RPKM就計算好了，那怎么可以驗證是不是正確的呢糕伐？
舉個栗子：
上面exprSet矩陣里第一行第4列砰琢，0610005C13Rik基因在SS2_15_0048_A4里原始count數(shù)是1。A4樣品的總count數(shù)是123748良瞧∨闫基因長度是3583。
那么它的RPKM應該是：

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

和我們用循環(huán)算出來的一模一樣褥蚯。好了挚冤，我人生中第一個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是變異系數(shù)。變異系數(shù)又稱離散系數(shù)或相對偏差肤京，這個相對偏差描述的是標準偏差與平均值之比颊艳，即：cv=sd/mean*100% 。CV的意義在哪里呢忘分？Genes which are stably expressed across replicates/experiments as the CV is low (0.5)棋枕。圖里的每個黑點代表每一個基因，紅點代表spike-in饭庞。說的明白點戒悠，這張圖就是用ERCC的數(shù)據(jù)做了一個技術(shù)誤差檢測，有點像我們熟悉的定量PCR里的標準曲線舟山，測得基因在ERCC以下绸狐，說明我們測得基因sd值小于ERCC標準的，說明基因的技術(shù)誤差也是在可接受范圍之內(nèi)累盗。

要重復這張圖寒矿，需要用到的數(shù)據(jù)是spike-in的數(shù)據(jù)，而我們在上面一部分RPKM計算的時候由于對基因進行annotation的時候沒有加入spike-in的基因若债，所以最后實際上計算出來的RPKM矩陣是沒有sike-in的符相。那么如何把spike-in的RPKM也計算出來并且和我們內(nèi)源基因合并呢？

#先用作者上傳的原始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數(shù)
> 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，也沒有給相關(guān)的信息傲须，于是我為了練手蓝牲，就暫時用ERCC92.gtf文件里的長度那一欄作為spike-in的基因長度，首先導入spike-in的gtf文件
> lengths_ERCC<- read.table("ERCC92.gtf",sep="\t")
#看一眼泰讽，我用的是第5列數(shù)據(jù)作為基因長度
> 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

下面就是把這個矩陣和我們之前得到的內(nèi)源基因的RPKM矩陣合并在一起～

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

記得把合并之后的矩陣保存下來已卸。

接下來就可以對新的矩陣進行變異系數(shù)相關(guān)性分析了：
因為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)#對表達矩陣每行求標準差

構(gòu)建數(shù)據(jù)框梦抢，計算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的數(shù)據(jù)框中添加兩列：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參數(shù)的意思就是：增加一個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,
+           )
+ )

這張圖跟原文不太一樣望几，可能是因為作者用的是經(jīng)過過濾后的rpkm矩陣進行畫圖的绩脆，所以他的圖里紅線并沒有接觸到x軸，而我的rpkm矩陣并沒有經(jīng)過過濾（我猜的橄抹。靴迫。。）我又用作者上傳的經(jīng)過過濾的rpkm矩陣進行了同樣的分析操作楼誓，得到了下面這個圖：

（三）看看兩個板有沒有批次效應
看批次效應玉锌，課程里用的是PCA的方法，兩個板的PCA如果能分開疟羹，說明這兩個板的批次效應很嚴重主守，如果吻合度很高，說明沒有批次效應榄融。
之前用到的幾個函數(shù)：
計算距離的dist()函數(shù)参淫，它是按行為操作對象；歸一化的scale()函數(shù)雖然是對列進行操作愧杯；現(xiàn)在PCA也是對行/樣本進行操作涎才，而我們現(xiàn)在的rpkm矩陣的樣品是列，所以需要先轉(zhuǎn)置力九。

#前面這些和之前的代碼是一樣的耍铜，就是提取板的信息和分組信息。只不過我用了這個新的有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是對行進行操作跌前，要求每一行是一個樣本棕兼，所以要先轉(zhuǎn)置一下
> 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, # 加環(huán)
  legend.title = "Groups")+xlim(-50,50)+ylim(-50,50) #給橫縱坐標給個limit

這個圖兩個板的點都重合在一起舒萎，沒有區(qū)分開程储，說明兩個板之間沒有批次效應

（四）探索分組
利用三種方法：logCPM、RPKM臂寝、logRPKM章鲤。
對于單細胞測序，每個細胞都是一個樣本咆贬，不像bulk-RNA-seq败徊，有對照組和處理組。所以單細胞測序不能直接進行差異分析掏缎，需要先分組皱蹦，看看哪些細胞離得更近煤杀，就劃分為一組，最后在組之間比較沪哺。那么如何進行分組呢颓遏？

> cpmmatrix<- read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
#用作者上傳的原始count矩陣進行CPM計算飘诗，再進行l(wèi)og2計算
> 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矩陣進行l(wèi)og2計算后進行可視化
> hc.logRPKM=hclust(dist(t(log(rpkmmatrix+0.1)))) 
> plot(hc.logRPKM,labels = F)

logCPM矩陣分組

rpkm矩陣分組

log-rpkm矩陣分組

可以看出來，logCPM和log-rpkm矩陣的分組比較像，使用table()函數(shù)看一看：

> 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方法繼續(xù)判斷：
畫tSNE豫尽，這里用的是Rtsne這個包炭晒。Rtsne函數(shù)是對行進行操作执庐，因此我們原來的表達矩陣需要轉(zhuǎn)置后運行贰镣。

> dat_matrix <- t(rpkmmatrix) # 矩陣轉(zhuǎn)置
> dat_matrix =log2(dat_matrix+0.01) 
> set.seed(42) #因為tsne函數(shù)每次運行都會繪制新的坐標，因此需要設(shè)定隨機種子孽惰，保證重復性
> 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

結(jié)果可以看出來晚岭，logCPM的群分的是最開的，其次是logRPKM勋功，最差的是rpkm坦报。以上是利用log2RPKM的矩陣對三種分組進行作圖。

（五）差異分析
在教程里（單細胞轉(zhuǎn)錄組學習筆記-13-差異分析及KEGG注釋簡介）狂鞋，大神用的是logCPM的矩陣燎竖，所以這里我也用同樣的方法先試一下，主要是先熟悉一下單細胞測序的差異分析流程要销。

#清空之前的變量
>rm(list = ls()) 
>options(stringsAsFactors = F)
#讀入作者上傳的原始count數(shù)據(jù)
>a<-read.table("GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt",header=T,sep="\t")
#對原始矩陣進行l(wèi)ogCPM的計算
>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進行差異分析的纤掸。根據(jù)教程里說的，單細胞測序的差異分析和常規(guī)的不一樣～
下面要想辦法得到差異基因：

# 剛才我們得到的logCPM矩陣要先對基因進行過濾浑塞，然后賦值給一個變量
>exprSet=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
# 然后因為是smart-seq2的數(shù)據(jù)借跪，會有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 
#調(diào)用我們要用的包
>library(edgeR)
> library(limma)

接下來就是一長串代碼掏愁，按照教程說的，我們只需要ctrl+C,ctrl+V就行了～

# 定義一個函數(shù)卵牍，輸入是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)
)

現(xiàn)在萬事俱備萝衩，可以準備畫圖了

> 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組的差異基因

這里教程是這樣說的：

看到這里火山圖的形狀和我們平常見到的不太一樣回挽，這是因為我們得到差異基因的方法存在問題，這里的單細胞數(shù)據(jù)不單單是原來bulk轉(zhuǎn)錄組的 3v3樣本這樣猩谊，每個細胞都是一個樣本千劈，而我們只是又將相似的細胞聚在一起當成一個大組，后來我們也是根據(jù)大組進行的差異分析（以deg1為例牌捷，就相當于312個樣本 vs 剩余的456個樣本）队塘。另外還是使用的limma包（原文用的ROTS包），于是產(chǎn)生的差異是可以理解的

總結(jié)：教程稱這個分析過程不是真正的單細胞測序流程宜鸯，算是入門了解。還需要繼續(xù)學習單細胞測序的R包使用遮怜。不過對于小白的我來說已經(jīng)有所收獲了～好好學習淋袖，天天向上～