有粉丝邮件求助,给了我两个vcf文件,旧的vcf文件走的是标准的bwa gatk流程,参考基因组是hg19,新的文件参考基因组是hg38,也是gatk标准流程。想有比较它们,首先得保证两个vcf文件的参考基因组一致,因为版本不一致,所以需要使用CrossMap等软件进行参考基因组版本转换,然后里使用 SnpSift 软件的 Concordance 命令比较它们。
如果简单看文件数量,当然是不行的。
代码语言:javascript复制137352 new.filter.sort.vcf
46976 old.filter.sort.vcf
首先看看两个vcf文件的染色体分布情况
代码语言:javascript复制cat old.filter.sort.vcf |grep -v '^#' |cut -f 1 |sort |uniq -c > old.chr
cat new.filter.sort.vcf |grep -v '^#' |cut -f 1 |sort |uniq -c > new.chr
paste new.chr old.chr
结果如下:
代码语言:javascript复制14341 chr1 4912 chr1
6562 chr10 2108 chr10
7427 chr11 2920 chr11
7436 chr12 2519 chr12
3200 chr13 856 chr13
4051 chr14 1408 chr14
4568 chr15 1520 chr15
5934 chr16 2223 chr16
7100 chr17 2702 chr17
2596 chr18 778 chr18
8020 chr19 3399 chr19
10160 chr2 3146 chr2
3349 chr20 1135 chr20
1908 chr21 666 chr21
3262 chr22 1095 chr22
7676 chr3 2708 chr3
6437 chr4 1990 chr4
6134 chr5 1975 chr5
6265 chr6 2718 chr6
7095 chr7 2115 chr7
4813 chr8 1595 chr8
6162 chr9 1915 chr9
31 chrM
2081 chrX 519 chrX
157 chrY 7 chrY
可以看到,新的vcf文件的突变位点数量远大于旧的vcf文件。仔细查看新vcf文件,**发现是没有做基本过滤,比如测序深度大于20等等指标。**所以我就顺便把它过滤一波,代码如下:
代码语言:javascript复制input=new.filter.sort.vcf
cat $input | java -jar ~/biosoft/snpEff/SnpSift.jar filter "( DP > 20 & FILTER = 'PASS' )" |
perl -alne '{print unless $F[0] =~ /_/}' |
awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1 -k2,2n"}' |
grep -v '1/2' > filtered.chr
# 然后再统计这个vcf文件的染色体发布情况
cat filtered.vcf |grep -v '^#' |cut -f 1 |sort |uniq -c > filtered.chr
paste new.chr old.chr filtered.chr
如下所示:
代码语言:javascript复制14341 chr1 4912 chr1 5731 chr1
6562 chr10 2108 chr10 2617 chr10
7427 chr11 2920 chr11 3419 chr11
7436 chr12 2519 chr12 2971 chr12
3200 chr13 856 chr13 1070 chr13
4051 chr14 1408 chr14 1636 chr14
4568 chr15 1520 chr15 1780 chr15
5934 chr16 2223 chr16 2555 chr16
7100 chr17 2702 chr17 3216 chr17
2596 chr18 778 chr18 957 chr18
8020 chr19 3399 chr19 3960 chr19
10160 chr2 3146 chr2 3822 chr2
3349 chr20 1135 chr20 1448 chr20
1908 chr21 666 chr21 722 chr21
3262 chr22 1095 chr22 1354 chr22
7676 chr3 2708 chr3 3060 chr3
6437 chr4 1990 chr4 2405 chr4
6134 chr5 1975 chr5 2334 chr5
6265 chr6 2718 chr6 2376 chr6
7095 chr7 2115 chr7 2615 chr7
4813 chr8 1595 chr8 1897 chr8
6162 chr9 1915 chr9 2497 chr9
2081 chrX 519 chrX 552 chrX
157 chrY 7 chrY 6 chrY
可以看到,经过质控后的两个vcf文件,至少是从染色体上**变异位点记录的数量上来说,非常类似的。**后面我们还需要进行更细致的探索。
再看基础的vcf统计命令
使用snpEff可以对vcf进行基础统计,它也会对染色体分布情况进行统计,更重要的有丰富的可视化图表
代码语言:javascript复制java -jar ~/biosoft/snpEff/snpEff.jar GRCh38.86
old.filter.sort.vcf > old.filter.sort.eff.vcf
java -jar ~/biosoft/snpEff/snpEff.jar GRCh38.86
new.filter.sort.vcf > new.filter.sort.eff.vcf
结果如下:
可以看到突变位点区域分类情况:
突变位点区域分类百分比
可以看到,两个vcf文件的变异位点在intron和exon区域的比例差异是最大的,其实是因为它们两个区域本来就长度很大。
另外一个统计指标
image-20200711195600818
最后看专业的软件进行两个vcf文件比较
这里使用 SnpSift 软件的 Concordance 命令,代码如下:
代码语言:javascript复制 java -Xmx1g -jar ~/biosoft/snpEff/SnpSift.jar concordance -v
new.filter.sort.vcf old.filter.sort.vcf 1>comp.results.txt 2>comp.log.txt
得到的结果文件不是很好理解:
代码语言:javascript复制2.2K Jul 11 11:57 comp.log.txt
4.2M Jul 11 11:57 comp.results.txt
590 Jul 11 11:57 concordance_new_old.by_sample.txt
89 Jul 11 11:57 concordance_new_old.summary.txt
其中 concordance_new_old.by_sample.txt 文件需要倒置查看:
代码语言:javascript复制 5 MISSING_ENTRY_new/MISSING_ENTRY_old 0
6 MISSING_ENTRY_new/MISSING_GT_old 0
7 MISSING_ENTRY_new/REF 0
8 MISSING_ENTRY_new/ALT_1 4842
9 MISSING_ENTRY_new/ALT_2 1971
10 MISSING_GT_new/MISSING_ENTRY_old 0
11 MISSING_GT_new/MISSING_GT_old 0
12 MISSING_GT_new/REF 0
13 MISSING_GT_new/ALT_1 0
14 MISSING_GT_new/ALT_2 0
15 REF/MISSING_ENTRY_old 0
16 REF/MISSING_GT_old 0
17 REF/REF 0
18 REF/ALT_1 0
19 REF/ALT_2 0
20 ALT_1/MISSING_ENTRY_old 7855
21 ALT_1/MISSING_GT_old 0
22 ALT_1/REF 0
23 ALT_1/ALT_1 22538
24 ALT_1/ALT_2 4
25 ALT_2/MISSING_ENTRY_old 7001
26 ALT_2/MISSING_GT_old 0
27 ALT_2/REF 0
28 ALT_2/ALT_1 20
29 ALT_2/ALT_2 17552
30 ERROR 0
可以看到 ALT_1/ALT_1有22538个,ALT_2/ALT_2有17552个,说明新旧vcf文件的一致性还挺好的!
然后ALT_1/MISSING_ENTRY_old是7855个,ALT_2/MISSING_ENTRY_old是7001个,说明新的vcf文件,多找出来了近1.5万个突变位点。
而MISSING_ENTRY_new/ALT_1的4842个,MISSING_ENTRY_new/ALT_2是1971个,说明旧的vcf文件,多找出来了不到7千个突变位点。
有意思的是ALT_1/ALT_1 22538
两个流程不可能完全一致,近4万个位点在两个vcf文件里面都有,超过80%的一致性了。挺好的。
如果要具体看哪些位置是一致的,哪些是冲突的,就需要去重点理解 comp.results.txt 文件。
其实可以分染色体查看
前面的 SnpSift 软件的 Concordance 命令,最后得到的是vcf文件全部染色体的,一致性的韦恩图罢了。
但是可以继续细致的探索 comp.results.txt 文件,拆分染色体后,继续统计上面提到的6种情况发生的频次。那就出一个学徒作业吧,比较两个vcf文件,然后区分染色体绘制韦恩图。
这两个vcf文件可以是不同人的,也可以是同一个人的不同批次测序或者不同数据分析流程拿到的vcf文件。
也有很多其它轮子
比如 vcf-compare 工具,bedtools等等
实际上考验的就是Linux知识
再怎么强调生物信息学数据分析学习过程的计算机基础知识的打磨都不为过,我把它粗略的分成基于R语言的统计可视化,以及基于Linux的NGS数据处理:
《生信分析人员如何系统入门R(2019更新版)》
《生信分析人员如何系统入门Linux(2019更新版)》
Linux的6个阶段需要跨越过去 ,一般来说,每个阶段都需要至少一天以上的学习:
- 第1阶段:把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅,主要目的就是去可视化,熟悉黑白命令行界面,可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。
- 第2阶段:做到文本文件的表格化处理,类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余,查找,切割,替换,合并,补齐,熟练掌握awk,sed,grep这文本处理的三驾马车。
- 第3阶段:元字符,通配符及shell中的各种扩展,从此linux操作不再神秘!
- 第4阶段:高级目录管理:软硬链接,绝对路径和相对路径,环境变量。
- 第5阶段:任务提交及批处理,脚本编写解放你的双手。
- 第6阶段:软件安装及conda管理,让linux系统实用性放飞自我。
很多粉丝发邮件询问我具体的软件,命令报错该如何解决,这些问题通常是输入文件,命令参数问题,还有java环境问题,太细致了,我没有空去帮大家debug哦。
如果是是Linux基础知识掌握的不牢固,其实拼命学习即可:
- 生信入门环境的:https://www.bilibili.com/video/BV1cJ411e7UH
- 生信服务器的:https://www.bilibili.com/video/BV1XW411d7Bp
- 生信软件的:https://www.bilibili.com/video/BV1ps411M7UZ