数据标准化（中心化+均一化）

Center，就是中心化。什么是中心化呢？the variables should be shifted to be zero centered，就是说中心化以后，原始数据点的中心移到与原点重合，为什么要中心化呢？其实为了，方便后续的运算，
Scale，也就是均一化，公式为(x-mean(x))/sd(x)
就是中心化数据再除以标准差以消除不同量级的影响，例如10、100、1000、10000等，若不对不同量级的数据进行均一化处理，原来数据，方差越大，对主成分的贡献越大，这样是不对的。那是不是所有的都应该方差化呢？并不是，因为均一化之后，会增大噪音的权重。噪音就是因为人为误差导致的样品状态有一些小的差异，但是并不会影响核心。
总结：当量级相差过大，我们应使用均一化；当差距不大，我们就不必使用均一化了。

相关性

Cor 就是协方差: 两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值时另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值；如果两个变量的变化趋势相反，即其中一个变量大于自身的期望值时另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。

今天的内容是将下载好的数据读入到R里面然后进行分组

options(stringsAsFactors = F)
# 读入featureCounts的表达矩阵
data <- read.table("counts.txt",sep = "\t",header = T,row.names = 1)
data[1:3,]

# 为`data`的列进行重新命名
str <- stringr::str_split(colnames(data),"\\.",simplify=T)
colnames(data) <- str[,1]
data[1:3,]
c=unlist(lapply(rownames(data),function(x){stringr::str_split(as.character(x),"[.]",)[[1]][1]}))
dat=t(data)
colnames(dat)=c
dataa=as.data.frame(t(dat))
dataa[1:3,]#检查数据
options(stringsAsFactors = F)
#读入基因名
library(readxl)
pdata=read_xls("41576_2014_BFnrg3813_MOESM25_ESM.xls",sheet = 2,col_names=T,na = "")
table(pdata$`consensus RNA target`)
y=as.data.frame(table(pdata$`consensus RNA target`))
colnames(y)=c('group','value')
library(ggplot2)
bp<- ggplot(y, aes(x="",y=value, fill=group))+
  geom_bar(width = 1, stat = "identity")+coord_polar("y")+
  labs(x = "", y = "Total RBPs", title = "") +  ## 将标签设为空
  theme(axis.ticks = element_blank())+  ## 把左上角多出来的“刻度线”去掉
  theme(axis.text.x = element_blank())+ ## 白色的外框即是原柱状图的X轴，把X轴的刻度文字去掉即可
  theme(panel.grid=element_blank())+    ## 去掉白色圆框和中间的坐标线
  theme(panel.border=element_blank()) +  ## 去掉最外层正方形的框框
    geom_text(aes(y = rev(c(0,cumsum(rev(y$value))[-length(y$value)])+rev(y$value/2)), x = sum(y$value)/925, label = group)) 
bp
ids=pdata[,c(1,3)]
colnames(ids)=c('symbol','probe_id')
ids=ids[ids$probe_id %in% rownames(dataa),]
ids$probe_id=as.character(ids$probe_id)
dataa[1:4,1:4]   
dataa=dataa[ids$probe_id,] 
ids$median=apply(dataa,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dataa=dataa[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
dataa[1:4,1:4] 
rownames(dataa)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dataa[1:4,1:4]  #保留每个基因ID第一次出现的信息
dim(dataa)
exprSet=dataa
group_li=c(paste0(rep("AT",2)),
           paste0(rep("10A",2)), 
           paste0(rep("CA1a",2)),
           paste0(rep("DCIS",2))) #根据样本信息写分组
group_list=factor(group_li,levels = c("AT","10A","CA1a","DCIS"))
save(exprSet,group_list,file = 'phf5aRNA-seq_exprSet.Rdata')
# 大家务必注意这样的代码技巧，多次存储Rdata文件。
# 存储后的Rdata文件，很容易就加载进入R语言工作环境。
load(file = 'phf5aRNA-seq_exprSet.Rdata')
table(group_list)
# 下面代码是为了检查
if(T){
  colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet))
  group_list
  tmp=data.frame(g=group_list)
  rownames(tmp)=colnames(exprSet)
  # 组内的样本的相似性理论上应该是要高于组间的
  # 但是如果使用全部的基因的表达矩阵来计算样本之间的相关性
  # 是不能看到组内样本很好的聚集在一起。
  pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
  dim(exprSet)
  # 所以我这里初步过滤低表达量基因。
  exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
  dim(exprSet)
  
  exprSet=log(edgeR::cpm(exprSet)+1)
  dim(exprSet)
  # 再挑选top500的MAD值基因
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  M=cor(log2(exprSet+1))
  tmp=data.frame(group=group_list)
  rownames(tmp)=colnames(M)
  pheatmap::pheatmap(M,annotation_col = tmp)
  # 现在就可以看到，组内样本很好的聚集在一起
  # 组内的样本的相似性是要高于组间
  pheatmap::pheatmap(M,annotation_col = tmp,filename = 'cor.png')
  
  
  
  library(pheatmap)
  pheatmap(scale(cor(log2(exprSet+1))))
  
}
dev.off()
# 以上代码，就是为了检查R包airway及其数据集airway里面的表达矩阵和表型信息

然后进行分组

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'phf5aRNA-seq_exprSet.Rdata')
group_li=as.character(c(paste0(rep("AT",2)),
           paste0(rep("A",2)), 
           paste0(rep("CA1a",2)),
           paste0(rep("DCIS",2))))
group_list1=as.factor(group_li[c(T,T,F,F)])
exprSet1=as.data.frame(exprSet[c(T,T,F,F)])
group_list2=as.factor(group_li[c(F,F,T,T,T,T,F,F)])
exprSet2=as.data.frame(exprSet[c(F,F,T,T,T,T,F,F)])
group_list3=as.factor(group_li[c(T,T,F,F,F,F,T,T)])
exprSet3=as.data.frame(exprSet[c(T,T,F,F,F,F,T,T)])
group_list4=as.factor(group_li[c(F,F,T,T)])
exprSet4=as.data.frame(exprSet[c(F,F,T,T)])
group_list1=relevel(group_list1,ref = 'AT')
group_list2=relevel(group_list2,ref = 'A')
group_list3=relevel(group_list3,ref = 'AT')
group_list4=relevel(group_list4,ref = 'A')

晚上的时候一直在理解这个分组代码。 今天还是在走RNA -seq的流程。