2021-05-20 获取数据、熟悉矩阵

作者: 学习生信的小兔子 | 来源:发表于2021-05-20 20:08 被阅读0次

2021-05-20 获取数据、熟悉矩阵
R语言-05根据数据值获取位置
OpenCV 获取像素值的几个方法
R语言操作——TCGA数据处理
GEO下载与分析
2021-05-20
加权基因共表达网络分析（WGCNA）实例
多项式回归
Octave快速入门（2）——数据加载、存储、修改
Andrew Ng机器学习算法入门:(矩阵和向量)

视频地址：https://www.bilibili.com/video/BV1dt411Y7nn?p=10
参考：公众号：单细胞天地

熟悉两个表达矩阵

Counts数据

rm(list = ls())  ## 魔幻操作，一键清空~
Sys.setenv(R_MAX_NUM_DLLS=999) ##Sys.setenv修改环境设置，R的namespace是有上限的，如果导入包时超过这个上次就会报错,R_MAX_NUM_DLLS可以修改这个上限
options(stringsAsFactors = F)

a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',
             header = T ,sep = '\t')  ##把表达矩阵文件载入R，header=T :保留文件头部信息，seq='\t'以tap为分隔符
a[1:6,1:4]
## 读取RNA-seq的 counts 定量结果，表达矩阵需要进行简单的过滤
x=a[1,]
sum(x>1)
table(x>1)
dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),] 
dim(dat)
# 12198   768

ncol()返回矩阵的列数值；floor()四舍五入取整数；
function() 定义一个函数；sum()求和
上面的apply()指令代表对矩阵a进行行计算，判断每行表达量>1的样本总个数，并筛选出细胞表达量合格的基因（行）
第一个参数是指要参与计算的矩阵——a
第二个参数是指按行计算还是按列计算，1——表示按行计算，2——按列计算；
第三个参数是指具体的运算参数,定义一个函数x（即表达量为x）
对每行中x>1的列（即样本数）求和，即得出的是每行中表达量大于1的样本数，
然后再筛选出大于floor(ncol(a)/50)的行，这样的行（基因）的细胞表达量才算合格
因为2%的细胞有表达量，所以对于768个细胞样本，每个行(基因)在细胞中的表达至少要有15.36（约等于15）个样本表达才算合格
2 % 的细胞有表达量

  dat=log2(edgeR::cpm(dat)+1) 
##归一化的一种选择，这里是CPM(count-per-million，每百万碱基中每个转录本的count值)
  ###CPM只对read count相对总reads数做了数量的均一化，去除文库大小差异。

未归一化

归一化

dist函数

dist函数计算行与行（样本）之间的距离
cor函数计算列与列（样本）之间的相关性
scale函数默认对每一列（样本）内部归一化
计算dist之前，最好是对每一个样本（列）进行scale一下

聚类

#层次聚类，因为近 800细胞，非常耗时。
  hc=hclust(dist(t(dat))) ##样本间层次聚类
# 原始表达矩阵转置后，细胞在行，所以计算的是细胞与细胞之间的距离。
## statquest 有详细讲解背后的统计学原理。
  class(hc)
  ?plot.hclust
 ## 查看说明书。
  plot(hc,labels = FALSE)
#被聚为四类
  
 #t:矩阵转置，行转列，列转行
 #分类时常常需要估算不同样本之间的相似性(Similarity Measurement)
# 这时通常采用的方法就是计算样本间”距离”(Distance)。
#dist函数是R语言计算距离的主要函数。dist函数可以计算行与行两两间的距离。
  # 所以之前的矩阵里面行是基因，转置后行是样本，因为我们要计算样本与样本之间的距离。
  # dist()函数计算变量间距离

聚类

查看聚类结果

clus = cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数

提取批次信息

colnames(dat) #取列名
#[1] "SS2_15_0048_A3"  "SS2_15_0048_A6"  "SS2_15_0048_A5"  "SS2_15_0048_A4" 
#[5] "SS2_15_0048_A1"  "SS2_15_0048_A2"  "SS2_15_0048_A8"  "SS2_15_0048_A9" 
#[9] "SS2_15_0048_A7"  "SS2_15_0048_A10" "SS2_15_0048_A11" "SS2_15_0048_A12"
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3] #取列名，以'_'号分割，提取第三列。
#str_split()函数可以分割字符串
table(plate)
plate
#0048 0049 
#384  384

n_g = apply(a,2,function(x) sum(x>1)) #统计每个样本有表达的有多少行（基因）
# 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。
# 而且大部分单细胞转录组技术通常超过75%的基因都没办法检测到。
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
 ##(样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】)
df$all='all' #添加列，列名为"all"，没啥意思，就是后面有需要
metadata=df
save(a,dat,df,file = '../input.Rdata') #保存a,dat,df这变量到上级目录的input.Rdata
# 因为另外一个项目也需要使用这个数据集，所以保存到了上级目录。

rpkm数据

 # 每次都要检测数据
  a[1:6,1:4] #对于a矩阵取第1~6行，第1~4列
  ## 读取RNA-seq的 counts 定量结果，表达矩阵需要进行简单的过滤
  dat=a[apply(a,1, function(x) sum(x>0) > floor(ncol(a)/50)),] #筛选表达量合格的行,列数不变
 
  dat[1:4,1:4]   
  #层次聚类
  hc=hclust(dist(t(log(dat+0.1)))) ##样本间层次聚类
  # 如果是基因聚类，可以选择 wgcna 等算法 
  ## statquest 
  plot(hc,labels = F)

clus = cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数
#1   2   3   4 
#344 189 217  18 
  
#提取批次信息
colnames(dat) #取列名
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3] #取列名，以'_'号分割，提取第三列。
#str_split()函数可以分割字符串
table(plate)
  
n_g = apply(a,2,function(x) sum(x>0)) #统计每个样本有表达的有多少行（基因）
# 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。
# 而且大部分单细胞转录组技术通常超过75%的基因都没办法检测到。
  
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
  
##(样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】)
  
df$all='all' #添加列，列名为"all"，没啥意思，就是后面有需要
metadata=df
save(dat,metadata,file = '../input_rpkm.Rdata') #保存a,dat,df这变量到上级目录的input.Rdata
# 因为另外一个项目也需要使用这个数据集，所以保存到了上级目录。