基因组的重复序列注释的个人心得

2022-07-20 第一版
2022-07-21 第二版: 群友建议,1)可以用BUSCO或者RNA-seq mapping的方式来评估重复序列屏蔽的效果。 2)用水稻人工构建的高质量TE数据库评估。

最近又折腾了下植物基因组的重复序列注释,这里做一点总结

  1. 我整理了许多文献(主要是植物)中关于TE占比的描述和对应的基因组大小。我发现,一般100-200Mb物种占比大概是30%以内(如拟南芥和A.lyrata),400Mb大概是40%(如水稻), 超过800Mb的物种占比50%以上。
  2. 我对比了RepeatModeler + RepeatMasker和EDTA + RepeatMakser 这两种分析组合的结果,前者注释的比例都会比后者高。以前我之前发现一个基因组,大小也就是200Mb左右,但是TE注释比例在45.78%,。文章描述方法中给的是RepeatModeler + RepeatProteinMask/RepeatMasker的组合。我用EDTA + RepeatMakser是40.53%。
  3. 然而在没有TEclass注释的前提下,RepeatModeler的TE库绝大部分注释的都是Unknown,因此RepeatMasker注释的结果中绝大部分的TE也是Unknown。我不确定这是不是因为我没有正确设置参数设置的原因。
  4. 基于上一条,RepeatModeler的输出结果需要注释TE才能用于RepeatMasker。我尝试用TEclass和TEsorter对这些Unknown进行注释,两者结果相差巨大,可能1000条未知序列,TEclass会注释出900条,TEsorter只能注释出90条。
  5. 我在EDTA构建的文库的基础上,分别用TEclass和TEsorter进行注释,差距会稍微小一点,但是依旧是TEclass注释的结果多。同时,对比TEclass的注释结果,和EDTA原先的注释结果,两者结果在大类上基本一致。
  6. 通过查阅TEsorter的文献,发现benchmark中的sensitivity这一项的确是低于其他软件,rice的other TE只有 16%. 但是TEsorter的precision则优于其他软件,几乎都是100%。也就是,如果你想结果更准,TEsorter是你的不二选择,但是你想要结果更多,TEsorter就未必满足你的需求了。
  7. 我发现EDTA在F.hispida会因为TIR-finder找不到结果导致无法正确构建TE文库,不过好消息,在我自己的测序物种里面,表现都挺好的。

对于我而言,以后在重复序列注释这一块,我采取的就是EDTA + RepeatMasker这个组合。

由于我对重复序列的研究比较少,以上观点仅仅是这几天跑代码的一点心得,抛转引玉,供大家参考。

「小技巧」给Seurat对象瘦瘦身

我们在使用Seurat处理单细胞数据的时候,会发现Seurat对象会不断变大,一不小心就成为内存无底洞,

例如我的一个Seurat对象就占了22.3的内存空间

1
2
3
old_size <- object.size(seu.obj)
format(old_size, units="Gb")
# 22.3 Gb

如果我中途需要关闭Rstudio,那么为了保证自己的工作连续性,我就需要将内存中的20多G的数据保存到磁盘上,并在下次分析加载会内存。这一来一回,考虑到磁盘的读写速度,耗时可能就需要10多分钟。

考虑到你可能要把数据上传到网盘或者GEO数据库,那么这20多G数据所需要花费的时间,就更加超出你的想象了。

那有没有办法给Seurat对象瘦身呢? 其实很简单,因为Seurat主要是在Scale这一步,将原本的稀疏矩阵变成了普通的矩阵,同时里面的元素都是浮点型,极其占用空间。只要我们在保存数据之前先把这个归一化的矩阵给清空,就可以让Seurat一下子瘦下来。

1
2
3
4
seu.obj@assays$RNA@scale.data <- matrix()
new_size <- object.size(seu.obj)
format(new_size, units="Gb")
# 5Gb

上面的操作,让内存占用降低到只有原来的20%左右。但是,问题来了,内存降低的代价是什么呢?代价就是,你需要对加载的数据进行scale,复原Seurat中的scale.data。

1
2
all.genes <- rownames(seu.obj)
seu.obj <- ScaleData(seu.obj, features = all.genes)

这就是计算机科学中常见思路,要么是空间换时间,要么是时间换空间。

这个小技巧除了能加速Seurat对象的保存和读写外,还有什么其他应用吗?因为scale后数据主要是给主成分分析(PCA)提供输入。后续的非线性降维(UMAP), 聚类分析(Cluster)都是基于PCA,而不是基于scale数据,因此,如果分析过程中发现内存空间吃惊,也可以先通过这个小技巧释放下内存空间。等跑完一些占内存的操作后,再恢复即可。

Arm架构的macOS安装conda

在arm架构的macOS系统安装conda实际上并没有特别需要注意的地方,我们只需要在最初下载的时候,选择conda-forge提供的miniconda,即Miniforge

Miniforge是针对conda-forge优化的conda,做了如下的预设

  • conda-forge作为默认channel
  • PyPy的可选支持,替代标准Python(CPython)
  • Mambda可选支持,替代conda
  • 多种CPU架构支持(x86_64, ppc64le, aarch64(M1))

安装过程代码如下

1
2
3
4
5
6
# 下载
cd ~/Download
wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-MacOSX-arm64.sh
# 安装
chmod +x ~/Downloads/Miniforge3-MacOSX-arm64.sh
sh ~/Downloads/Miniforge3-MacOSX-arm64.sh

运行时会有一些提示信息,输入yes或者回车就行,完成后会有如下信息

安装信息

注意我在途中标识的部分,这里有两个信息我们需要指导

  • 我们需要重启终端才能调用conda
  • conda之后会默认启动base环节,除非你用 conda config --set auto_activate_base false 取消这个行为。

不过,遗憾的是许多生信软件还没有arm64版本,例如bwa, samtools 等,但是没有关系这些软件可以通过 homebrew 进行安装。

上手了一个自然语言模型BLOOM

在2年前,OpenAI搞了一个1750亿个参数的神经网络模型, 即GPT-3(是它的前辈GPT-2, 约15亿个参数, 的100多倍)。可惜的是, GPT-3并没有开源它预训练的参数数据集,使用者只能调用他们提供的API,由于这个API必须要申请,所以我还没机会用到。不过,已经有很多公司基于GPT-3开发了一些应用,比如说GitHub就开放一个GitHub Copilot用于智能补全代码,大家感兴趣的可以去试试。

通常而言,GPT-3这种级别的神经网络模型只有OpenAI, 谷歌这种巨头科技公司能够开发和训练,毕竟1750个参数可不是一个小数目。但是,一个由约1000多个学术志愿者组成的国际志愿者不信邪,他们正在用价值700万美元的公共资金赞助的计算时间去训练一个具有1760亿个参数的自然语言模型,BLOOM。

出于兴趣爱好,我在自己的Macbook Pro上测试了这个BLOOM-1b3模型,顾名思义就是该模型有13亿参数。为什么不测试完整的1750亿个参数呢?因为光读取13亿参数到Python中,就需要占用6G左右的内存,完整的1750亿个参数,起码得是一个1T内存的服务器才行。而且完整参数版还没有完工,仍在训练中。

BLOOM是基于Transformers,而Transformers的安装要求满足如下三点

  • Python 3.6+
  • Flax 0.3.2+/ PyTorch 1.3.1+ / TensorFlow 2.3+ 任一框架
  • Rust(用于编译tokenizers)

在Mac上,我们可以通过homebrew来安装rust

1
brew install rust

在深度学习框架上,我选择了PyTorch, 因为近期支持使用M1芯片的GPU进行加速

1
2
3
4
5
6
7
8
9
10
# 安装虚拟环境
python3 -m pip install --user --upgrade pip
python3 -m pip install --user virtualenv
# 创建虚拟环境
python3 -m venv pytorch
source pytorch/bin/activate
# 安装pytorch
pip3 install --pre torch torchvision torchaudio --extra-index-url https://download.pytorch.org/whl/nightly/cpu
# 安装transformers
pip3 install transformers

之后打开python, 加载我们的模型(第一次运行时, transformers会自动帮忙下载模型)

1
2
3
4
5
6
from transformers import AutoTokenizer, AutoModelForCausalLM

# 分词
tokenizer = AutoTokenizer.from_pretrained('bigscience/bloom-1b3')
# 模型
model = AutoModelForCausalLM.from_pretrained('bigscience/bloom-1b3')

由于我对自然语义处理一窍不通,所以只能根据示例代码,使用transformres的pipeline, 构建一个文本生成工具。

1
2
3
4
5
6
7
8
9
10
11
12
13
from transformers import pipeline
generator = pipeline(task="text-generation", model=model, tokenizer=tokenizer)

# 英文输入
generator("bioinformatics is ", max_length=30, num_return_sequence=5)

# 输出结果
[{'generated_text': 'bioinformatics is a branch of computer science that deals with the analysis of data and the design of algorithms that can be used to solve problems in'}]

# 中文输入
generator("生物信息学是一门")
# 输出结果
[{'generated_text': '生物信息学是一门新兴的交叉学科,它涉及生物信息学、计算机科学、数学'}]

从结果来看,13亿参数的BLOOM模型的表现已经非常符合我的预期,远超出之前摸索Transformers时测试的默认GPT-2模型。

可惜我的能力有限,对BLOOM的探索就止步于此了,等到后续我对Transformers有更多的掌握后,再更新相关的内容吧。

参考资料

服务器上R调用png显示x11报错怎么办?

太长不读版

  • 治本的方法,服务器安装pango, 之后重新编译R语言
  • 治标的方法,在R的配置文件中增加options(bitmapType='cairo')

服务器上装完R语言之后,发现自己的PNG函数无法正常调用,始终会出现一个X11连接失败的错误。

1
2
3
4
5
> png('ab.png')
Error in .External2(C_X11, paste0("png::", filename), g$width, g$height, :
unable to start device PNG
In addition: Warning message:
In png("ab.png") : unable to open connection to X11 display ''

使用 capabilities查询之后,的确是没有启用X11。

1
2
3
4
5
6
7
8
9
> capabilities()
jpeg png tiff tcltk X11 aqua
TRUE TRUE TRUE TRUE FALSE FALSE
http/ftp sockets libxml fifo cledit iconv
TRUE TRUE FALSE TRUE TRUE TRUE
NLS Rprof profmem cairo ICU long.double
TRUE TRUE FALSE TRUE TRUE TRUE
libcurl
TRUE

回溯到我源代码的配置过程,发现配置项目包括X11的界面支持,Interfaces supported:X11, tcltk. 但是我思索了下,这个X11支持指的应该是Linux的图形界面下的R,而不是终端的R,我的直觉告诉我,这个信息应该是和我们的报错无关了。

之后,我又去查了下其他资料,谢益辉提到这可能是png函数的type参数选择出现了问题,也就是下面这段函数

1
2
if (missing(type))
type <- getOption("bitmapType")

进一步这个bitmapType选项的设置又和 grDevices:::.onLoad函数中的.Call(C_cairoProps, 2L)有关。其中.Call是R语言调用C语言编写函数的一种方式,回溯到对应的C代码(代码位置为src/library/grDevices/src/init.c)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#ifndef _WIN32
/* This really belongs with the X11 module, but it is about devices */
static SEXP cairoProps(SEXP in)
{
int which = asInteger(in);
if(which == 1)
return ScalarLogical(
#ifdef HAVE_WORKING_CAIRO
1
#else
0
#endif
);
else if(which == 2)
return ScalarLogical(
#ifdef HAVE_PANGOCAIRO
1
#else
0
#endif
);
return R_NilValue;
}
#endif

还好,我略懂一些C语言,发现这个输出结果又和变量 HAVE_PANGOCAIRO 有关。 搜索R源代码中所有和HAVE_PANGOCAIRO有关的内容,发现在configure有相关记录,对应的逻辑判断语句如下

1
2
3
4
5
if test "x${r_cv_has_pangocairo}" = xyes; then

printf "%s\n" "#define HAVE_PANGOCAIRO 1" >>confdefs.h

fi

进一步,我们在configure中查找 r_cv_has_pangocairo, 找到如下的判断语句

1
2
3
4
5
6
7
8
else $as_nop
if "${PKG_CONFIG}" --exists pangocairo; then
r_cv_has_pangocairo="yes"
else
r_cv_has_pangocairo="no"
fi

fi

那看来有可能是pangocairo这个库出现了问题。使用下列代码在我的服务器上测试,发现输出是no, 说明的确没有安装。

1
/bin/pkg-config --exists pangocairo && echo "yes" || echo "no"

我尝试着用下述代码安装pango

1
2
3
4
# CentOS
sudo yum install pango-devel pango libpng-devel
# ubuntu
sudo apt install libpango-1.0-0

对R进行重新编译安装,发现 getOption("bitmapType") 得到的结果是cairo,符合预期了。

当然,你可能也不需要那么麻烦的方法,还有一种简单的思路是在 ~/.Rprofile中加入这一行,手动指定cairo。

1
options(bitmapType='cairo')

注意: 直接用apt/yum安装的R, getOption("bitmapType") 输出直接是 cairo.

参考资料

「macOS」为什么我的时间机器的备份磁盘无法被推出?

今天我用时间机器备份我的Mac上的数据后,打算推出我的磁盘,但是遇到了如下提示

image.png

一般这种情况,我们有一个国际通用的解决方案,那就是关机。但是,我不想关机啊,我想找到到底是什么程序占用了这个磁盘,但是它不能顺利的被推出。

作为一个略懂linux的半个码农,我熟练的使用lsof 查看了磁盘的占用进程

1
2
3
$ sudo lsof  /Volumes/Time-Machine
COMMAND PID USER FD TYPE DEVICE SIZE/OFF NODE NAME
mds 334 root 27r DIR 1,29 160 2 /Volumes/Time-Machine

结果显示,有一个PID为334的mds的进程正在使用硬盘。那么问题来了,这个mds是啥呢? 使用百度加上 man mds, 我发现这是MacOS上的一个索引工具,同时我还找到一位跟我有类似经历的文章。简单说,就是因为时间机器备份需要进行文件的比对,为了实现这种快速的文件件对比,那么我们就需要建立文件索引。

强制干掉mds进程会影响我们索引的完整性,所以,最后我还是关机了。

「流畅的Python」关于函数参数,自省和闭包的一些记录

「流畅的Python」英文原版出版于2015年,2017年出版了中文版,我看到一些人推荐,于是就买了。但是遗憾的是那时候Python的道行太浅,发现每个字都认识,但是连起来就是读不懂。现在是2022年,由于疫情缘故,我终于有时间去深入学习一些Python知识。并且,经过几年的编程积累,我也终于有能力阅读这本书了。

“在Python中,函数是一等对象”。(最初看到这句话时,我脑海里想的是,编程语言还分个三六九等吗?)一等对象(first-class object, 又称第一类对象)是编程语言理论家定义的,符合下述要求的程序实体

  • 在运行时创建
  • 能赋值给变量或数据结构的元素
  • 能作为参数传给函数
  • 能作为函数的返回结果

Python的函数之所以是一等对象,就是因为它满足上述条件,也就是能够实现下面这些操作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def foo():
print("foo")

#作为函数的参数
def bar(func):
func()
print("bar")

#作为函数的返回
def zoo(func):
return func

f = foo #能赋值给变量
f()

ff = zoo(foo)
ff()

当一个函数能够接受函数作为参数,或者能够将函数作为结果进行返回,那么它就是一个高阶函数(higher-order function).因此上面代码的 bar, zoo 就是高阶函数。

我们常用的sorted函数也是高阶函数, 我们可以通过传入自定义的函数来改变默认行为,下面的案例代码使得我们可以根据value对字典的key排序。

1
2
3
d = {'a': 100, 'b': 1, 'c': 2}
# lambda是一个匿名函数, 这里是用来获取字典的值,
sorted(d, key = lambda x : d[x] )

sorted中的key后面用到了lambda匿名函数,是一种非常简单的表达式,不支持try/catch, while,for等语句。由于lambda还会带来阅读上的困难,我们会更建议定义一个有名字的函数来替代。

Python函数支持非常灵活的参数处理机制,我们可以使用 ***来展开可迭代对象,映射到单个参数中。例如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#案例代码
def test_arguments(first, *middle, **last):
print(f"first is: {first}")

for i,j in enumerate(middle):
print(f"the {i} of middle is {j}")

for k,v in last.items():
print(f"key {k} is value {v}")

# middle捕获了b,c,d, 存成一个列表, last将k=1,k2=2存成一个字典, {"k":1,"k2":2}
test_arguments("a", "b", "c", "d", k=1, k2=2)

#first is: a
#the 0 of middle is b
#the 1 of middle is c
#the 2 of middle is d
#key k is value 1
#key k2 is value 2

因为Python一切皆对象,因此函数也是对象,既然是对象,我们就可以查看其属性(函数自省)。 Python使用 dir函数可以查看函数的所有属性, 如 dir(sorted)

这里的自省是一个计算机科学术语,是在程序运行时检查对象类型的一种能力,通常也称之为’运行时类型检查‘。我一开始还以为函数能自我反省,所以面对不理解的名词(包括后面出现的闭包),要果断去查,而不是由着自己乱想,。

log化的TPM能做差异分析吗

有人提了一个问题: 我用log2TPM做差异,log2本来就是对TPM标准化的嘛,我看样品表达分布还是不是特别齐,还需要再进一步标准化吗

对于这个问题,我的回答如下:

首先,我们如果要对转录组做差异分析,那么最好的数据应该就是raw count,而不是这些转换后的数据。这是因为差异分析的本质是做统计检验,而我们应该都知道统计检验都会对数据有一些基础假设,例如,我们平常用的t-test,默认假设是数据服从正态分布。常用的DESeq2, limma(voom), edgeR之所以都要求输入数据是count,就是因为它们的统计假设是基于count。

其次,如果我们只能拿到TPM或log2TPM的数据,那在这个基础上做差异分析所能用的统计检验方法就建议是无参方法,例如wilcoxon秩和检验。尽管log2转换会让偏倚的数据更符合钟形曲线(正态分布),但你不能在基础上用t-test。 举个例子,

1
2
3
4
5
6
7
8
9
x <- c(100, 200, 300) 
y <- c(1000, 2000, 3000)

t.test(x,y) # p=0.08787
t.test(log(x), log(y)) # p=0.0071

wilcox.test(x,y) # p=0.1
wilcox.test(log(x), log(y)) # p=0.1

在我构建的两个数据中,经log转换,t-test从p>0.05变成了p<0.01, 而秩和检验没有变化。

最后,我们在转录组分析中会用到TPM归一化和log2转换的目的并不是为了差异分析,而是为了后续的可视化。比如说绘制热图,我们就不希望数据里面不要有太高和太低的值,使得整体的颜色极端化。

In order to test for differential expression, we operate on raw counts and use discrete distributions as described in the previous section on differential expression. However for other downstream analyses – e.g. for visualization or clustering – it might be useful to work with transformed versions of the count data.
Analyzing RNA-seq data with DESeq2

综上,我们用log2TPM做差异应该采用无参的方法,而无参的方法并不考虑样本的分布。

推荐阅读:

使用非负最小二乘回归进行细胞类型转移

2019年发表在Nature上的文章【The single-cell transcriptional landscape of mammalian organogenesis】在方法部分提到,使用NNLS(non-negative linear-square)回归的方法分析两个细胞图谱数据集中相关细胞类型。

这个方法,在我搜索的中文教程中都没有出现过,所以这里以两个pbmc的数据集为例进行介绍,如何复现文章的方法。

10x的细胞数据集的预处理部分不做过多介绍, 如下代码以10x官网提供的数据为例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(Seurat)
pbmc.data <- Read10X_h5("./data/10x/10k_PBMC_3p_nextgem_Chromium_X_filtered_feature_bc_matrix.h5")
seu.obj <- CreateSeuratObject(counts = pbmc.data,
project = "pbmc10k",
min.cells = 3,
min.features = 200)
library(dplyr)
seu.obj <- seu.obj %>%
Seurat::NormalizeData(verbose = FALSE) %>% #归一化
FindVariableFeatures(selection.method = "vst") %>% #筛选特征,找高变基因
ScaleData(verbose = FALSE) %>% #标准化
RunPCA(pc.genes = seu.obj@var.genes, npcs = 100, verbose = FALSE) %>% #降维,降噪
FindNeighbors( dims = 1:10) %>% #建立SNN
FindClusters( resolution = 0.5) %>% #分群
RunUMAP( dims = 1:10) # 非线性降维

相同的处理流程我们应用到Seurat 3k教程里. 最终我们得到两个Seurat对象

  • seu.obj1: 10x的10k pmbc
  • seu.obj2: 10x的3k pbmc

结果对比

基本上看多了PBMC数据的UMAP结果,大概就能通过UMAP的拓扑结构推断不同的细胞类型,比如说两个数据集的cluster3, 肯定都是B细胞了。

接着,我们基于文章的描述来实现代码

第一步是对数据按照cluster进行合并,并进行标准化

To identify correlated cell types between two cell atlas datasets, we first aggregated the cell-type-specific UMI counts, normalized by the total count, multiplied by 100,000 and log-transformed after adding a pseudocount

1
2
3
4
5
6
7
8
9
10
11
12
13
# 确保两个数据集的特征数相等
shared_gene <- intersect(rownames(seu.obj), row.names(seu.obj2))

## 计算各个cluster之和并标准化
seu.obj.mat <- Seurat::AggregateExpression(seu.obj, assays = "RNA", features = shared_gene,slot = "count")$RNA
seu.obj.mat <- seu.obj.mat / rep(colSums(seu.obj.mat), each = nrow(seu.obj.mat))
seu.obj.mat <- log10(seu.obj.mat * 100000 + 1)


seu.obj.mat2 <- Seurat::AggregateExpression(seu.obj2, assays = "RNA", features = shared_gene, slot = "count")$RNA
seu.obj.mat2 <- seu.obj.mat2 / rep(colSums(seu.obj.mat2), each = nrow(seu.obj.mat2))
seu.obj.mat2 <- log10(seu.obj.mat2 * 100000 + 1)

第二步是筛选基因. 文章为了提高准确性,先根据目标数据集中给定细胞类型和整体细胞类型中位数的倍数变化,选择前200个;然后给根据目标数据集中给定细胞类型和其他细胞类型最大值的倍数变化,选择前200个,对两个结果进行合并。(我使用差异富集的marker的前100个,发现效果也差不多)

To improve accuracy and specificity, we selected cell-type-specific genes for each target cell type by (1) ranking genes on the basis of the expression fold-change between the target cell type versus the median expression across all cell types, and then selecting the top 200 genes; (2) ranking genes on the basis of the expression fold-change between the target cell type versus the cell type with maximum expression among all other cell types, and then selecting the top 200 genes; and (3) merging the gene lists from steps (1) and (2).

1
2
3
4
5
6
7
8
9
10
11
12
13
# 以cluster3为例
cluster <- 3

seu.obj.gene <- seu.obj.mat[, cluster + 1]

gene_fc <- seu.obj.gene / apply(seu.obj.mat, 1, median)
gene_list1 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_fc <- seu.obj.gene / apply(seu.obj.mat[,-(cluster+1)], 1, max)
gene_list2 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_list <- unique(c(gene_list1, gene_list2))

第三步: 是在相关系数不为负的限制下,求解 $T_a = \beta_0a + \beta_{1a}M_b$. 其中 $M_b$ 是预测数据集的基因表达量, $\beta_{1a}$就是对应系数。

1
2
3
4
5
6
7
Ta <- seu.obj.mat[gene_list, cluster+1]

Mb <- seu.obj.mat2[gene_list,]

library(lsei)
solv <- nnls(Mb, Ta)
corr <- solv$x

这里算出预测数据集的cluster3系数最高,为0.93, 其他都是0.

第四步: 调换预测数据集和目标数据集,重新进行第二步和第三步

Similarly, we then switch the order of datasets A and B, and predict the gene expression of target cell type (Tb) in dataset B with the gene expression of all cell types (Ma) in dataset A.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#  predict a with b
seu.obj.gene <- seu.obj.mat2[, cluster]

gene_fc <- seu.obj.gene / apply(seu.obj.mat2, 1, median)
gene_list1 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_fc <- seu.obj.gene / apply(seu.obj.mat2[,-cluster], 1, max)
gene_list2 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_list <- unique(c(gene_list1, gene_list2))

Ta <- seu.obj.mat2[gene_list, cluster]
Mb <- seu.obj.mat[gene_list,]
solv <- nnls(Mb, Ta)
corr2 <- solv$x

第五步: 经过上面运算后,我们就得到对于数据集A中的各个细胞类型,数据集B各个细胞类型对应的回归系数,记作 $\beta_{ab}$ 以及对于数据集B中的各个细胞类型,数据集A各个细胞类型对应的回归系数, $\beta_{ba}$, 通过如下公式整合两个参数, $\beta = 2(\beta_{ab} + 0.01)(\beta_{ba} + 0.01)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
# 批量计算各个类群的回归系数
list1 <- list()
for (cluster in seq(1, ncol(seu.obj.mat))) {

# predict a with b
seu.obj.gene <- seu.obj.mat[, cluster]

gene_fc <- seu.obj.gene / apply(seu.obj.mat, 1, median)
gene_list1 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_fc <- seu.obj.gene / apply(seu.obj.mat[,-cluster], 1, max)
gene_list2 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_list <- unique(c(gene_list1, gene_list2))

Ta <- seu.obj.mat[gene_list, cluster]
Mb <- seu.obj.mat2[gene_list,]
solv <- nnls(Mb, Ta)
corr <- solv$x

list1[[cluster]] <- corr

}
list2 <- list()
for (cluster in seq(1, ncol(seu.obj.mat2))) {

# predict a with b
seu.obj.gene <- seu.obj.mat2[, cluster]

gene_fc <- seu.obj.gene / apply(seu.obj.mat2, 1, median)
gene_list1 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_fc <- seu.obj.gene / apply(seu.obj.mat2[,-cluster], 1, max)
gene_list2 <- names(sort(gene_fc, decreasing = TRUE)[1:200])

gene_list <- unique(c(gene_list1, gene_list2))

Ta <- seu.obj.mat2[gene_list, cluster]
Mb <- seu.obj.mat[gene_list,]
solv <- nnls(Mb, Ta)
corr <- solv$x

list2[[cluster]] <- corr

}

# 合并结果
mat1 <- do.call(rbind, list1)
mat2 <- do.call(rbind, list2)

# 计算系数
beta <- 2*(mat1 + 0.01) * t(mat2+0.01)

row.names(beta) <- paste0("C", 1:nrow(beta)-1)

colnames(beta) <- paste0("C", 1:ncol(beta)-1)
# 可视化
pheatmap::pheatmap(beta, cluster_rows = FALSE,cluster_cols = FALSE)

结果如下:

image.png

观察发现: 3k数据集UMAP上邻居的C4,6对应10k数据集的C6和C8也是临近关系, 说明基于NNLS的预测方法,也是相对比较准确。

当然文章也通过对多个数据集的相互比较验证了这个分析方法。

For validation, we first applied cell-type correlation analysis to independently generated and annotated analyses of the adult mouse kidney (sci-RNA-seq component of sci-CAR19 versus Microwell-seq). We subsequently compared cell subclusters from this study (with detected doublet-cell ratio ≤10%) to fetus-related cell types (those with annotations including the term ‘fetus’) from the Microwell-seq-based MCA. A similar comparison was performed against cell types annotated in BCA.

但是,文章并没有解释他可以用NNLS回归的方法做标记的迁移,鉴于他没有引述参考文章,我们可以认为这是作者第一次提出的标记迁移思路。所以,不妨让我做一些推测吧。

我们比较耳熟能详的分析是,线性回归里面的最小二乘法,也就是给定两组变量a和b, 通过最小二乘的方法求解a = bx + c中的系数x和c. 那么,假设存在三个细胞类型A, B,C. 其中A和B是同一类型,A和C不是同一类型。对于同一种细胞类型的两个数据,只不过由于实验,测序和分析等各种因素,导致基因表达值并不相同,但是存在着一定关联。我们可以通过求解 A = Bx + Cy 中的系数x, y,通过比较系数来评估A和哪个细胞类型最类似。进一步, 显然我们不能让一个细胞类型的贡献为负,最少是无作用,也就是0, 因此我们就得 限定系数不为负,也就是非负最小二乘(non-negative linear square, NNLS)回归。

苹果的通用控制初体验

今天(3.15)我更新了我MacOS系统到12.3, iPadOS系统到13.4, 用上了苹果的通用控制功能. 经过几个小时的使用之后, 我先简单聊聊自己使用体验.

打开通用控制非常简单, 只需要在升级后的MacOS上以此点开【设置】->【显示器】,就会看到多出了一个通用控制的按钮。

显示器

点击通用控制后,勾选前两个选项就算启用了通用控制(目前还是BETA版本)

通用控制

之后你就可以在将鼠标移动到屏幕边缘,就会惊喜的发现鼠标居然移动到另一台苹果设备上,实现对另一台设备的控制。接着,你就可以直接将各种文件拖动到另一台设备上,仿佛在使用同一台设备一样.

由于我只有一台MacBook Pro和一台iPad Pro,所以交互仅限于MacOS系统和iPadOS系统。我尝试了如下的操作

  • 使用MacBook的键盘进行交互
  • 使用外接键盘进行交互
  • 使用iPad的妙控键盘进行交互
  • 将PDF文件拖动到iPad的应用上,如MarginNote3, GoodReader
  • 将图片拖动到iPad的应用上, 如备忘录, keynote,

在键盘使用上,我感觉到很强的割裂感。为iPad Pro配备妙控键盘能够跨设备编辑,能够使用cmd + c/v 在不同设备间复制粘贴文字,使用shift 输入标点符号,使用地球键切换输入法,使用Caps切换搜狗输入法的大小写,唯一遗憾的是 没有Esc键(除非用地球键替换成Esc). 但是当你用Macbook自带的键盘或者外接的键盘时,这些常用的快捷键都失效了。此外,如果你将妙控键盘地球键替换成Esc键,你会发现MacOS的键盘的CTRL和CAPS LOCK都会变成中英文切换键。

在触控板和鼠标的使用上,妙控键盘自带触控板,MacBook的触控板和鼠标体验上差不多。

在文件交互上,通用控制表现的相对不错,我可以直接从Mac往iPad上拖动文件,也可以从iPad上【文件】应用里将文件拖动到Mac上。你可以认为是非常方便的隔空投递。

当然也有几点不足: 第一,iPad将文件拖动到屏幕边缘的时候,逻辑上很有可能会变成iPad的分屏,所以一定要拖动的快一点。第二,文件传输速度虽然有80m/s, 但是我在拖动iPad上图片到Mac上时能看到读条,有一些延迟。第三,部分应用还没有跟上这套逻辑,比如说,我没能把iPad的图片拖拽到Mac浏览器里面实现图片上传;

总的来说,我感觉目前的通用控制更像是简化版的隔空投递,以及能够自动切换设备的无线键鼠套装。如果后续能够更进一步,能够交互MacOS和iPad之间或者MacOS和MacOS之间的应用,那可能才会让我更加觉得有用,举个简单的例子,我一台设备装了Photoshop,另外的设备没有安装,最好的情况当然你直接把应用拖到另一个设备的时候,直接能在另一台设备上查看了,而不是再装一个PS。