基因组数据降采样
我的学习目的：开发工具的过程中需要去对深测序和浅测序样本进行转换。
介绍
通过阅读文献，看到了两种方法去对基因组数据降采样：
-
对reads数据进行降采样，通常输入是BAM或SAM文件，比如使用Picard这样的工具去做。
-
对mutation数据进行降采样，文献中没有给出方法[1]。
以及一个问题：为什么我们需要降采样？可能是：
- 这次使用深度测序，下次是否通过更浅的测序可以得到类似的结果，从而降低成本？
Picard下载
git clone https://github.com/broadinstitute/picard.git
cd picard/
使用gradle来build Picard：
./gradlew shadowJar
在安装过程中，我最初使用的是java11，而Picard更新后要求使用java17：
sudo apt-get install openjdk-17-jdk
调用：
java -jar build/libs/picard.jar
使用Picard
使用Picard对SAM或BAM文件进行降采样，随机保留reads的子集，从相同的模板得到的reads同时被保留或者丢弃，来满足reads是从一定概率的输入模板中得到的，这个条件。该工具提供不同的降采样策略：
-
ConstantMemory：采用hash-project strategy来限制内存，由于降采样是随机的，所以真实保留比例是围绕给定限定比例变化的。由于内存固定，适合大的输入，由于随机性，在数量多的地方准确，少的地方准确度下降。
-
HighAccuracy：尽可能接近指定的比例，需要和传入reads流中模板名称输入成比例的内存，对于大的输入，需要很大的内存。
-
Chained：综合上述两个工具的优点。对于数据量比较大的情况（比如上百万的reads）很准确，对于小的输入，推荐HighAccuracy。
例子1: 保留2%的reads
java -jar picard.jar DownsampleSam \
I=input.bam \
O=downsampled.bam \
STRATEGY=Chained \
P=0.02 \
ACCURACY=0.0001
例子2: 默认参数保留2%的reads
java -jar picard.jar DownsampleSam \
I=input.bam \
O=downsampled.bam \
P=0.02
例子3: 精确度很高地默认参数保留0.001%的reads（需要很大的内存）
java -jar picard.jar DownsampleSam \
I=input.bam \
O=downsampled.bam \
STRATEGY=HighAccuracy \
P=0.00001 \
ACCURACY=0.0000001
由作者进行了比较，和samtools相比，Picard节约了80%的时间。
我的任务：把一个平均深度为50X的pcawg样本降到2X。downsample.sh：
#!/bin/bash
# 记录开始时间
start_time=$(date +%s)
input=/home/data/sda/tzy/pcawg/73c4a73954614dda9be097dfb8305bd8.bam
output=/home/data/sda/tzy/pcawg/downsamp_73c4a73954614dda9be097dfb8305bd8.bam
java -jar /home/tzy/tools/picard/build/libs/picard.jar DownsampleSam \
I=$input \
O=$output \
P=0.04
# 记录结束时间
end_time=$(date +%s)
# 计算总时间
elapsed_time=$((end_time - start_time))
echo "进程执行时间为 $elapsed_time 秒"
run:
nohup bash downsample.sh > nohup.downsample.out &
总时长花了2个小时左右。
批量并行
任务目标：
- depth_level记录了降采样目标：0.5X, 1X, 2X, 5X, 10X
0.5
1
2
5
10
- ERR174341.ave_depth.txt 记录的是samtools获得的深度，文件内容是
Average depth: 11.5458
- 想通过depth_level除以ERR174341.ave_depth.txt中的深度，来分别降采样，并以新的降采样深度命名。
data="/home/data/sdb/tzy/NA12878/raw_data/ERR174341"
depth="/home/tzy/project/CNbenchmark/script/ERR174341/depth_level"
average_file="/home/data/sdb/tzy/NA12878/raw_data/ERR174341/ERR174341.ave_depth.txt"
export data
export depth
export average_file
# 提取 ERR174341.ave_depth.txt 文件中的数字部分
average=$(grep -oE '[0-9]+([.][0-9]+)?' "${average_file}")
export average
function parallel_down {
start_time=$(date +%s)
echo "This is $1"
number=$(echo "$1" | awk -v num="${average}" '{ if (num != 0) printf "%.2f\n", $1/num }')
echo "$number"
java -jar /home/tzy/tools/picard/build/libs/picard.jar DownsampleSam \
I=$data/ERR174341.hg19.bam \
O=$data/ERR174341.hg19.$1.bam \
P=${number}
# 记录结束时间
end_time=$(date +%s)
# 计算总时间
elapsed_time=$((end_time - start_time))
echo "进程执行时间为 $elapsed_time 秒"
}
export -f parallel_down
parallel --link parallel_down :::: ${depth}
最终并行生成以下文件，且不需要人为指定DownsampleSam中的P参数（P=目标深度/当前深度）
ERR174341.hg19.0.5.bam
ERR174341.hg19.1.bam
ERR174341.hg19.2.bam
ERR174341.hg19.5.bam
ERR174341.hg19.10.bam