R Notebook

数据整理和mixOmics的使用

文章：Rui, J. , Chunming, Z. , Binbin, G. , Jiawen, X. , Wei, S. , & Peng, S. . (2018). Mixomics analysis of breast cancer: long non-coding rna linc01561 acts as cerna involved in the progression of breast cancer. International Journal of Biochemistry & Cell Biology, 102, S1357272518301377-. 复现Fig1.A.The SPLSDA sample plots:miRNA

library(data.table)
library(dplyr)
BrCa_miR_Origin<-data.table::fread(file="1-s2.0-S1357272518301377-mmc5.csv",stringsAsFactors = F)#数据来自文章supplement tables
BrCaID<-colnames(BrCa_miR_Origin)
TumorID<-grep(pattern = '-01',x = BrCaID,value = T)
NormalID<-grep(pattern = '-11',x = BrCaID,value = T)
as.data.frame(t(BrCa_miR_Origin %>% select(V1,TumorID)))->BrCa_miR_Tumor0
colnames(BrCa_miR_Tumor0)<-BrCa_miR_Tumor0[1,]
BrCa_miR_Tumor0[-1,]->BrCa_miR_Tumor1
BrCa_miR_Tumor1 %>% mutate(SampleID=row.names(BrCa_miR_Tumor1),Tissue="Tumor") ->BrCa_miR_Tumor2

as.data.frame(t(BrCa_miR_Origin %>% select(V1,NormalID)))->BrCa_miR_Normal0
colnames(BrCa_miR_Normal0)<-BrCa_miR_Normal0[1,]
BrCa_miR_Normal0[-1,]->BrCa_miR_Normal1
BrCa_miR_Normal1 %>% mutate(SampleID=row.names(BrCa_miR_Normal1),Tissue="Normal") ->BrCa_miR_Normal2

dplyr::bind_rows(BrCa_miR_Tumor2, BrCa_miR_Normal2) -> BrCa_miR_Table
head(BrCa_miR_Table)

BrCa_miR_Table %>% select(-c(SampleID, Tissue))->M1
as.matrix(M1)->M2
GeneInfo <- matrix(as.numeric(M2),nrow=nrow(M2), dimnames = list(rownames(M2),colnames(M2)))

BrCa_miR_Table %>% select(Tissue) ->T1
TissueInfo<-as.factor(T1$Tissue)
class(GeneInfo)

[1] "matrix" "array" 

class(TissueInfo)

[1] "factor"

library(mixOmics)
#载入mixOmics

BrCamiRNA.plsda <- plsda(GeneInfo, TissueInfo, ncomp = 10)
set.seed(2523) # for reproducibility, only when the `cpus' argument is not used
perf.plsda.BrCamiRNA <- perf(BrCamiRNA.plsda, validation = "Mfold", folds = 5, 
                  progressBar = FALSE, auc = TRUE, nrepeat = 20) 
plot(perf.plsda.BrCamiRNA, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")
#作出的图可以看到ncom=6之后曲线变化变小了，后面就用ncomp=6进行


list.keepX <- c(10,15)  #非常不明白这个keepX的参数是怎么取的
#做出来的结果不一样可能也是因为这个keepX的参数不一样
#以下tuning sPLS-DA分析的参数
tune.splsda.BrCamiRNA <- tune.splsda(GeneInfo, TissueInfo, ncomp = 6, validation = 'Mfold', folds = 5, progressBar = F, dist = 'max.dist', measure = "BER", test.keepX = list.keepX, nrepeat = 10, cpus = 4)

error <- tune.splsda.BrCamiRNA$error.rate
ncomp <- tune.splsda.BrCamiRNA$choice.ncomp$ncomp
select.keepX <- tune.splsda.BrCamiRNA$choice.keepX[1:ncomp] 

#进行sPLS-DA分析
splsda.BrCamiRNA <- splsda(GeneInfo, TissueInfo, ncomp = ncomp, keepX = select.keepX)
#绘制The SPLSDA sample plots
plotIndiv(splsda.BrCamiRNA, comp = c(1,2),
          group = TissueInfo, ind.names = FALSE, 
          ellipse = TRUE, legend = TRUE,
          title = 'Block:miRNA')