ScRNA_xuexi
yangxinyu
2022-10-11
a=read.table('/home/data/t040334/R_project/R_scRNA_xuexi/GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',
header = T ,sep = '\t')把表达矩阵文件载入R，header=T :保留文件头部信息，seq=’以tap为分隔符 每次都要检测数据
对于a矩阵取第16行，第14列 行是基因 列是细胞
a[1:6,1:4]## SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
## 0610005C13Rik 0 0 0 1
## 0610007P14Rik 0 0 18 11
## 0610009B22Rik 0 0 0 0
## 0610009L18Rik 0 0 0 0
## 0610009O20Rik 0 0 1 1
## 0610010B08Rik 0 0 0 0读取RNA-seq的 counts 定量结果，表达矩阵需要进行简单的过滤
dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
dat[1:4,1:4]## SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
## 0610007P14Rik 0 0 18 11
## 0610009B22Rik 0 0 0 0
## 0610009L18Rik 0 0 0 0
## 0610009O20Rik 0 0 1 1筛选表达量合格的行(基因), 列（细胞）数不变 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 % 的细胞有表达量 是一种过滤标准，过滤没有绝对标准，具体如何过滤数据看自己的需要 以下为网上查到的资料： 这个对矩阵的过滤，没有统一的标准。我查询了几篇文章，几篇文章里的过滤标准都不一样。有的是按照每一个细胞里的基因数过滤的（有的人要求一个细胞里要有7000个基因表达，有的人100就够了），有的是按照每一个基因在一定的样品数里有表达，这个就不一定了，要看你的样品量和你的分析目的。我特意贴了三篇对矩阵过滤的文章，可以看出这一步完全看心情？？？比如下面这篇： 使用limma、Glimma和edgeR，RNA-seq数据分析易如反掌http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html “在任何样本中都没有足够多的序列片段的基因应该从下游分析中过滤掉。这样做的原因有好几个。 从生物学的角度来看，在任何条件下的表达水平都不具有生物学意义的基因都不值得关注，因此最好忽略。 从统计学的角度来看，去除低表达计数基因使数据中的均值 - 方差关系可以得到更精确的估计，并且还减少了在观察差异表达的下游分析中需要进行的统计检验的数量。”
还有这篇： 单细胞RNA-seq分析 https://blog.csdn.net/myy3075/article/details/86248939 1.去除在任何细胞中都不表达的基因。 2.具有少量reads/molecules，很可能是已经被破坏或捕获细胞失败，应该移除这类细胞 3.手动过滤细胞：大多数细胞检测到的基因在7000-10000之间，这对于高深度scRNAseq来说是正常的。但是这也取决于实验协议和测序深度。比如说基于droplet或其他测序深度较低的方法，通常其每个细胞检测到的基因较少。如果细胞检出率相等，则分布应近似正常，因此我们移除分布尾部的那些细胞（检测到的基因少于7000的细胞） 4.自动过滤细胞：自动异常检测来识别可能存在问题的细胞：scater包创建一个矩阵（行代表cell，列代表不同的QC度量），然后PAC提供了按照QC度量排序的cells的2D表示，然后使用来自mvoutlier包的方法检测异常值。在这里插入图片描述 5.过滤基因：移除表达水平被认为是undetectable的基因。We define a gene as detectable if at least two cells contain more than 1 transcript from the gene. If we were considering read counts rather than UMI counts a reasonable threshold is to require at least five reads in at least two cells. 然而很多时候阈值取决于测序深度，而且很重要的一点是基因过滤一定要在细胞过滤之后，因为一些基因可能仅仅被检测到只存在在低质量的细胞里。
还有一篇是这样写的： 单细胞测序（scRNA-seq）通关||数据处理必知必会https://www.jianshu.com/p/8cee8bd4ad6f 根据基因的表达量等特征，对细胞进行过滤，通常的做法就是指定一个阈值，比如要求一个细胞中检测到的基因数必须大于100，才可以进入到下游分析，如果小于这个数字，就过滤掉该细胞。需要强调的是，在设定过滤的阈值时，需要人为判断，这样的设定方式会受到主观因素的干扰，所以往往都会指定一个非常小的过滤范围，保证只过滤掉极少数的离群值点。 综上所述，矩阵过滤没有统一标准，请自行尝试。
a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),] 对于a这个数据框，逗号左边是筛选行
dat[1:4,1:4]## SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
## 0610007P14Rik 0 0 18 11
## 0610009B22Rik 0 0 0 0
## 0610009L18Rik 0 0 0 0
## 0610009O20Rik 0 0 1 1sum(dat[,3])## [1] 206831#0610007P14Rik in SS2_15_0048_A5
log2(18*1000000/sum(dat[,3])+1)## [1] 6.459884## 18 -- > 6.459884 ## SS2_15_0048_A5
dat=log2(edgeR::cpm(dat)+1) sum(dat[,3] 相当于看一下该细胞测的数据大小 计算18条read在总文库的比例，再乘以一百万，在一百万的大小的文库上中对应的read数，最后的结果需要再log2一下 归一化的一种选择，这里是CPM(count-per-million，每百万碱基中每个转录本的count值) CPM只对read count相对总reads数做了数量的均一化，去除文库大小差异。
dist 函数
x=1:10
y=2*x
z=rnorm(10)
x## [1] 1 2 3 4 5 6 7 8 9 10 y## [1] 2 4 6 8 10 12 14 16 18 20 z## [1] 0.7559805 -2.0104416 -1.0037621 -1.7108411 -0.5407896 1.0336956
## [7] -0.5151514 -0.6181044 2.4207210 -1.3824054tmp=data.frame(x,y,z)
tmp## x y z
## 1 1 2 0.7559805
## 2 2 4 -2.0104416
## 3 3 6 -1.0037621
## 4 4 8 -1.7108411
## 5 5 10 -0.5407896
## 6 6 12 1.0336956
## 7 7 14 -0.5151514
## 8 8 16 -0.6181044
## 9 9 18 2.4207210
## 10 10 20 -1.3824054dist(tmp)## 1 2 3 4 5 6 7
## 2 3.557118
## 3 4.805902 2.452224
## 4 7.147392 4.482160 2.345200
## 5 9.037788 6.867305 4.496036 2.523692
## 6 11.183789 9.448109 7.010794 5.247140 2.734777
## 7 13.476490 11.279889 8.957608 6.813932 4.472209 2.720097
## 8 15.712674 13.488462 11.186989 9.010775 6.708649 4.767436 2.238437
## 9 17.965839 16.267612 13.846555 11.919304 9.421812 6.850098 5.349705
## 10 20.237902 17.899565 15.657055 13.420427 11.211972 9.264855 6.764032
## 8 9
## 2
## 3
## 4
## 5
## 6
## 7
## 8
## 9 3.772858
## 10 4.536977 4.411776只计算了行与行之间的距离，dist函数是针对行来处理的
dist(t(tmp))## x y
## y 19.62142
## z 20.58511 39.98958转置之后计算 x-x，x-y之间的距离
cor(tmp)## x y z
## x 1.000000 1.000000 0.251452
## y 1.000000 1.000000 0.251452
## z 0.251452 0.251452 1.000000 dist(t(scale(tmp)))## x y
## y 0.000000
## z 3.670676 3.670676可以看到dist函数计算样本直接距离和cor函数计算样本直接相关性，是完全不同的概念。虽然我都没有调它们两个函数的默认的参数。 总结： - dist函数计算行与行（样本）之间的距离 - cor函数计算列与列（样本）之间的相关性 - scale函数默认对每一列（样本）内部归一化 - 计算dist之前，最好是对每一个样本（列）进行scale一下
层次聚类，因为近 800细胞，非常耗时。##样本间层次聚类
hc=hclust(dist(t(dat))) 原始表达矩阵转置后，细胞在行，所以计算的是细胞与细胞之间的距离。 statquest 有详细讲解背后的统计学原理。 思考：什么是层次聚类？为什么要用hclust这种方法进行层次聚类，背后的统计学方法以后有机会再研究吧 以下为网上搜到的一些学习文档 你知道如何聚类吗？层次聚类与聚类树： https://zhuanlan.zhihu.com/p/107631101 三种聚类方法的比较： https://www.freesion.com/article/9257733489/
class(hc)## [1] "hclust"?plot.hclustplot(hc,labels = FALSE)
t:矩阵转置，行转列，列转行
分类时常常需要估算不同样本之间的相似性(Similarity Measurement)
这时通常采用的方法就是计算样本间”距离”(Distance)。
dist函数是R语言计算距离的主要函数。dist函数可以计算行与行两两间的距离。
所以之前的矩阵里面行是基因，转置后行是样本，因为我们要计算样本与样本之间的距离。
dist()函数计算变量间距离 hclust函数用来层次聚类
得到上面那个图之后，大概分成4类，现在要提取出样本名
clus = cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数## group_list
## 1 2 3 4
## 312 300 121 35提取批次信息
#colnames(dat) #取列名，列名太多了，我这里只显示几个
colnames(dat[,1:5])## [1] "SS2_15_0048_A3" "SS2_15_0048_A6" "SS2_15_0048_A5" "SS2_15_0048_A4"
## [5] "SS2_15_0048_A1"分割：
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3] #取列名，以'_'号分割，提取第三列。
#str_split()函数可以分割字符串
table(plate)## plate
## 0048 0049
## 384 384n_g = apply(a,2,function(x) sum(x>1))
head(n_g)## SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 SS2_15_0048_A1
## 2624 2664 3319 4447 4725
## SS2_15_0048_A2
## 5263统计每个样本有表达的有多少行（基因） 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。 而且大部分单细胞转录组技术很烂，通常超过75%的基因都没办法检测到。
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
head(df)## g plate n_g
## SS2_15_0048_A3 1 0048 2624
## SS2_15_0048_A6 1 0048 2664
## SS2_15_0048_A5 2 0048 3319
## SS2_15_0048_A4 3 0048 4447
## SS2_15_0048_A1 2 0048 4725
## SS2_15_0048_A2 3 0048 5263样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】
df$all='all' #添加列，列名为"all"，没事意思，就是后面有需要hist(n_g,breaks = 30)
对细胞进行过滤：根据自己实际研究需要可以选择 删除小于多少基因数的细胞
如果需要可以：save(a,dat,df,file = ‘../input.Rdata’)
#保存a,dat,df这变量到上级目录的input.Rdata
如果是RPKM值的矩阵：
a=read.table('/home/data/t040334/R_project/R_scRNA_xuexi/GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',
header = T ,sep = '\t')每次都要检测数据
a[1:6,1:4]## SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
## 0610005C13Rik 0 0 0.00000 6.966712
## 0610007P14Rik 0 0 74.95064 69.326130
## 0610009B22Rik 0 0 0.00000 0.000000
## 0610009L18Rik 0 0 0.00000 0.000000
## 0610009O20Rik 0 0 2.19008 3.314831
## 0610010B08Rik 0 0 0.00000 0.000000dat=a[apply(a,1, function(x) sum(x>0) > floor(ncol(a)/50)),]
dim(dat)## [1] 12689 768对于RPKM值的数据 筛选选择了x＞0，因为x＞1很难满足 过滤后也是12000多个基因 后续CPM值的那些步骤就不需要了，因为RPKM值已经normalization过了
直接进行层次聚类：
hc=hclust(dist(t(log(dat+0.1)))) 要保证是样本之间的层次聚类，如果是基因间的层次聚类 那就相当于是WGCNA了
这里代码hc=hclust(dist(t(log(dat+0.1)))) 和上文 counts矩阵的 hc=hclust(dist(t(dat))) 有所差异 区别在于 加了0.1后取log做聚类 但是上文中的dat=log2(edgeR::cpm(dat)+1) 表达矩阵的数据本身经过log之后的 所以这样看是没有去别的
但是为什么要将数据取log后再进行聚类？ 这一点我没有查到答案 相关搜索： log与否会改变rpkm形式表达矩阵top的mad基因列表 https://zhuanlan.zhihu.com/p/108672519 在统计学中为什么要对变量取对数？ https://www.zhihu.com/question/22012482?from=singlemessage
https://wenku.baidu.com/view/25513c373269a45177232f60ddccda38376be17b.html
理解一点点：取log是为了去除数据的离散程度，使数据变得集中 但是RPKM值取log合适吗，这个问题没查到，继续往下走吧，甚至RPKM和FPKM本身也是存在争议，那么这个问题不是本次学习的重点，先搁置。
plot(hc,labels = F)
clus = cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数## group_list
## 1 2 3 4
## 344 189 217 18数值和刚才基于CPM进行聚类有所差异
后续代码和上面类似
#提取批次信息
#colnames(dat) #取列名，列名太多了，我这里只显示几个
colnames(dat[,1:5])## [1] "SS2_15_0048_A3" "SS2_15_0048_A6" "SS2_15_0048_A5" "SS2_15_0048_A4"
## [5] "SS2_15_0048_A1" 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>0)) #统计每个样本有表达的有多少行（基因）
# 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。
# 而且大部分单细胞转录组技术很烂，通常超过75%的基因都没办法检测到。
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
##(样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】)
df$all='all' #添加列，列名为"all"，没事意思，就是后面有需要
metadata=df