vcf文件提取信息及排序

vcf文件示例

VCF文件简介:
VCF(Variant Call Format)文件是一种常用的存储基因组变异信息的文本文件格式。它由一系列的元数据行(以"#"开头)和数据行组成。下面是VCF文件的主要组成部分:

  1. 元数据行(Metadata lines):以"#"开头的行,用于描述VCF文件的相关信息,如文件版本、参考基因组等。其中,以"##"开头的行为注释行,用于描述或定义特定的元数据。例如,##fileformat用于指定文件的格式,##INFO用于定义变异的相关信息。
  2. 列标题行(Header line):以"#"开头的非注释行,用于描述数据列的名称和格式。其中,前九列(CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO、FORMAT)为必需的列,其他列用于存储样本相关的信息。
  3. 数据行(Data lines):由列标题行之后的每一行组成,用于存储具体的变异信息。每一列对应于一个样本,其中前九列是固定的,后续列用于存储每个样本的基因型信息。
  4. INFO字段:位于每一行的INFO列中,用于存储与变异相关的附加信息。它以分号分隔的键值对形式表示,如"DP=20;AF=0.5",其中DP表示测序深度,AF表示突变等位基因的频率。
  5. FORMAT字段:位于列标题行的FORMAT列中,用于定义每个样本的基因型信息的格式。它以冒号分隔的关键字表示,如"GT:DP:GQ",其中GT表示基因型,DP表示测序深度,GQ表示基因质量。
    VCF文件是广泛用于存储和共享基因组变异数据的标准格式,具有可读性和灵活性,方便进行基因组变异的注释和分析。

因为我的有不同的变异类型,所有用正则把它们区分出来

open IN,"<$ARGV[0]";
open OU,">$ARGV[1]";
while (){
if($_=~/^#/){next};
chomp;
my @list=split;
my $chr=@list[0],my $pos=@list[1];

my @r=@list[0]=~/(/d+)/g;
my $r=@r[0];


  # 使用正则表达式匹配冒号前的字符
my @match=@list[7] =~ /;(.*?)=/g;
my $match=@match[0];
# 替换字符VDB为SNP
$match =~ s/VDB/SNP/g;  
# 替换字符IS为INDEL
$match =~ s/IS/INDEL/g;

print OU "$chr/t$r/t$pos/t$match/n";
}
close IN;
close OU;
匹配的结果

这里我想看变异来自哪套染色体,所以要排序,来自A和B染色体的区分开,所以这里只排序第二和第三列的内容:

sort -k2,2n -k3,3n t1 >t2
排序后的结果
import os

# 设置目录路径
directory = './'

# 统计第1列含有字母A到I的数量的字典
letter_count = {'A': 0, 'B': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0}

# 统计第一列相同时第四列中INDEL和SNP的数量的字典
indel_snp_count = {}

# 遍历目录中的文件
for filename in os.listdir(directory):
    filepath = os.path.join(directory, filename)
    
    # 只处理bed文件
    if filename.endswith('.bed'):
        # 统计第1列含有字母A到I的数量
        with open(filepath, 'r') as file:
            for line in file:
                fields = line.strip().split('/t')
                first_column = fields[0]
                fourth_column = fields[3]
                
                # 判断第1列是否含有字母A到I
                for letter in letter_count.keys():
                    if letter in first_column:
                        letter_count[letter] += 1
                        break
                    
                # 统计第一列相同时第四列中INDEL和SNP的数量
                if first_column in indel_snp_count:
                    if fourth_column == 'INDEL':
                        indel_snp_count[first_column]['INDEL'] += 1
                    elif fourth_column == 'SNP':
                        indel_snp_count[first_column]['SNP'] += 1
                else:
                    indel_snp_count[first_column] = {'INDEL': 0, 'SNP': 0}
                    if fourth_column == 'INDEL':
                        indel_snp_count[first_column]['INDEL'] = 1
                    elif fourth_column == 'SNP':
                        indel_snp_count[first_column]['SNP'] = 1

# 输出第1列含有字母A到I的数量统计结果
print("第1列含有字母A到I的数量统计:")
for letter, count in letter_count.items():
    print(f"{letter}: {count}")

# 输出第一列相同时第四列中INDEL和SNP的数量统计结果
print("第一列相同时第四列中INDEL和SNP的数量统计:")
for first_column, counts in indel_snp_count.items():
    print(f"第一列为{first_column}时:")
    print(f"INDEL数量:{counts['INDEL']}")
    print(f"SNP数量:{counts['SNP']}")
    print()

版权声明:
作者:admin
链接:https://www.techfm.club/p/99219.html
来源:TechFM
文章版权归作者所有,未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>