【原】没有生物学重复怎么办，pseudo replicates了解一下

对于ATAC_seq, chip_seq等蛋白富集型实验而言，设置生物学重复是非常有必要的，通过IDR软件合并生物学重复的peak calling结果，可以得到更加稳定，更具代表性的peak。生物学重复的必要性不言而喻，但是对于某些特殊样本，确实没有生物学重复该怎么办呢？

此时可以参考Encode的ATAC分析流程，从数据层面来构建一个’假’的生物学重复。基本思想是随机抽样，从总体中随机抽取一定比例的reads。比如随机抽取50%的reads, 抽取两次就可以生成两个生物学重复，然后进行下游分析。具体的实现过程如下，需要借助两个linux下的工具

1. shuf

该命令用于随机乱序显示文件中的内容，比如一个文件中的内容如下

cat a.txt
1
2
3
4
5
6
7
8
9
10

通过shuf命令可以打乱顺序

shuf a.txt
4
1
6
10
7
8
2
3
9
5

需要注意的是，默认情况下，由于是随机打乱，每次运行结果都不相同

shuf a.txt
2
3
8
7
5
6
9
1
10
4

为了保证多次运行结果的一致性，对于随机抽样的软件而言，都有一个随机数发生器的概念，其取值相同时，可以保证每次抽样的结果一致。在shuf命令中，通过以下方式设置随机数发生器

shuf --random-source seed.txt a.txt

random-source参数的值为一个文件，利用该文件中的内容作为随机数发生器，只要每次使用同一个文件，就可以保证多次运行的结果完全一致。

2. split

该命令将一个大的文件分割成多个小文件，可以按照行数进行拆分, 基本用法如下

split -l 5 a.txt

-l参数指定每个小文件中包含的行数，运行完成后，会生成两个文件，其内容分别如下

cat xaa
1
2
3
4
5cat xab
6
7
8
9
10

通过管道将上述两个命令组合起来，就可以实现随机抽取，生成pseudo replicates的过程，具体的用法可以参考Encode的流程

https://github.com/ENCODE-DCC/atac-seq-pipeline/blob/master/src/encode_task_spr.py

核心代码如下

cmd1 = 'zcat {} | shuf --random-source=<(openssl enc '
cmd1 += '-aes-256-ctr -pass pass:$(zcat -f {} | wc -c) '
cmd1 += '-nosalt </dev/zero 2>/dev/null) | '
cmd1 += 'split -d -l {} - {}.'
cmd1 = cmd1.format(
ta,
ta,
nlines,
prefix)
run_shell_cmd(cmd1)

提前计算好需要抽取的行数，通过shuf和split的组合，就可以构建出pseudo replicates。

·end·