生信分析笔记

V1

2022/11/18阅读:26主题:山吹

Linux环境下利用perl脚本批量筛选VCF文件指定样本

今天分享一个实用的数据提取思路:从VCF文件中提取指定列,比如有1000个样本,需要从中提取出200个制定样本,形成新的vcf文件。需要在Linux环境下进行操作,主要用到perl脚本和shell脚本。

原始数据

  1. 原始文件vcf :开始为若干行注释信息以##开头,然后是表头信息由#开始(该行包括样本ID)

  2. 待提取的样本ID:txt文件,每个样本ID一行(list.txt)

什么是VCF文件?

伴随着大规模的基因分型及测序数据的产生,迫切需要一种新的格式来记录高效的记录这些变异信息。VCF(Variant Call Format)就是这样一种用来贮存基因序列变异信息的文本文件。

VCF 格式文件信息:

  1. header section:元信息(meta-information),以‘##’为前缀,通常包含fileformat、fileDate、reference等信息。头行信息( header line ),以‘#’为前缀
  2. data section:该部分为主题部分,记录了每个样品每个位点处的基因分型信息。可以看做是一个大的Excel表格,每行是一个位点,每列是一个样本。

创建一个工作文件夹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文件的每一列为一个样本,为了进行筛选,先将其转置,变为每行一个样本的格式,使用如下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(<IN>){
        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(<IN1>){
chomp;
my @a=split/\s/,$_;
$hash{$a[0]}=1;
}
close IN1;
while(<IN2>){
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

分类:

后端

标签:

后端

作者介绍

生信分析笔记
V1

欢迎关注公众号:生信分析笔记