文章目錄

• 前期了解
• 文件分割(樣本少可以不做)
• 數(shù)據(jù)注釋
• 文件整合(樣本少可以不做)
• 注釋的過濾

前期了解

VCF是測序文件的一種格式，詳細(xì)記錄測序樣本的SNP信息，有很多介紹vcf文件內(nèi)容的，在這里就不詳細(xì)描述了。以下的工作主要對人類樣本的SNP過濾出SNV的操作進(jìn)行描述。

我收到的樣本是這樣子的：

共有24個(gè)，分別是24條染色體的.vcf.gz文件，每一條染色體的vcf文件都包含多個(gè)人的樣本。

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

VCF測序文件可能包含有多個(gè)樣本，所以需要執(zhí)行bcftools腳本進(jìn)行分割，首先把每一條染色體的單獨(dú)文件放入到獨(dú)立文件夾之中，類似于這樣：

我們查看一下第一個(gè)文件夾的內(nèi)容：

接下來就需要把這個(gè)文件分開成每個(gè)個(gè)體的1號(hào)染色體vcf文件：

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

你需要在這個(gè)腳本里更改的內(nèi)容有：
n：使用線程數(shù)量
Sample_name: vcf所包含的一個(gè)樣本的名稱
WGS_1.vcf.gz ：你輸入的文件名稱
Sample_name.vcf.gz ：你要輸出的文件名稱

這樣你就可以在一號(hào)染色體的文件夾內(nèi)看到被分開的多個(gè)樣本信息。

數(shù)據(jù)注釋

我使用的注釋軟件是Annovar，具體的使用詳情可以另行查找，基于篩選的注釋。
對一個(gè)樣本的一條染色體的注釋如下：

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:
這幾個(gè)都是提前在Annovar下載好的數(shù)據(jù)庫，g,r,f,f,f,f是注釋的類型，具體需要看一下Annovar的使用方法。

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

注釋的結(jié)果很簡單，對每一個(gè)樣本的一條染色體會(huì)生成這樣的三個(gè)文件：

我們需要.txt結(jié)尾的文件整合到一個(gè)文件夾里面：
sample.list存放所有的樣本名稱，chrom.list存放1~24條染色體的編號(hào)。

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

注意找好路徑，最終這些文件會(huì)被存放在一個(gè)叫做single_person的文件夾里面，每一個(gè)樣本的所有染色體存放為一個(gè)單獨(dú)的文件夾。

注釋的過濾

我們的注釋是全部的注釋，但是有很多是不可靠的，需要按照一定的標(biāo)準(zhǔn)進(jìn)行過濾。
過濾執(zhí)行的是awk腳本，可以進(jìn)行聯(lián)合過濾。
批量執(zhí)行腳本如下：

cat final.list |while read input; doawk ' ($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}.txtdone

$21指的是第21列，awk的正則化表達(dá)可以上網(wǎng)搜一下，對每一個(gè)樣本的全套染色體進(jìn)行過濾操作，最終得到一個(gè)樣本一個(gè)txt文件。打開可以看到每一個(gè)樣本過濾后的信息。

總結(jié)

以上是生活随笔為你收集整理的从SNP_VCF文件提取SNV的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。