samtools 统计重测序数据深度 depth、depth |
您所在的位置:网站首页 › samtools参数 › samtools 统计重测序数据深度 depth、depth |
测试数据链接: 链接: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 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |