VCF变异统计及绘制染色体分布图

Posted 2023-04-21

参考技术A VCF文件的变异分布图，除了使用Circos图之外，统计特定材料间的变异分布情况时可能会用到染色体分布图，为了节省时间，实现统计及绘图的自动化，现提供统计数量用的Python脚本和绘图使用的R脚本，仅供参考。如有其它更快捷的方法，或觉得脚本太拙劣，欢迎留言批评交流。

PYTHON脚本，使用脚本整理100Kb滑窗中变异位点数量：

下面画图，使用R脚本如下：

GATK3.2.2小结(转载）

http://blog.csdn.net/skenoy/article/details/38346489

经过几天的摸索和网上资料的查询对GATK软件有点小心得，现总结如下：

 

1. fasta文件最好用定位到染色体上的数据，可以不用注释VCF文件（GVF），但如果用VCF文件保证以下几个条件：

1）VCF染色体必须和fasta的染色体数目一致，顺序一致

2）VCF的位点必须从小到大排序

3）VCF的碱基有可能有其他符号，如“~”等，要去除干净

 

2. 做之前分别使用bwa index，picard中的CreateSequenceDictionary.jar和samtools中的faidx对fasta文件建立索引，且最好在fasta同一个文件夹下面

 

3. bwa做比对时，最好加入-r参数："@RG\tID:name\tLB:name\tPL:ILLUMINA\tSM:name"，为了以后不再加入头文件

 

4. picard中ReorderSam.jar是为了矫正你的sam文件的头文件与fasta相一致，如果一致，可以不用做这一步

 

5. 使用picard处理bwa的paired的sam或bam的任意程序，最好加入VALIDATION_STRINGENCY=LENIENT，因为paired reads有一条比对到染色体的末端时，另外一条picard无法识别就会报错终止运行

 

6. 如果说合并样本call variant，GATK的多线程有两个，nt代表几个样本使用一个CPU；ncr代表一个样本使用几个CPU

 

7. GATK 3.0以后不再支持ReduceReads这个程序

 

最新补充：

8. 有时候reads的cigar值会出问题、或者质量值和碱基对不上、又或者reads出现其他符号，加入下列参数：-filterRNC -filterMBQ -filterNoBases -rf UnmappedRead -rf BadMate -rf DuplicateRead -rf NotPrimaryAlignment -rf MappingQualityUnavailable

 

现阶段没有做质量值矫正和变异矫正，一是要求数据量比较大，如果小于100M的reads就不要做了；二是目前的商业项目很难做如此麻烦的处理，除了人的项目，因为有相应的很多的注释文件

当然还有其他方法进行矫正，比如跟samtools mpileup的结果相一致的才认为是可靠的