Linux环境下利 用perl脚本批量筛选VCF文件指定样本
发布时间：2023-02-09 11:10:25 所属栏目：Linux 来源：
导读： 今天分享一个实用的数据提取思路：从VCF文件中提取指定列，比如有1000个样本，需要从中提取出200个制定样本，形成新的vcf文件。需要在Linux环境下进行操作，主要用到perl脚本和shell脚本。
原始数据
原始数据
|
今天分享一个实用的数据提取思路：从VCF文件中提取指定列，比如有1000个样本，需要从中提取出200个制定样本，形成新的vcf文件。需要在Linux环境下进行操作，主要用到perl脚本和shell脚本。 原始数据 原始文件vcf ：开始为若干行注释信息以##开头，然后是表头信息由#开始（该行包括样本ID） 待提取的样本ID：txt文件，每个样本ID一行（list.txt） 什么是VCF文件？ 伴随着大规模的基因分型及测序数据的产生，迫切需要一种新的格式来记录高效的记录这些变异信息。VCF(Variant Call Format)就是这样一种用来贮存基因序列变异信息的文本文件。 VCF 格式文件信息：header section：元信息(meta-information)，以‘##’为前缀，通常包含fileformat、fileDate、reference等信息。头行信息( header line )，以‘#’为前缀data section：该部分为主题部分，记录了每个样品每个位点处的基因分型信息。可以看做是一个大的Excel表格，每行是一个位点，每列是一个样本。 linux 压缩文件夹命令_linux压缩命令tar.gz_linux压缩 linux 压缩文件夹命令_linux压缩命令tar.gz_linux压缩 创建一个工作文件夹work，将第一个文件放在其data子文件夹下，第二个文件放在工作文件夹work中。 一般每个染色体有一个VCF文件，压缩为xxx.vcf.gz，将全部压缩包保存到一个目录data下。 操作步骤批量解压缩 在data文件夹中，输入gunzip *.vcf.gz解压所有文件，获得很多个VCF文件 保存注释信息 由于VCF文件刚开始几十行为注释信息（假设为78行），可以先把注释信息截取出来单独存放为一个zhushi.vcf，然后再对剩下的数据矩阵进行转置筛选，最后将注释信息加上。 head -n 78 ./data/xxx.vcf > zhushi.vcf 转置VCF文件 由于VCF文件的每一列为一个样本，为了进行筛选，先将其转置linux压缩，变为每行一个样本的格式，使用如下perl脚本（s1.pl）： #! /usr/bin/perl -w use strict; die "perl $0 \n" unless @ARGV==1; my $lst=shift; open IN,$lst; my (@a,@b); my $len; my $max=0; while(){ chomp; next if /^##/; # 上面这行代码表示以##开头的行被忽略 @b=split/\t/,$_; $len=@b; $max=$max > $len ? $max:$len; push @a,[@b]; } close IN; for my $i(0..$max-1){ for(@a){ @$_[$i]||=""; print "@$_[$i]\t"; } print "\n"; #换行显示 } 执行脚本： perl s1.pl ./data/xxx.vcf > ./data/xxx.vcf_1 筛选指定样本 转置之后的vcf_1文件每一行是一个样本，根据所需的样本待提取编号进行提取，提取所用到的perl脚本（s2.pl）如下： #! /usr/bin/perl -w use strict; open IN1,"$ARGV[0]" or die; open IN2,"$ARGV[1]" or die; open OUT,">$ARGV[2]" or die; #需要三个参数，需要提取的样本ID编号序列信息、待提取的原始文件名、提取后新生成的文件名 my %hash; while(){ chomp; my @a=split/\s/,$_; $hash{$a[0]}=1; } close IN1; while(){ chomp; print OUT $_,"\n" if /^#/; print OUT $_,"\n" if /POS/; print OUT $_,"\n" if /ID/; print OUT $_,"\n" if /REF/; print OUT $_,"\n" if /ALT/; print OUT $_,"\n" if /QUAL/; print OUT $_,"\n" if /FILTER/; print OUT $_,"\n" if /INFO/; print OUT $_,"\n" if /FORMAT/; #如果改行一以上述关键字出现，则全部保留，因为这些是文件的有用信息 my @b=split/\s/,$_; if(exists $hash{$b[0]}){print OUT join("\t",@b),"\n"} else{next} } close IN2;close OUT; 执行脚本： perl s2.pl list.txt ./data/xxx.vcf_1 ./data/xxx.vcf_2 重新转置数据 筛选完成后，需要对数据进行转置，还原为之前的数据格式，直接调用s1.pl即可完成。 perl s1.pl ./data/xxx.vcf_2 > ./data/xxx.vcf_3 转置后需要将刚开始提取的表头注释信息再添加上去，用cat命令将两个文件合并到一起。 cat zhushi.vcf ./data/xxx.vcf_3 > ./data/xxx.vcf_4 如上得到的vcf_4文件就是筛选得出的最终结果，对其改名，然后使用gzip压缩为原来的压缩包格式。流程结束！ mv ./data/*.vcf_4 `sed 's/.vcf_4/.vcf/g'` gzip *.vcf 批量进行操作 以下内容为bash脚本，用于批量处理多个vcf文件。 批量解压、转置、提取 #!/bin/bash gunzip ./data/*.gz for i in $(ls ./data/*.vcf) do echo $i perl s1.pl $i > ${i}_1 perl s2.pl list.txt ${i}_1 ${i}_2 perl s1.pl ${i}_2 > ${i}_3 cat zhushi.vcf ./data/xxx.vcf_3 > ./data/xxx.vcf_4 echo "ok" done 批量改文件名 #!/bin/bash for i in `ls ./data/*.vcf_4`;do echo $i | mv $i `sed 's/.vcf_4/.vcf/g'` echo ${i}ok done 本文由 mdnice 多平台发布 （编辑：财气旺网 - 财气网） 【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容! |
站长推荐