samtools 统计重测序数据深度 depth、depth -a、depth -aa的联系与区别
测试数据链接：
链接：https://pan.baidu.com/s/18_R3W1K7DYdUA9DxGZ0jxw
提取码：nklh
--来自百度网盘超级会员V6的分享
背景知识：
将质控后的fastq文件比对到参考基因组之后，有可能会有部分染色体根本比对不上。 比如参考基因组中一共有10条染色体，结果fastq的reads仅比对到8条染色体，那么就是说有两条染色体压根就没有比对上。
用一个实际的测试数据来看一下：
这里准备了两个文件：一个参考基因组fasta文件，和一个按照一般流程将fastq的reads数据比对到参考基因组后生成的排序后的bam文件，利用这两个文件可以查看参考基因组一共有多少条染色体， fastq数据一共比对上多少条染色体。
统计参考基因组的染色体数目：
(base) [b20223040323@admin1 test]$ ls GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ grep "^>" GCF_000005845.2_ASM584v2_genomic.fna ## 统计参考基因组的染色体数目 >chr1 >chr2 >chr3
可见该参考基因组一共有三条染色体。
利用bam文件统计一共比对上多少条染色体：
(base) [b20223040323@admin1 test]$ ls GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ samtools view SRR1770413.sorted.bam | cut -f 3 | uniq ## 利用bam文件统计一共比对上多少条染色体 chr1 chr2 *
可见一共比对上的染色体数有两条，也就是说chr3并没有比对上。
比较samtools计算测序深度时：depth、depth -a、depth -aa的区别。
001、首先先统计一下参考基因组中每一条染色体的长度，可以使用samtools软件实现：
(base) [b20223040323@admin1 test]$ ls GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ samtools faidx GCF_000005845.2_ASM584v2_genomic.fna ### 统计每一条染色体的长度 (base) [b20223040323@admin1 test]$ ls GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam GCF_000005845.2_ASM584v2_genomic.fna.fai (base) [b20223040323@admin1 test]$ cat GCF_000005845.2_ASM584v2_genomic.fna.fai ## 查看，第二列是每一条染色体的长度 chr1 4641652 6 80 81 chr2 39862 4699685 80 81 chr3 398 4740052 80 81
可见chr1的长度为：4641652。。。。。
002、利用samtools的depth、depth -a、depth -aa分别对同一个bam文件进行测序深度统计
(base) [b20223040323@admin1 test]$ ls GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam GCF_000005845.2_ASM584v2_genomic.fna.fai (base) [b20223040323@admin1 test]$ samtools depth SRR1770413.sorted.bam > depth01.txt ## depth (base) [b20223040323@admin1 test]$ samtools depth -a SRR1770413.sorted.bam > depth02.txt ## depth -a (base) [b20223040323@admin1 test]$ samtools depth -aa SRR1770413.sorted.bam > depth03.txt ## depth -aa (base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ wc -l depth0* ## 统计生成的depth文件的行数，可见三个文件行数不一致 4652444 depth01.txt 4681514 depth02.txt 4681912 depth03.txt 14015870 总用量
003、samtools的depth参数统计的是比对到参考基因组的染色体上所有位点测序深度大于0的所有测序深度数据，用一下脚本简单验证：
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ head -n 5 depth01.txt ## depth01.txt文件记录的是位点的测序深度，每行一个位点 chr1 1 12 chr1 2 13 chr1 3 14 chr1 4 15 chr1 5 15 (base) [b20223040323@admin1 test]$ awk '{print $3}' depth01.txt | sort -n | head -n 3 ## 对第三列深度信息进行排序，可见最小深度为1 1 1 1 (base) [b20223040323@admin1 test]$ awk '{print $3}' depth01.txt | sort -n | tail -n 3 ## 输出最大测序深度 509 510 510
可见 samtools depth参数输出测序深度的结果最小的测序深度为1.
004、samtools depth -a参数输出的结果为fastq比对到参考基因组的染色体上的所有位点的测序深度，即在depth参数输出结果的基础上，也输出测序深度为0的位点。可用如下脚本进行验证：
a、如果输出fastq数据比对到参考基因组的染色体上所有位点的测序深度， 那么depth02.txt的行数应该等于比对到染色体长度的总和，即chr1和chr2的总长度。
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ wc -l depth02.txt 4681514 depth02.txt (base) [b20223040323@admin1 test]$ cat GCF_000005845.2_ASM584v2_genomic.fna.fai chr1 4641652 6 80 81 chr2 39862 4699685 80 81 chr3 398 4740052 80 81 (base) [b20223040323@admin1 test]$ awk 'NR <= 2 {sum += $2} END { print sum }' GCF_000005845.2_ASM584v2_genomic.fna.fai ##统计比对到的染色体的总长度 4681514
可以看到depth02.txt的行数等于比对到的染色体的总长度。
b、如果depth -a参数是在depth的基础上增加输出了测序深度为0的位点，那么depth02.txt 和 depth01.txt相比，相同的部分和depth01.txt完全一致， 不同的位点的测序深度都为0.
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ sort depth01.txt | md5sum ## depth01.txt排序，并计算MD5 53c631aeb9024f330280bb3762310189 - (base) [b20223040323@admin1 test]$ cat depth01.txt depth02.txt | sort | uniq -d | md5sum ## 去depth01.txt和depth02.txt的交集并计算MD5 53c631aeb9024f330280bb3762310189 - (base) [b20223040323@admin1 test]$ cat depth01.txt depth01.txt depth02.txt | sort | uniq -u | awk '{print $3}' | uniq -c ## 查看depth02.txt中多处的行，并查看测序深度 29070 0
可以看到depth01.txt和depth02.txt中的重合部分和depth01.txt完全一致，同时depth02.txt中多出的部分的测序深度都为0.
005、samtools depth -aa参数输出fasta文件中所有染色体上的位点的测序深度，包括没有比对上的染色体的位点上的测序深度。
a、如果depth -aa输出fasta文件中所有染色体上的位点的测序深度，那么depth03.txt的行数等于fasta文件中所有染色体的总长度：
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ wc -l depth03.txt ## 统计行数 4681912 depth03.txt (base) [b20223040323@admin1 test]$ cat GCF_000005845.2_ASM584v2_genomic.fna.fai chr1 4641652 6 80 81 chr2 39862 4699685 80 81 chr3 398 4740052 80 81 (base) [b20223040323@admin1 test]$ awk '{sum += $2} END {print sum}' GCF_000005845.2_ASM584v2_genomic.fna.fai ## 统计所有染色体的总长度 4681912
可以看到depth03.txt的长度等于fasta文件中所有染色体的总长度。
b、depth03.txt和depth02.txt相比，如果多出的是没有比对上的染色体的总长度，那么depth03.txt比depth02.txt多出没有比对到的染色体的长度，且多出位点的测序深度均为0.
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ echo "$(wc -l < depth03.txt) - $(wc -l < depth02.txt)" | bc ## 计算行数只差 398 (base) [b20223040323@admin1 test]$ cat GCF_000005845.2_ASM584v2_genomic.fna.fai chr1 4641652 6 80 81 chr2 39862 4699685 80 81 chr3 398 4740052 80 81 (base) [b20223040323@admin1 test]$ cat depth02.txt depth02.txt depth03.txt | sort | uniq -u | awk '{print $3}' | uniq -c ## 输出depth03.txt相对于depth02.txt多出的位点的测序深度信息 398 0
c、depth03.txt 和 depth01.txt相对，多出额外的一条染色体的测序深度信息，且测序深度均为0.
(base) [b20223040323@admin1 test]$ ls depth01.txt depth03.txt GCF_000005845.2_ASM584v2_genomic.fna.fai depth02.txt GCF_000005845.2_ASM584v2_genomic.fna SRR1770413.sorted.bam (base) [b20223040323@admin1 test]$ cut -f 1 depth01.txt | uniq ## 输出染色体 chr1 chr2 (base) [b20223040323@admin1 test]$ cut -f 1 depth03.txt | uniq ## 输出染色体 chr1 chr2 chr3 (base) [b20223040323@admin1 test]$ cat depth01.txt depth01.txt depth03.txt | sort | uniq -u | awk '{print $3}' | uniq -c ## 输出depth03.txt相对于depth01.txt额外的位点的深度信息 29468 0
可以看到depth03.txt先对于depth01.txt多出了一条染色体，且多出的位点的测序深度均为0.
小结：
samptools depth： 输出fastq数据比对到染色体的位点深度大于0的深度信息。
samtools depth -a：输出fastq数据比对到染色体的所有位点的测序深度信息，包含0深度位点。
samtools depth -aa：输出fasta文件中所有染色体上位点的深度信息，包含fastq数据没有比对到的染色体。
