摘要

接到一个wgs项目，要帮助客户统计vcf文件中snp突变频率，比较两个样品的突变位点。这个工作在上一个项目中是手动处理的，当时参考序列短，突变位点少。这次经过比对后，发现了有个样品有上万个snp位点，肯定不能用手动处理的方式。因此，写了一个脚本来统计各个样品的突变频率。
需要统计的信息
包括染色体，突变位置，参考位点，各样品突变位点，突变率（AD杂合位点覆盖度/DP总覆盖度）

环境与方法

python 3.7
R version 3.6.1

使用代码

统计各样品各位点突变频率
我们输入的vcf文件为GATK经过去重过滤后生成的，

bash python snp_frequency.py XX.snp.vcf

  import csvimport sysvcf_file = sys.argv[1] #文件调用vcf文件genome_file = open(vcf_file,'r')gtf_file = vcf_file.replace('_snp.vcf','')# 去除文件名后缀genome_line = genome_file.readlines()#genome_list = genome_line.split()newgtf_name = '%s_stat'%gtf_file #保留样品名desktop_path = './'file_path = desktop_path+newgtf_name+'.csv' #确定文件路径和文件名，生成文件file = open(file_path,'w')for i in genome_line:snp_pos = i.split('\t')if "#" in snp_pos[0] and len(snp_pos) > 8:title = snp_pos[0]+'\t'+snp_pos[1]+'\t'+snp_pos[3]+'\t'+snp_pos[4]+'\t'+'%s_FRE'%vcf_file+'\n' file.write(title)if "#" not in snp_pos[0] and len(snp_pos) > 8:info = snp_pos[9].split(':')AD = info[1].split(',')AD_0 = int(AD[0])#转化为数字，后面需要加和AD_1 = int(AD[1])AD_total = int(AD[0])+int(AD[1]) #杂合子覆盖度求和DP = info[2]out_line = snp_pos[0]+'\t'+snp_pos[1]+'\t'+snp_pos[3]+'\t'+snp_pos[4]+'\t'+'%.0f'%AD_total+'/'+DP+'\n' #提取信息列#print(out_line)file.write(out_line)file.close()genome_file.close()

将比较样品进行合并，比较不同样品的突变情况。

WT=read.csv("USA_300_WT_stat.csv",header=T,sep='\t',check.names=F) S1=read.csv("USA_300_S1_stat.csv",header=T,sep='\t',check.names=F)
S2=read.csv("USA_300_S2_stat.csv",header=T,sep='\t',check.names=F)
WT_1 <- merge(WT,S1, all = TRUE) #合并WT和S1样品突变位点信息，不会对已有列进行重复（chr、pos、REF、ALT等）
WT_2 <- merge(WT,S2, all = TRUE) #同上
write.csv(WT_2,"USA_WT_vs_S2_snpstat.csv")
write.csv(WT_1,"USA_WT_vs_S1_snpstat.csv")

分析结果

得到.csv格式文件，这里用editplus打开

总结

WGS的最初的目的是为了寻找突变位点，这里是进行样品间比较，客户更需要的数据不是都发生突变的位点，而是某一个位点在其中一个样品发生突变（比如S1有S2无，反之亦然）。这样可以找到不同处理下的突变位点，深入研究。

2021.06.08|提取、比较各样品vcf文件中snp突变频率相关推荐

Python 导出手机通讯录文件 VCF 文件中的手机号码
文章目录 Python 导出手机通讯录文件 VCF 文件中的手机号码 1.代码 Python 导出手机通讯录文件 VCF 文件中的手机号码 1.代码 if __name__ == '__main__' ...
python提取文本中的手机号_Python从vcf文件中读取手机号并进行去重操作
文章目录 1. Python代码 file = open('test.vcf', 'r', encoding='utf-8') tels = [] for line in file: line = l ...
VCF文件中QUAL和GQ的区别
最近开始分析vcf文件, 于是去搜了相关VCF格式解读的博客. 大部分关于这两个指标的解读如下,都是描述质量值的,但也没说具体啥区别. QUAL:Phred格式(Phred_scaled)的质量值,表 ...
【最新实用版】Python批量将pdf文本提取并存储到txt文件中
#注意:笔者在2021/11/11当天调试过这个代码是可用的,由于pdfminer版本的更新,网络上大多数的语法没有更新,我也是找了好久的文章才修正了我的代码,仅供学习参考. 1.把pdf文件移动到本 ...
shell高级技巧：提取vcf文件中一个contig
这是一个很小众的需求.大部分变异检测都是基于组装质量比较高的基因组,而不是那种初步拼接的contig. 由于初步拼接的参考序列通常会有成千上万个contig序列,也就导致在VCF的头文件的##cont ...
使用bcftools提取指定样本的vcf文件（extract specified samples in vcf format）
1.下载安装bcftools. 2.准备样本ID文件,这里命名为samplelistname.txt,一个样本一行,如下所示: sample1 sample2 sample3 3.输入命令: bcft ...
基因数据处理77之从vcf文件中提取某条染色体的数据
1.代码: /*** @author xubo*/ package org.gcdss.cli.vcfimport org.apache.spark.{SparkConf, SparkContext} ...
提取Insight-MVT_Annotation_Train 数据集标签xml文件中的信息
Insight-MVT_Annotation_Train 数据集标签xml文件中的信息从xml文件中解析出所要的信息 type height width top left 写成Pasc ...
利用PyVCF模块处理VCF文件
利用PyVCF模块处理VCF文件转载自:微信公众号生信说欢迎大家扫码关注工欲善其事,必先利其器. VCF,全称Variant Call Format,是生物信息学领域最常用的遗传突变存储格式. ...

2021.06.08|提取、比较各样品vcf文件中snp突变频率

目录

摘要

环境与方法

使用代码

分析结果

总结

2021.06.08|提取、比较各样品vcf文件中snp突变频率相关推荐

最新文章

热门文章