从SNP_VCF文件提取SNV

前期了解

VCF是测序文件的一种格式，详细记录测序样本的SNP信息，有很多介绍vcf文件内容的，在这里就不详细描述了。以下的工作主要对人类样本的SNP过滤出SNV的操作进行描述。

我收到的样本是这样子的：

共有24个，分别是24条染色体的.vcf.gz文件，每一条染色体的vcf文件都包含多个人的样本。

文件分割(样本少可以不做)

VCF测序文件可能包含有多个样本，所以需要执行bcftools脚本进行分割，首先把每一条染色体的单独文件放入到独立文件夹之中，类似于这样：

我们查看一下第一个文件夹的内容：

接下来就需要把这个文件分开成每个个体的1号染色体vcf文件：

bcftools view --threads n -s Sample_name WGS_1.vcf.gz -Oz -o Sample_name.vcf.gz

你需要在这个脚本里更改的内容有：
n：使用线程数量
Sample_name: vcf所包含的一个样本的名称
WGS_1.vcf.gz ：你输入的文件名称
Sample_name.vcf.gz ：你要输出的文件名称

这样你就可以在一号染色体的文件夹内看到被分开的多个样本信息。

数据注释

我使用的注释软件是Annovar，具体的使用详情可以另行查找，基于筛选的注释。
对一个样本的一条染色体的注释如下：

perl Your_PATH_1/table_annovar.pl -thread 1 Sample_name.vcf.gz Your_PATH_1/annovar/humandb/ -buildver hg19 -out Your_PATH_2/Sample_name -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a,EUR.sites.2015_08 -operation g,r,f,f,f,f -nastring . -vcfinput

Your_PATH_1：Annovar的路径
Your_PATH_2：输出的路径在哪里
refGene,cytoBand,exac03,avsnp147,dbnsfp30a,EUR.sites.2015_08:
这几个都是提前在Annovar下载好的数据库，g,r,f,f,f,f是注释的类型，具体需要看一下Annovar的使用方法。

文件整合(样本少可以不做)

注释的结果很简单，对每一个样本的一条染色体会生成这样的三个文件：

我们需要.txt结尾的文件整合到一个文件夹里面：
sample.list存放所有的样本名称，chrom.list存放1~24条染色体的编号。

cat sample.list |
  while read input; 
    do
      mkdir -p single_person/${input}
      cat chrom.list | 
        while read reads; 
          do 
             cp -p ${reads}/out/${input}_${reads}.hg19_multianno.txt single_person/${input}
          done
    done

注意找好路径，最终这些文件会被存放在一个叫做single_person的文件夹里面，每一个样本的所有染色体存放为一个单独的文件夹。

注释的过滤

我们的注释是全部的注释，但是有很多是不可靠的，需要按照一定的标准进行过滤。
过滤执行的是awk脚本，可以进行联合过滤。
批量执行脚本如下：

cat final.list |
  while read input; 
    do
        awk ' ($21<=0.05 && ($22 ~/D|\./) &&  $23>=0.453 &&  ($24 ~ /D|P|\./) &&  $25>=0.447 &&  ($26 ~ /D|P|\./) && ($28 ~ /D|\./) &&  ($30 ~ /A|D|\./) && ($32 ~ /H|M|\./) && ($34 ~ /D|\./) \
        && ($36 ~/D|\./) &&  $39>=20 &&  ($42 ~/D|\./) && ($44 ~/D|\./) &&  ($46 ~/D|\./)) { print $0} ' ${input}/*.hg19_multianno.txt >> PATH/${input}.txt
    done

$21指的是第21列，awk的正则化表达可以上网搜一下，对每一个样本的全套染色体进行过滤操作，最终得到一个样本一个txt文件。打开可以看到每一个样本过滤后的信息。