PCA

#加载R包

> library(psych)
> library(reshape2)
> library(ggplot2)
> library(factoextra)

##提前把名为ehbio_salmon.DESeq2.normalized.symbol.txt的文件放在工作目录，以制表符分割的.txt文件

> exprData <- "ehbio_salmon.DESeq2.normalized.symbol.txt"  ##这里只是定义exprData和sampleFile两个变量，文本类型的（一开始我还以为就把文件内容给赋值
> sampleFile <- "sampleFile"  ##了）

> data <- read.table(exprData,header = T,row.names = NULL,sep ="\t")  ###read.table()读入文件。
> class(data$id)  ##判断某一列的数据类型
[1] "character

> data <- read.table(exprData,header = T,row.names = "id",sep ="\t")  ##一开始不明白row.names=NULL是怎么回事，并且我试了不加这个参数时data是一样的。
Error in read.table(exprData, header = T, row.names = "id", sep = "\t") :   ##看文档知道row.names是指定哪一列为行名的，
'row.names'里不能有重复的名字   ##看了错误提示有些懵懂，之所以不能以id列为行名，是因为这列id有相同的。我下面验证了一下

> rownames(data) <- data$id
Error in `.rowNamesDF<-`(x, value = value) : 不允许有重复的'row.names'
此外: Warning message:
non-unique values when setting 'row.names': ‘ALG1L9P’, ‘ATXN7’, ‘BMS1P21’, ‘BMS1P4’, ‘CCDC39’, ‘CYB561D2’, ‘DIABLO’, ‘DNAJC9-AS1’, ‘DUXAP8’, ‘GOLGA8M’, ‘HSPA14’, ‘IPO5P1’, ‘ITFG2-AS1’, ‘LINC-PINT’, ‘LINC00484’, ‘LINC00941’, ‘LINC01238’, ‘LINC01297’, ‘LINC01422’, ‘LINC01481’, ‘MATR3’, ‘OR7E47P’, ‘PKD1P1’, ‘POLR2J3’, ‘POLR2J4’, ‘RAET1E-AS1’, ‘RF00012’, ‘RF00017’, ‘RF00019’, ‘RF02271’, ‘RGS5’, ‘RMRP’, ‘SCO2’, ‘SNHG28’, ‘SNX29P2’, ‘SPATA13’, ‘TBCE’, ‘TMSB15B’, ‘ZNF503’

# 处理重复名字，谨慎处理，先找到名字重复的原因再决定是否需要按一下方式都保留

> rownames_data <- make.names(data[,1],unique=T)  ##有效名称由字母 数字及dot和下划线组成，字符串中的空格或者-等“非法”字符都会被转为.。
> length(rownames_data)  ##上面讲第一列导出到rownames_data，遇上重复名称，按顺序加上.1 .2 ...等后缀。
[1] 27186

> data <- data[,-1,drop = F]  ##[1,]访问第一行or 1:n；[,1]访问第一列or 1:n。这里是得到去掉第一列的剩余数据，drop=F表示更改是余下数据类型不变。

> rownames(data) <- rownames_data  ##将没有重复的第一列数据保存到data中，作为行名。

****

> write.table(rownames_data,file = "rownames_data.csv",quote = F)  ##不加quote参数，字符串及行号都会加上引号。行号默认有1 2 3。。列名x
> write.table(rownames_data,file = "rownames_data_1.csv",quote = F,row.names = F)  ##
> write.table(rownames_data,file = "rownames_data_1.csv",quote = F,row.names = F,col.names = F)  ##不要引号，不要行列名

****

> data <- data[rowSums(data)>0,] #去掉和为零的行

> nrow(data)
[1] 27186

> data <- data[apply(data, 1,var)!= 0,] ##!=表示不等于；去掉方差为0 的行，这些本身没有意义，也妨碍后续运算；1表示对阵列的行进行操作

> ?rowSums  ##计算一个阵列的行的和，返回的是一个向量，长度为这个阵列的行数。
> x <- cbind(x1 = 3, x2 = c(4:1, 2:5))
> x
x1 x2
[1,] 3 4
[2,] 3 3
[3,] 3 2
[4,] 3 1
[5,] 3 2
[6,] 3 3
[7,] 3 4
[8,] 3 5
> rowSums(x)
[1] 7 6 5 4 5 6 7 8
> colSums(x)
x1 x2
24 24
> x1 <- x[rowSums(x)>5,]
> x1
x1 x2
[1,] 3 4
[2,] 3 3
[3,] 3 3
[4,] 3 4
[5,] 3 5

## apply(X, MARGIN, FUN, ...) 函数：Returns a vector or array or list of values obtained by applying a function to margins of an array or matrix.

##

X
an array, including a matrix.

MARGIN
a vector giving the subscripts which the function will be applied over. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. Where X has named dimnames, it can be a character vector selecting dimension names.

****

> mads <- apply(data,1,mad)  ##计算中值绝对偏差 (MAD, median absolute deviation)度量基因表达变化幅度

> data <- data[rev(order(mads)),] ##将数据先按索引排序-升序order()，然后在反过来，然后在按照反过来的索引顺序形成data。

> a <- c(1,4,0,9,57,32,6)
> order(a)
[1] 3 1 2 7 4 6 5  ##0在a中的索引值是3，故3在order()之后的第一位，order()返回的是索引值，以此类推

> rev(order(a))
[1] 5 6 4 7 2 1 3
> b <- a[rev(order(a))]
> b <- a[rev(order(a)),]
Error in a[rev(order(a)), ] : 量度数目不对
> b <- a[rev(order(a))]
> b
[1] 57 32 9 6 4 1 0

 

> abc <- data.frame(a = c(1:3),b = c(4:6),c = c(7:9))
> abc
  a b c
1 1 4 7
2 2 5 8
3 3 6 9
> mads <-  apply(abc,1,mad)
> mads
[1] 4.4478 4.4478 4.4478
> cbind(abc,d = c(1,1,1))
  a b c d
1 1 4 7 1
2 2 5 8 1
3 3 6 9 1
> abc
  a b c
1 1 4 7
2 2 5 8
3 3 6 9
> abc <-  cbind(abc,d = c(1,1,1)) ##给加一列
> abc
  a b c d
1 1 4 7 1
2 2 5 8 1
3 3 6 9 1
> mads <- apply(abc,1,mad)
> abc
  a b c d
1 1 4 7 1
2 2 5 8 1
3 3 6 9 1
> mads
[1] 2.2239 2.9652 3.7065
> abc  <- abc[rev(order(mads)),] ##从结果看，按行的mad的降序重新排列了数据框
> abc
  a b c d
3 3 6 9 1
2 2 5 8 1
1 1 4 7 1

##data[x:y]，访问数据框时，如果不加标点访问的是啥呢？

> data_t[1:5]
[1] 245667.7 427435.1 221687.5 371144.2 240187.2
> data_t[,1:5]
                   FN1      DCN     CEMIP    CCDC80    IGFBP5
untrt_N61311  245667.7 212953.1  40996.34 137229.15  77812.65
untrt_N052611 427435.1 360796.2 137783.10 232772.17 288609.20
untrt_N080611 221687.5 258977.3  53813.92  86258.13 210628.87
untrt_N061011 371144.2 408573.1  91066.80 212237.32 168067.42
trt_N61311    240187.2 210002.2  62301.12 136730.76  96021.74
trt_N052611   450103.2 316009.1 223111.85 226070.89 217439.21
trt_N080611   280226.2 225547.4 212724.84 124634.56 162677.38
trt_N061011   376518.2 393843.7 157919.47 236237.81 168387.36

 

> data_t <- t(data)  ##t()函数对数据框进行转置，行变列+列变行，对角线不变。

> abc_t <- t(abc)
> abc
  a b c d
3 3 6 9 1
2 2 5 8 1
1 1 4 7 1
> abc_t
  3 2 1
a 3 2 1
b 6 5 4
c 9 8 7
d 1 1 1

##文中说rows are samples and columns are variables，这句话不太懂。

> variableL <- ncol(data_t)
> variableL
[1] 27186

> if(sampleFile != ""){  ##这里的if()是何意???
+     sample <- read.table(sampleFile,header = T,row.names = 1,sep ="\t")  
+     data_t_m <- merge(data_t,sample,by = 0)
+     rownames(data_t_m) <- data_t_m$Row.names
+     data_t <- data_t_m[,-1]
+ }
> data_t[,1:5]
                   FN1      DCN     CEMIP    CCDC80    IGFBP5
trt_N052611   450103.2 316009.1 223111.85 226070.89 217439.21
trt_N061011   376518.2 393843.7 157919.47 236237.81 168387.36
trt_N080611   280226.2 225547.4 212724.84 124634.56 162677.38
trt_N61311    240187.2 210002.2  62301.12 136730.76  96021.74
untrt_N052611 427435.1 360796.2 137783.10 232772.17 288609.20
untrt_N061011 371144.2 408573.1  91066.80 212237.32 168067.42
untrt_N080611 221687.5 258977.3  53813.92  86258.13 210628.87
untrt_N61311  245667.7 212953.1  40996.34 137229.15  77812.65

## row.names = 1?

## merge()?

## rownames(data_t_m) <- data_t_m$Row.names?

row.names = 1?

> sample_2 <- read.table("new_1.txt",header = T,row.names = 2,sep ="\t")
> sample_2
                Samp
untrt1  untrt_N61311
untrt2 untrt_N052611
untrt3 untrt_N080611
untrt4 untrt_N061011
trt5      trt_N61311
trt6     trt_N052611
trt7     trt_N080611
trt8     trt_N061011
> sample_2 <- read.table("new_1.txt",header = T,row.names = NULL,sep ="\t")
> sample_2
           Samp conditions
1  untrt_N61311     untrt1
2 untrt_N052611     untrt2
3 untrt_N080611     untrt3
4 untrt_N061011     untrt4
5    trt_N61311       trt5
6   trt_N052611       trt6
7   trt_N080611       trt7
8   trt_N061011       trt8

##看来参数row.names是指定哪一列作为行名。

##所以上面的row.names = 1表示的是以第一列为行名，所以原数据框中第一行第一列的Samp就去掉了。

 

merge()?

> data_t[,1:5]
                   FN1      DCN     CEMIP    CCDC80    IGFBP5
untrt_N61311  245667.7 212953.1  40996.34 137229.15  77812.65
untrt_N052611 427435.1 360796.2 137783.10 232772.17 288609.20
untrt_N080611 221687.5 258977.3  53813.92  86258.13 210628.87
untrt_N061011 371144.2 408573.1  91066.80 212237.32 168067.42
trt_N61311    240187.2 210002.2  62301.12 136730.76  96021.74
trt_N052611   450103.2 316009.1 223111.85 226070.89 217439.21
trt_N080611   280226.2 225547.4 212724.84 124634.56 162677.38
trt_N061011   376518.2 393843.7 157919.47 236237.81 168387.36
> sample <- read.table(sampleFile,header = T,row.names = 1,sep ="\t")
> sample
              conditions
untrt_N61311       untrt
untrt_N052611      untrt
untrt_N080611      untrt
untrt_N061011      untrt
trt_N61311           trt
trt_N052611          trt
trt_N080611          trt
trt_N061011          trt ##这里可以看出data_t与sample的行名是一致的。隐隐约约觉得有个问题：这俩只是行名一致，又不是某列，这也能合并啊？？
> data_t_m <- merge(data_t,sample,by = 0)  ##merge()中参数all默认为FALSE，即合并时取交集--返回匹配的行；all = TRUE时，取按相同列取并集。
> data_t_m[,1:5]  ##data_t_m的前5列    ###为何相同的列合并之后有了行名，还叫Row.names？？
      Row.names      FN1      DCN     CEMIP    CCDC80
1   trt_N052611 450103.2 316009.1 223111.85 226070.89
2   trt_N061011 376518.2 393843.7 157919.47 236237.81
3   trt_N080611 280226.2 225547.4 212724.84 124634.56
4    trt_N61311 240187.2 210002.2  62301.12 136730.76
5 untrt_N052611 427435.1 360796.2 137783.10 232772.17
6 untrt_N061011 371144.2 408573.1  91066.80 212237.32
7 untrt_N080611 221687.5 258977.3  53813.92  86258.13
8  untrt_N61311 245667.7 212953.1  40996.34 137229.15
> ncol(sample)
[1] 1
> ncol(data_t)
[1] 27186
> ncol(data_t_m)
[1] 27188
> data_t_m[,27182:27188] ##后7列

     FSTL4   COL23A1    BEST2      IBSP      LTF     TAC1 conditions
1 1.515631 0.0000000 0.000000 0.0000000 1.722211 0.000000        trt
2 0.000000 1.7616370 1.181002 0.0000000 0.000000 0.000000        trt
3 0.000000 0.7821581 1.935371 0.7084139 2.398343 0.000000        trt
4 0.000000 3.7734772 2.274773 0.0000000 0.000000 3.956904        trt
5 5.018991 0.0000000 0.000000 4.1753179 0.000000 0.960596      untrt
6 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.000000      untrt
7 2.512400 0.0000000 0.000000 0.0000000 2.629531 0.000000      untrt
8 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.000000      untrt

####不明白的地方是data_t和sample合并之后为何行的顺序会发生变化，且变化之后的顺序又是什么顺序呢？

rownames(data_t_m) <- data_t_m$Row.names?  

> rownames(data_t_m) <- data_t_m$Row.names  ##把Row.names这一列（是第一列）内容设置为行名，但是这一列不变！
> data_t <- data_t_m[,-1]  ##去掉第一列，

 

 

> pca <- prcomp(data_t[,1:variableL],scale = T)

PCA结果展示  ##有限展示，因为刚接触，这部分图的逻辑不太懂

> fviz_eig(pca, addlabels = TRUE) ##addlabel参数为T，显示百分比。

> fviz_pca_ind(pca, repel=T) ##repel=T 自动调整文本位置
> fviz_pca_ind(pca)

 

 

后面略过，没跟着学了

参考https://mp.weixin.qq.com/s/4R14xJkQVPtaufaoXOcPIw