不知道从什么时候开始，我的脑海里就植入了一个想法，“我要学会C语言”。虽然我在大学时学过C语言，还参加过C语言等级考试，但是我现在能写的代码也就是”hello world!”而已。你让我用C语言具体的完成一些事情，比如说读取一个FASTQ文件将其转成一个FASTA，我甚至都不会打开文件。

我买过很多C语言相关的书，比如下面这些

• 笨方法学C语言
• 啊哈C
• C和指针
• C专家编程
• C陷阱和缺陷
• C语言程序设计
• C语言从入门到精通
• C程序设计语言
• C primer Plus

有些书太难，我看着看着就困了，有些书自我感觉太简单，我看着看着就无聊了，到最后我一本书都没有看完，每当处理数据的时候，还是掏出我的Python和R吧。

你说，Python库和R包，它不香吗？为啥要折腾地去学C语言，何苦呢，何必呢？

但是我还是不甘心，还是会去看C语言相关的书，忍不住点开bwa的源代码（然后自闭）。每次都感觉自己啥都没有学进去，但其实这些内容都在潜意识中不断的加工积累。终于在不久前，我有一种感觉，我站到了C语言的目前。

这感觉就像多年前我刚开始接触RNA-seq，看书都是似懂非懂（就是那本「RNA-seq best practice」)。就跟段子写的一样，打开书，马冬梅，关上书，马什么梅？打开书，马冬梅，关上书，什么冬梅？考试，孙红雷。直到某一次生物统计课后，我在回去的路上，突然感觉一切都连接在一起，整个大脑都兴奋了起来。那一天，我才感觉自己站到了生物信息学的大门前。

当然光看到门是不够的，我还需要不断强化这种知识的联结，最好的方式就是通过写作的方式倒闭自己输入。因此，我将会更新一系列和C语言有关的内容，把自己对C语言的理解写下来。

最后，不是所有人都需要学C语言。我学C语言是为了让自己心安，只不过在学习过程中，我开始思考如何编写更高效的Python和R代码，也能解决和C语言相关的报错（比如说段错误和编译失败）。

以下，是我本次学习C语言的一些计划

• 掌握文件的读写
• 彻底掌握指针
• 学会对C语言代码进行调试
• 学会使用结构体
• 学会使用动态内存
• 在实际项目中使用C语言编写小工具
• 学习C语言库，包括不限于
  • sds
  • klib
  • htslib
• 学习使用C处理生信数据
• 学会使用多线程
• 学会和R/Python进行交互

在组装基因组之前一定要先对要组装的物种有一个大致的了解，判断其复杂程度, 标准如下

• 基因组大小：基因组越大，测序花的钱越多
• 简单基因组: 杂合度低于0.5%, GC含量在35%~65%, 重复序列低于50%
• 二倍体普通基因组: 杂合度在0.5%~1.2%中间，重复序列低于50%。或杂合度低于0.5%，重复序列低于65%
• 高复杂基因组: 杂合度>1.2% 或 重复率大于65%

k-mers估计法

最简单的策略就是基于k-mer对基因组做一个简单的了解, 使用jellyfish统计k-mers，然后作图

1
2
3
4
5
6
7

jellyfish count  -m 21 -s 20G -t 20 -o 21mer_out  -C  <(zcat test_1.fq.gz) <(zcat test_2.fq.gz)
# -m k-mers的K
# -s Hash大小, 根据文件大小确定
# -t 线程
# -o 输出前缀
# -C 统计正负链
jellyfish histo -o 21mer_out.histo 21mer_out

一些注意事项：

1. 绝对不要用--min-qual-char或其他参数，它们会将低质量的碱基替换成N
2. 在测序时由于不知道测得到底是DNA的哪一条链，因此k-mer及其互补链其实是等价的，所以一定要用-C参数

将数据导入R语言中，进行作图

1
2
3
4

pdf("21_mer.out.pdf")
dataframe19 <- read.table("21mer_out.histo")
plot(dataframe19[1:200,], type="l")
dev.off()

由于只有一个主峰，说明该物种的杂合度并不高，基本上也就是二倍体。如果图中出现多个峰，说明它可能是多倍体或者是基因组杂合度高。

基因组大小(G)估计算法为:

$$
G= K_{num} / K_{depth}
$$

其中 $K_{depth}$ 为K-mer的期望测序深度, $K_{num}$ 为K-mer的总数。 通常将K-mer深度分布曲线的峰值作为其期望深度。

基因组的杂合性和使得来自杂合片段的K-mer深度较纯合区段降低50%。如果目标基因组有一定的杂合性，会在k-mer深度分布曲线主峰位置(c)的1/2处(c/2)出现一个小峰。杂合度越高，该峰越明显。

推荐文献: Genomic DNA k-mer spectra: models and modalities

基于组装

基于K-mers可以较好的预测基因组大小，并定性的了解基因组的复杂情况，如果想更具体的了解基因组的复杂度，可以先将50X以上的段片段进行组装，然后进行分析。

组装的工具比较多，推荐用SOAPdenovo，因为速度快。

新建一个contig.config, 增加如下内容

1
2
3
4
5
6
7
8
9
10
11

max_rd_len=150
[LIB]
avg_ins=200
reverse_seq=0
asm_flags=3
rd_len_cutoff=100
rank=1
pair_num_cutoff=3
map_len=32
q1=read_1.fq
q2=read_2.fq

组装出参考序列

1

~/opt/biosoft/SOAPdenovo2/SOAPdenovo-63mer all -s contig.config -R -K 63 -p 30 -o assembly/graph

最后graph.scafSeq是拼接后的序列, 提取出大于300bp的序列.

1
2

# adjust format
bioawk -c fastx -v name=1 '{if(length($seq)>300) print ">"name "\n" $seq;name+=1}' assembly/graph.scafSeq >contig.fa

杂合度估计

将原来的序列回贴到contig上，并用samtools+bcftools进行snp calling.统计变异的碱基占总体的比例。

1
2
3
4

mkdir -p index
bwa index contig.fa -p index/contig
bwa mem -v 2 -t 10 index/contig read_1.fq read_2.fq | samtools sort -n > align.bam
samtools mpileup -f contig align.bam | bcftools call -mv -Oz -o variants.gz

一方面由于SOAPdenovo组装过程中会出错, 另一方面samtools在变异检测上也存在很高的假阳性, 所以总得先按照深度和质量过滤一批假阳性。

1
2

bcftools view -i ' DP > 30 && MQ > 30' -H variants.vcf.gz | wc -l
# 325219, 无过滤是445113

变异数目占基因组大小的比例就是杂合度。我的contig大概是200M，找到0.3M左右的变异，也就是0.0015，即0.15%.

重复序列估计

基于同源注释，用RepeatMasker寻找重复序列. 这里要注意分析的fasta的ID不能过长，也就是最好是>scaffold_1这种形式，不然会报错。

1
2
3
4
5
6

~/opt/biosoft/RepeatMsker/RepeatMasker -e ncbi -species arabidopsis -pa 10 -gff -dir ./ contig.fa
# -e ncbi
# -species 选择物种 用~/opt/biosoft/RepeatMasker/util/queryRepeatDatabase.pl -tree 了解
# -pa 并行计算
# -gff 输出gff注释
# -dir 输出路径

输出结果中主要关注如下三个

• output.fa.masked, 将重复序列用N代替
• output.fa.out.gff, 以gff2形式存放重复序列出现的位置
• output.fa.tbl, 该文件记录着分类信息

1
2
3
4
5
6
7

==================================================
file name: anno.fasta
sequences:         62027
total length:  273135210 bp  (273135210 bp excl N/X-runs)
GC level:         36.80 %
bases masked:   79642191 bp ( 29.16 %)
==================================================

也就是说我们的物种有30%的重复序列，作为参考，拟南芥125Mb 14%重复序列, 水稻389M，36%重复

附录：软件安装

安装RepeatMasker

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

cd ~/src
wget http://tandem.bu.edu/trf/downloadstrf409.linux64
mv trf409.linux64 ~/opt/bin/trf
chmod a+x ~/opt/bin/trf
# RMBlast
cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-src.tar.gz
wget http://www.repeatmasker.org/isb-2.6.0+-changes-vers2.patch.gz
tar xf ncbi-blast-2.6.0+-src
gunzip isb-2.6.0+-changes-vers2.patch.gz
cd ncbi-blast-2.6.0+-src
patch -p1 < ../isb-2.6.0+-changes-vers2.patch
cd c++
./configure --with-mt --prefix=~/opt/biosoft/rmblast --without-debug && make && make install
# RepeatMasker
cd ~/src
wget http://repeatmasker.org/RepeatMasker-open-4-0-7.tar.gz
tar xf RepeatMasker-open-4-0-7.tar.gz
mv RepeatMasker ~/opt/biosoft/
cd ~/opt/biosoft/RepeatMasker
## 解压repbase数据到Libraries下
## 配置RepatMasker
perl ./configure

今天在写代码的时候，发现项目中出现了一些重复的代码，所以要把他们封装成一个单个函数。在封装的过程中，我遇到了一个让我头疼的问题。

在使用dplyr的时候，你可能会注意到一个非常有趣的细节，那就是你不用""来区别变量和字符串，dplyr能够帮你好这个事情。举个例子，下面的代码都是让iris数据集按照”Sepal.Length”进行排序。

1
2

group_by(iris, Sepal.Length)
group_by(iris, "Sepal.Length")

这时候，让我们思考一个问题，如果在之前命名了一个group.by <- "Sepal.Length"，那么运行group_by(iris, group.by)的时候， 这个group.by会被替换成Sepal.Length吗？下面的代码会报错吗？大家可以思考一下，然后往下看。

1
2

group.by <- "Sepal.Length"
group_by(iris, group.by)

实际上，运行上面的代码，你会得到一个报错”Error: Column group.by is unknown”. group_by没有替换掉你的变量名。

为什么会出现这个情况？这个就涉及到dplyr编程的内容，具体可以参考https://cran.r-project.org/web/packages/dplyr/vignettes/programming.html

参考dplyr的教程，如果要让上面的代码能够运行，我们需要需要在变量名前加上!!或者调用UQ函数

1
2

group_by(iris, !!group.by)
group_by(iris, UQ(group.by))

假如你要创造一个函数，要用到group_by，那么你应该怎么写呢？我们的直觉就是下面的代码

1
2
3

my_group_by <- function(df, group.by){
  df <- group_by(df,group.by )
}

根据前面的铺垫，你应该知道，运行my_group_by(iris, Sepal.Length)会出现报错，报错信息为” Error: Column group.by is unknown”。 于是你试着之前的解决方法加上了”UQ”

1
2
3
4

my_group_by <- function(df, group.by){
  df <- group_by(df, UQ(group.by) )
  return(df)
}

思考下，my_group_by(iris, Sepal.Length)能够得到结果吗？

很遗憾，代码报错了

1
2

Error in splice(dot_call(capture_dots, frame_env = frame_env, named = named,  : 
 object 'Sepal.Length' not found 

正确的调用方法是my_group_by(iris, "Sepal.Length"). 当然由于你用习惯了dplyr，你希望是my_group_by(iris, Sepal.Length)调用代码，那么你的函数需要怎么写呢？为了解决这个问题，你可能要仔细阅读https://cran.r-project.org/web/packages/dplyr/vignettes/programming.html才行，理解什么叫做”Quasiquotation”

经过一波折腾之后，我受够了这种将dplyr代码改成函数时的不一致性（也可能是我不熟练），我决定还是只用R基础代码来实现我的功能，不就是分组排序吗，为啥一定要用group_by 和 arrange呢.

无非就是先利用因子将数据库分成多个列表(split)，然后对每个列表按照某一列进行排序(lapply)，而这里排序过程就是获取从最大到最小的索引(order)，最后按行进行合并(do.call, rbind)而已呀。如下是实现的代码

my.group.by <- function(df, group.by, sort.by,
                        decreasing = TRUE){
  
  df.split <- split(df, df[[group.by]])
  
  df.split.sort <- lapply(df.split, function(x){
    x.order <- order(x[[sort.by]],decreasing = decreasing)
    x <- x[x.order,]
    x
  })
   df <- do.call(rbind, df.split.sort)
  return(df)
}

my.group.by(iris, group.by = "Species", sort.by = "Sepal.Length")

一句话评价：重复序列注释用EDTA就完事了。

简介

EDTA, 全称是 Extensive de-novo TE Annotator, 一个综合性的流程工具，它整合了目前LTR预测工具结果，TIR预测工具结果，MITE预测工具结果，Helitrons预测工具结果, 从而构建出一高可信，非冗余的TE数据库，用做基因组的注释。流程图如下

安装

如果没有管理员权限，可以用conda进行安装。如果有管理员权限，可以尝试用docker或者singularity进行安装。

PS: 如果之前用过EDTA的话，可以更新一下版本，因为从1.7.0开始，EDTA把intact TE和homology based注释结合在一起，最终产生了很高质量的gff，合并了所有注释；1.7.1版把包含基因的intact也去掉了，进一步过滤。（来自于作者的建议）

conda

使用conda的安装方法如下

1
2
3
4
5
6
7

conda create -n EDTA
conda activate EDTA
python2 -m pip install --user numpy==1.14.3 biopython==1.74 pp
conda config --env --add channels anaconda --add channels conda-forge --add channels biocore --add channels bioconda --add channels cyclus
conda install -n EDTA -y cd-hit repeatmodeler muscle mdust repeatmasker=4.0.9_p2 blast-legacy java-jdk perl perl-text-soundex multiprocess regex tensorflow=1.14.0 keras=2.2.4 scikit-learn=0.19.0 biopython pandas glob2 python=3.6 trf
git clone https://github.com/oushujun/EDTA
./EDTA/EDTA.pl

需要注意的一点是，bioconda建议添加国内镜像站点, 否则可能会下载失败。

singularity

方法1: 使用EDTA上提供的docker镜像，以singularity进行安装

1

singularity build edta.sif docker://kapeel/edta

在使用时，有一点需要注意，需要用-B将外部的RepeatMasker的Libraries绑定的Libraries，否则可能会在检查依赖这一步失败。

这里，我参考的是LoReAN的方法

1
2
3
4
5
6

cd ~ # 切换到家目录
mkdir -p LoReAn && cd LoReAn
wget https://github.com/lfaino/LoReAn/raw/noIPRS/third_party/software/RepeatMasker.Libraries.tar.gz && tar -zxvf RepeatMasker.Libraries.tar.gz
singularity exec \
    -B /home/xzg/LoReAn/Libraries/:/opt/conda/share/RepeatMasker/Libraries/ \
    /home/xzg/edta.sif /EDTA/EDTA.pl -h

上面的/home/xzg/是我的家目录，需要根据实际情况进行选择

方法2: 使用 @wangshun1121 构建的docker镜像, 他解决了需要-B进行挂载的问题。

1

singularity build edta.sif docker://registry.cn-hangzhou.aliyuncs.com/wangshun1121/edta

实战

我们以拟南芥的第一条染色体为例，进行介绍

1
2
3
4
5
6

# singularity
singularity exec ~/LoReAn/edta.sif /EDTA/EDTA.pl \
    -genome chr1.fa -species others -step all -t 20
# conda
# 我的EDTA在我的家目录下
~/EDTA/EDTA.pl -genome chr1.fa -species others -step all -t 20

这里的参数比较简单,

• -genome: 输入的基因组序列
• -species: 物种名，Rice, Maize和others三个可选
• -step: 运行步骤, all|filter|final|anno, 根据具体情况选择
• -t: 线程数，默认是4

此外还有几个参数可以关注下

• -cds: 提供已有的CDS序列（不能包括内含子和UTR），用于过滤.这个值也比较重要，建议提供下，否则会降低busco值（来自于作者的推荐）

• -sensitive: 是否用RepeatModeler分析剩下的TE，默认是0，也就是不要。RepeatModeler运行时间比较久，量力而信。

• -anno: 是否在构建TE文库后进行全基因组预测，默认是0.

• -evalues 1: 默认是0，需要同时设置-anno 1才能使用。建议加上，它能够查看注释质量，是非常不错的功能哦（来自于作者的推荐）
  运行结束之后，会在当前目录下留下运行时的中间文件，保证你程序中断之后，能够断点续跑

• xxx.EDTA.raw

• xxx.EDTA.combine

• xxx.EDTA.final

以及你关注的xxx.EDTA.TElib.fa, 这就是最终的TE文库。

需要注意的是，在实际运行的时候，你不能单条染色体的运行，这不是程序设计的目的，我们这里用一条染色体仅仅是为了演示，测试程序能否顺利运行。

而在实际项目中，一定要用所有的染色体或者scaffold

可能问题

我在使用EDTA时，就遇到了两个问题。一个是singularity的EDTA直接使用时无法通过依赖检测，解决方法已经在安装部分提过，这里不在赘述。

另一个问题我在”Identify TIR candidates from scratch”这一步出现下面的报错

1
2
3

what():  terminate called after throwing an instance of 'Resource temporarily unavailable std::system_error'
what():  Resource temporarily unavailable
terminate called after throwing an instance of 'std::system_error'

我对这个报错进行了分析，找到了对应代码，即sh $TIR_Learner -g $genome -s $species -t $threads -l $maxint. 用实际内容替换变量后，即下面这行代码

1

sh /EDTA/bin/TIR-Learner2.4/TIR-Learner2.4.sh -g chr.fa -s others -t 20 -l 5000

更具体一点，可以将问题定位到脚本的Module 3, Step 3: Get dataset

1
2
3
4
5
6
7

genomeFile=/data/xzg_data/1800_assembly/annotation/repeatAnnotation/chr.fa #基因组文件的实际路径
genomeName=TIR-Learner
tir_path=/EDTA/bin/TIR-Learner2.4 # TIR-Learner2.4的路径
t=1
dir=`pwd`
export OMP_NUM_THREADS=1
python3 $tir_path/Module3_New/getDataset.py -g $genomeFile -name $genomeName -p $tir_path -t $t -d $dir"/Module3_New"

将线程数设置为1后，该代码顺利跑通。进一步，我定位getDataset.py的出问题的地方实际是predict函数。当然接着执行后续的代码，发现改动这一参数并不影响下面代码的运行。

1
2
3
4

echo "Module 3, Step 4: Check TIR/TSD"
python3 $path/Module3_New/CheckTIRTSD_M3.py -name $genomeName -p $path -t $t -d $dir"/Module3_New"
echo "Module 3, Step 5: Write to Gff"
python3 $path/Module3/WriteToGff_M3.py -name $genomeName -p $path -t $t -d $dir"/Module3_New"

我发现predict函数涉及到了Python的多进程调用，最终在偶然间找到问题真正所在，即Linux系统对用户的资源限制，可以通过ulimit -a查看。

最终我通过设置ulimit -u 9000，提高允许运行的总程序数，将问题解决。

参考资料

• EDTA官方文档
• RMblast的问题
• 一个关于fork资源不够的解决过程
• Ou, S., Su, W., Liao, Y., Chougule, K., Agda, J.R.A., Hellinga, A.J., Lugo, C.S.B., Elliott, T.A., Ware, D., Peterson, T., et al. (2019). Benchmarking transposable element annotation methods for creation of a streamlined, comprehensive pipeline. Genome Biology 20, 275.

segmentation fault (core dumped)是一个曾经让我非常头疼的报错，当然它还有一种中文形式，**段错误(核心已转移)**。

这种错误可能会在各种程序中出现，每次我都用软件名加报错信息去搜索，对于一些常见软件，我能够找到一些解决方案，于是我能按照它们提供的方式进行解决。而对于一些比较小众的软件，就几乎无解了。

直到我学习「c和指针」，我自己创造出了一个段错误，我终于不再害怕这个报错（当然也仅限于不在害怕而已）。

下面这行代码，我觉得和printf("hello world!\n")一样重要。”hello world!”让你开始入门编程，而理解下面的代码则让你深入编程。

1
2
3
4
5
6
7

#include <stdio.h>
int
main()
{
    int *a ;
    *a = 12;
}

我们可以将其复制到到一个error.c中，然后编译运行

1
2
3

gcc -o error error.c
./error
[2]    110633 segmentation fault (core dumped) ./error

在代码中，我们创建一个叫做a的指针变量，后面赋值语句则是把12存储在a所指向的内存地址，从逻辑上看好像没啥问题。其背后的问题，且听我慢慢道来。

a是一个指针变量，它和普通变量一样，它也有特定的内存地址，在那个内存递上也存放着值，该值是一个地址。对于一个普通变量而言，如果它们在声明的时候没有初始化，那么它们会有一个默认值，比如说int i;在我的机器上的结果是0。那么对于一个指针变量而言，如果没有在声明的时候没有初始化，它的值是什么呢？我用下面的代码测试了下

1
2
3
4
5
6
7
8
9
10

#include <stdio.h>

int
main(){
    int i;
    int *a;

    printf("i:%d\n",i);
    printf("a:%p\n",a);
}

在三次运行过程中，i的值始终为0，而a的值则是0x7ffee20e7dc0, 0x7ffc1d202d10,0x7fff50972560.这些都对应着内存上的一个地址。对*a进行赋值的过程，也就是对它对应的内存地址进行赋值，这就带来了一个问题。如果这些内存是合法能被我们操作的，那么这个操作将会顺利进行，不会出现报错，但是可能会出现一些你完全想象不到的问题。但如果内存是受到保护或者对应地址不存在（大部分情况），那么在Linux系统上就会跳出错误。

现在我们知道了，段错误(core dumped)是程序访问/使用它不能使用的内存地址所导致的，那么今后面对程序出现这个报错时，我们应该做呢？

1. 使用软件名 + 段错误对应的英文(core dumped)进行检索, 比如说 “STAR core dumped”, 相对而言，用英文检索的资料会比较
2. 检查自己的服务器内存是否充足，对于一些比对软件或者组装软件而言，它们比较占用内存，因此内存不够就会出现错误
3. 检查自己的输入是否合法，这可能是程序没有合理的错误捕捉机制，导致遇到错误格式的输入依旧运行了代码，结果中间指针没有正确初始化，最终出错
4. 更新软件：新的软件可能会修复之前代码错误
5. 换个软件：有些软件的版本可能就停留在它发表文章那一天，那么这个工具可能也不太可

当我开始写2019年总结的时候，我一时间不知道自己在这一年中做了些什么，总感觉自己忙忙碌碌，但是一无所获。

好在，我记录了，虽然不是以日记的形式，而是以一篇篇学习笔记的形式发布在网络上。当我将时间线拉回到2019年1月1日，然后慢慢向现在拉近时，我知道我这一年还是有一些进步的。

第一件值得开心的事，应该是我的个人博客xuzhougeng.top在7月25日的时候又上线了。这是我在2019年年初就想做的事情，我想要把自己简书上的内容完全由自己掌控，不再被无故404，目前已经有3万多人次访问。

其实搭建网站是我很小时候就想做的一件事情，那时候看到一篇文章写到用花生壳1元做自己的网站，就觉得这是一件很酷的事情。只不过无人指导，所以我一直在走弯路，但终于是把路走完了。

第二个事情，应该是我的编程水平的提高。我从年初就决定去学习数据结构和算法，只不过断断续续，最终也就是知道了大概，LeetCode上面的题目刷的还是不够。我尝试在R语言中使用Cpp提高代码运行速度，希望克服对C++的恐惧，从熟悉的语言切入慢慢掌握这门复杂的语言。我还挣扎在R语言绘图系统中，不断的看Y叔代码和ggplot2的源码，最终在翻译了ggproto相关内容后，对ggplot2和图形语法有所理解。

第三个应该是我的数据处理能力上的进步，因为踩了更多的坑，所以有了更多的项目经验

• 学习使用Nanopore做基因组组装，把相关的工具都学了一遍，Canu/NECAT/NextDenovo/FLYE
• 学习用HiC辅助基因组组装，相关工具为ALLHiC/3D-DNA/SALSA
• 学习BioNano相关的知识点，发现了BioNano在处理杂合基因组上的不足，并且尝试去解决它。
• 从年初开始的miRNA-seq项目，终于在年末和师兄的不断讨论中，理解了原理，并发现之前存在的问题。
• 学了Shiny制作网页工具，用在了和师兄合作的单细胞文章数据展示上
• 在简化基因组鉴定SNP和使用SNP进行下游分析上，我还在不断的思考
• 学习了使用ANOVA进行QTL定位，当然还有质量性状的定位方法。也写了几篇BSA的教程。

除了学习上和技术外，我还坚持锻炼，去健身房撸铁。这需要感谢我的是一个师弟，在他的督促下，我保证了每周起码两次的撸铁，否则我更喜欢在电脑前发呆。目前肚子上还有一层肉，在新的一年里可能要多跑跑步，少喝奶茶才行了。

今年差不多就是这些，对于明年，我估计也会做上面这些事情，不断的学习，不断的写作。

硬盘配额

Quota在是针对一个文件系统进行限制，而不是针对某个具体目录。也就是如果设置了/dev/sda1的quota，那么一个用户在这个/dev/sda1使用的数据量就会受限。如果你还有另外的/dev/sda2, 用户可以在超过/dev/sda1后继续使用/dev/sda2。

Quota针对整个文件系统的限制项目分为如下几个部分

• 容量限制或文件数目限制(block或inode)
• soft/hard: 超过soft时，用户会被警告，但是依然能操作，但是一旦超过hard，那么直接无法操作硬盘
• 宽限时间: 超过soft后，给用户一定时间进行数据转移

对于比较新的系统（例如CentOS7），建议在创建文件系统的时候就使用xfs，采用最新的quota策略。

Ext4文件系统的Quota配置

第一步。增加文件系统支持

1
2
3
4

mount | grep /data2
# dev/sdb2 on /data2 type ext4 (rw,relatime,seclabel,stripe=64,data=ordered)
mount -o remount,usrquota,grpquota /data2
# /dev/sdb2 on /data2 type ext4 (rw,relatime,seclabel,quota,usrquota,grpquota,stripe=64,data=ordered)

第二步: 扫描文件系统并新建Quota的配置文件

1
2
3
4
5

quotacheck -avug
# -a: 扫描所有在`/etc/mtab`内含有quota支持的文件系统
# -v: 显示扫描过程
# -u: 新建aquota.user，按照用户进行扫描记录
# -g: 新建aquota.group, 按照用户组进行扫描记录

这一步可能出现的提示信息和解决方法

如果出现下面提示，说明内核支持xfs, 建议使用xfs的quota管理。

1

quotacheck: Your kernel probably supports journaled quota but you are not using it. Consider switching to journaled quota to avoid running quotacheck after an unclean shutdown.

如果出现下面提示，说明该磁盘正在独写，使用lsof/fuser查找可能的进程并关闭，之后重新运行。

1
2

quotacheck: Cannot remount filesystem mounted on /data2 read-only so counted values might not be right.
Please stop all programs writing to filesystem or use -m flag to force checking.

下面的提示会在首次执行时，可以忽略

1
2
3
4
5
6
7
8

quotacheck: Scanning /dev/sdb1 [/data3] done
quotacheck: Cannot stat old user quota file /data3/aquota.user: No such file or directory. Usage will not be subtracted.
quotacheck: Cannot stat old group quota file /data3/aquota.group: No such file or directory. Usage will not be subtracted.
quotacheck: Cannot stat old user quota file /data3/aquota.user: No such file or directory. Usage will not be subtracted.
quotacheck: Cannot stat old group quota file /data3/aquota.group: No such file or directory. Usage will not be subtracted.
quotacheck: Checked 44 directories and 23 files
quotacheck: Old file not found.
quotacheck: Old file not found.

第三步： 启动/关闭quota

1
2

quotaon -vug /data3
quotaoff -ug /data3

第四步: 编辑账户/用户组的限值和宽限时间

1

edquota -u 用户名

编辑表分为7个字段，修改的是第三个和第四个用于限制使用数据量(单位为KB)，修改第六个和第七个限制使用的文件数目(不常用)，0表示无限制。

和quota相关的Linux命令如下:

• quotactl
• quotaon/quotaoff
• edquota
• quotacheck
• repquota
• warnquota
• setquota

XFS文件系统的Quota配置

如果使用xfs作为文件系统，那么就不能或者说没有必要采用上面的方法。xfs文件系统下，quota是文件系统的元信息，使用journaling提供更高层次的一致性保证，因此在管理上就会存在差别。

1. quotacheck对XFS系统没有效果。在mount时打开quota后，XFS会自动运行quotacheck做配额核算（quota accounting）
2. 没有必要在XFS文件系统下创建quota文件。aquota.user, aquota.group
3. 配额核算和限制增强之间存在区别。配额核算必须在mount的时候打开，但是限制增强可以在配额核算打开之后随意关闭和开启
4. 建议使用state监控XFS quota系统

挂载的时候打开quota

1
2
3

umount /mnt/disk1
mount -o uquota,gquota /dev/sbd1 /mnt/disk1
chmod 1777 /mnt/disk1

为了保证开机运行，还得修改/etc/fstab的第四列的defaults成defaults,uquota,gquota

ext4可以通过fstransform转换成xfs

直接输入xfs_quota会进入交互模式

• help: 查看帮助信息, 例如help quota
• quota: 查看配额信息
• print: 列出设备
• df/free: 查看-bir(block,inode,realtime)信息，显示方式为 -hN, -f指定文件
• quit/q： 退出交互模式

使用-x会进入专家模式， 之后通过-c "命令"执行命令, 举几个例子

1
2
3
4
5
6

# 显示使用情况
xfs_quota -x -c 'report -ugibh' /mnt/disk1
# 限制用户的文件数
xfs_quota -x -c 'limit -u isoft=5 ihard=10 user1' /mnt/disk1
# 限制用户的数据量
xfs_quota -x -c 'limit -u bsoft=1m bhard=2m user1' /mnt/disk1

其他可用命令

• path: 显示设备
• timer: 宽限时间，默认7天
• enable/disable: 开启或关闭配额增强
• off: 永久性关闭

我们可能都知道，计算机以二进制的方式存储数字，举两个例子：

对于整型的125，在十进制里是1*100 + 2*10+ 5*1，而对应的二进制形式是1111101(64 + 32 + 16 + 8 + 4 + 0 + 1)。

对于浮点型的0.125，在十进制里是1/10 + 2/100 + 5/1000，对应二进制下的0.001，即0/2 + 0/4 + 1/8。

在不溢出的情况下，基本上所有十进制上的整型都能无损转成二进制。但是对于浮点型，就不能保证所有十进制下的小数都能无损转成二进制的表示方式。举个例子，假如你要在二进制下表示0.1(1/10)，已知0.1小于我们之前的0.125，因此它不能是0.001, 那么它是0.0001吗？然而0.0001对应的是1/16，也就是0.0625，比0.1小。0.00011对应0.9375， 稍微靠近了0.1，0.000111对应0.109375比0.1大，0.0001101对应0.1015625依旧比0.1大，0.00011001对应0.09765625又靠近了0.1。也就是说，0.1在二进制中这是一个无限循环小数。

1

0.0001100110011001100110011001100110011001100110011...

这就类似于试图用小数形式描述十进制下的1/3，只可能是0.3333333…。

循环是无限的，但是计算机存储是有限的，不能用有限的空间来存放无限的循环。因此对于不同编程语言，它们都会为不同的数据类型分配不同的内存大小。对于浮点数运算而言，基本上所有的编程语言都用的是IEEE 754标准。

例如Python而言的PEP 754就是关于它的浮点型存放方式。在R里面用?as.double查看帮助文档时也能发现R也遵从IEEE 754。这里不在拓展介绍这个标准，只是为了说明我们只能用有限的空间来存储无限循环的小数。因此，大部分的小数都可能不是你看到的样子。大部分的编程语言，其实都会在浮点型运算结果中给你体现出这种差异，比方说Python

1
2

>>> 0.1 + 0.2
0.30000000000000004

看到这个结果，你就会好奇为啥0.1+0.2的结果不是0.3，不会好奇为啥03不等于0.3。但是在R语言中，你会好奇为啥0.3不等于0.3， 因为明明看起来一样啊

1
2
3
4

> 0.1 + 0.2
[1] 0.3
> (0.1+0.2) == 0.3
[1] FALSE

如果你想看到它的实际值，就需要借助于sprint这个函数了

1
2

sprintf("%.20f", 0.1 + 0.2 )
[1] "0.30000000000000004441"

因此，在R语言做浮点运算时，一定要谨慎，可以考虑用Hadley开发的dplyr

1
2

library(dplyr)
near(0.1+0.2, 0.3)

在数据分析的时候，我们的目标是找到样本之间真实的生物学差异。但是这种真实的生物学因素往往会受到各种因素影响，举几个场景

• 不同样本
• 同一样本的生物学重复
• 同一样本的技术重复
• 同一样本在同一个实验室由同一团队在不同时间点处理
• 同一细胞系/小鼠在不同实验室
• 不同建库策略，10X平台，Drop-seq, SMART2-seq
• 不同测序平台，BGI/Illumina
• 不同分析流程（甚至一个工具的多个版本，如salmon，CellRanger）

这些因素之间有些是生物学真实的差异，有些是抽样时的随机波动。有些是系统性因素，比如说批次效应(batch effect)。

从直觉上讲，最好是分析的过程中过滤掉这种批次效应。但是即便是同一个人对同一个样本做的相同实验，也有可能因为时间差异导致批次效应，我们需要对这种数据集进行批次效应校正吗？我们对批次效应进行校正的同时也会引入新的问题，它很有可能将生物学本身的差异视为批次效应，然后将其去除。因此解决同一种实验的批次效应最好的方法就是搞一套比较好的实验设计（例如将样本分开在不同的实验批次中，(Kang et al., 2018)）。但是受限于实验条件，基本上做不到这些，而且这类批次效应可能也没有那么大的影响。

我们更关注的一类批次效应，其实是不同实验室，不同建库手段，不同测序平台所引起的批次效应。当我们希望通过合并同一组织数据挖掘出更有意义的信息时，就不可避免的会发现，明明是同个组织的数据，表达量就是存在明显的差异(PCA, t-SNE降维可视化)。

之前有人用bulk RNA-seq的方法（limma, ComBat,RUVseq, svaseq）对单细胞数据进行校正，但是这些工具的基本假设都是”bulk RNA-seq数据中的细胞组成相似”，可能适用于一些数据集，但是可推广性不强(Haghverdi et al., 2018)。于是就有一些专门用于单细胞转录组批次校正的工具，这里对这些软件做一个罗列

• Seurat/Integration
• batchelor/FastMNN
• scran/mnnCorrect
• Python/BBKNN
• BBER
• Conos
• LIGER
• Harmony

最后对于不同的实验设计，仅从主观上推荐的批次矫正方法

• 技术重复: ComBat
• 细胞系，生物学重复: ComBat
• 同一个人的癌症和癌旁组织: 不校正/Harmony
• 不同实验室的同一组织: Harmony
• 同一个实验室做的不同人的样本: 不校正/Harmony

Harmony就是给soft k-mean（GMM）加了一个KL divergence，然后这个惩罚项最小化可以使得不同batch的cells尽可能分配到每个cluster内部，然后可以给每个cell分配一个到cluster的probability，那么每个cluster都会有一个和细胞数维度一样的细胞分配到这个cluster的probability，其实就是后续一个（batch number+1）维度的向量到embedding出来向量的一个回归方程的样本权（权重线性回归），简单而言它认为将不同batch的cell embedding的向量是由batch的variation和cell type的variation所决定，只要把cell type这部分找到（其实就是残差）找到就可以去除batch的影响。但是要定义样本在每个cluster里面的entropy，而且要使这个值尽可能大，所以我不知道这个文章有没有做rare celltype相关的一些测试
-@汪伟旭(复旦大学)

对于批次效应的检查，可以试试kBET，但是它的内存和运算时间和数据集大小密切相关，很有可能跑不下去。

参考文献

• Haghverdi, L., Lun, A.T.L., Morgan, M.D., and Marioni, J.C. (2018). Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat Biotechnol 36, 421–427.
• Büttner, M., Miao, Z., Wolf, F.A., Teichmann, S.A., and Theis, F.J. (2019). A test metric for assessing single-cell RNA-seq batch correction. Nat Methods 16, 43–49.
• Kang, H.M., Subramaniam, M., Targ, S., Nguyen, M., Maliskova, L., McCarthy, E., Wan, E., Wong, S., Byrnes, L., Lanata, C.M., et al. (2018). Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat Biotechnol 36, 89–94.

Seurat的分析流程有两步, 对数据的normalization和scaling. 两种的作用不同，前者是为了处理每个细胞的总count不同的问题，而后者则是让每个基因的表达量的均值为0，方差为1.

normlization对应的函数是NormalizeData，通过数据进行一些列变换，消除文库大小的影响。 它有三种方法, LogNormalize, CLR, RC

默认方法方式是LogNormalize， 即对于每个细胞，将每个基因的count除以总数，然后乘以一个scale.factor, 之后以自然对数进行转换。为了提高效率，Seurat编写了C++代码用于加速

1
2
3
4
5
6
7
8
9
10
11
12

// [[Rcpp::export]]
Eigen::SparseMatrix<double> LogNorm(Eigen::SparseMatrix<double> data, int scale_factor, bool display_progress = true){
  Progress p(data.outerSize(), display_progress);
  Eigen::VectorXd colSums = data.transpose() * Eigen::VectorXd::Ones(data.rows());
  for (int k=0; k < data.outerSize(); ++k){
    p.increment();
    for (Eigen::SparseMatrix<double>::InnerIterator it(data, k); it; ++it){
      it.valueRef() = log1p(double(it.value()) / colSums[k] * scale_factor);
    }
  }
  return data;
}

也就是说，如果你不追求效率，我们是可以用纯R代码实现。

1
2
3
4

mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat_norm <- LogNormalize(data = mat)
# LogNormalize等价于
log1p(t(t(mat) / colSums(mat)) * 10000)

Scale这一步对应的函数是ScaleData， 在它处理之后，使得每个基因在所有样本的均值是0，而方差是1。这让每个基因在下游分析中的具有相同的权重，使得高表达基因不那么显著。

而如果需要移除数据集中不需要的变异来源(unwanted sources of variation), ScaleData需要设置额外的参数vars.to.regress，但是作者更加推荐使用SCTransform。

尽管ScaleData对应的源代码非常的长，但是和数据处理相关的就是如下这几条

1
2
3
4
5
6
7
8

# Currently, RegressOutMatrix will do nothing if latent.data = NULL
data.scale <- scale.function(
  mat = object[features[block[1]:block[2]], split.cells[[group]], drop = FALSE],
  scale = do.scale,
  center = do.center,
  scale_max = scale.max,
  display_progress = FALSE
)

其中scale.function是根据数据类型来决定是FastRowScale还是FastSparseRowScale, 而这两个代码Seurat也是通过C++进行提速了。因此，如果不追求效率，还是可以用纯R代码实现。

1
2
3
4
5

pbmc_small <- ScaleData(pbmc_small)
GetAssayData(pbmc_small, "scale.data")[c("IGLL5"),1:2]
# 等价于
scaled_data <- t(scale(t(GetAssayData(pbmc_small, "data"))))
scaled_data[c("IGLL5"),1:2]

默认R语言的scale是按列处理，对于行为基因，列为样本的数据，可以直接套用。但是Seurat的输入数据是行为基因，列为样本，那么就需要按行scale。换句话说，如果你需要对行进行scale，那么你可以通过Seurat:::FastRowScale()的方式调用Seurat写的FastRowScale函数。此外Seurat编写了大量的Fast*函数，都可以尝试用在代码中。