vcf文件需要用perl处理取出DP>3的点

Tracy_Panda

全基因組數據，經過SNP calling產生了很大的vcf文件，但是很多點實際上質量很差沒有計算的意義，所以想把DP>1的數據篩出來進行後續計算，請教大家，不知道如何處理比較好，因爲是全基因組數據，文件過大利用excel根本打不開，所以excel的方法暫時不考慮了。文件格式如下：
#CHROM        POS        ID        REF        ALT        QUAL        FILTER        INFO        FORMAT        M69B.q30.srt.rmp.bam
chr1        824423        .        G        A        7.8        .        DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=37;FQ=-30        GTL:GQ        1/1:37,3,0:4
chr1        824425        .        G        A        7.8        .        DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=37;FQ=-30        GTL:GQ        1/1:37,3,0:4
chr1        837554        .        G        A        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        838570        .        G        A        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        840821        .        G        A        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        846005        .        C        T        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        847156        .        C        T        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        847807        .        C        T        7.8        .        DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=37;FQ=-30        GTL:GQ        1/1:37,3,0:4
chr1        850705        .        G        A        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        854250        .        A        G        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GTL:GQ        1/1:74,6,0:10
chr1        856420        .        G        A        4.77        .        DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=33;FQ=-30        GT:PL:GQ        0/1:33,3,0:3
chr1        856427        .        G        A        4.13        .        DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=33;FQ=-30        GT:PL:GQ        0/1:32,3,0:3
chr1        857025        .        G        T        42.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-33        GT:PL:GQ        1/1:74,6,0:10
chr1        860251        .        G        A        22.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=32;FQ=-33        GT:PL:GQ        1/1:54,6,0:10
chr1        860252        .        G        A        33.8        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=32;FQ=-33        GT:PL:GQ        1/1:65,6,0:10
chr1        860273        .        G        A        15.9        .        DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=32;FQ=-33        GT:PL:GQ        1/1:47,6,0:10
chr1        860288        .        CG        CGG        29.2        .        INDEL;IS=2,1.000000;DP=2;VDB=5.960000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=37;FQ=-40.5        GT:PL:GQ        1/1:68,6,0:10



DP是INFO列下的第一個參數；
因爲文件很大，所以無法上傳，（壓縮之後還是不行），所以只節選了文件中的一小部分，希望大家能幫個忙，謝謝了先！！！！！

sunzhiguolu

文件有多大?

sunzhiguolu

1. perl -nle 'next if(/\A\s*#/);print if(/\bdp=(\d+)/i and $1>1)' f

复制代码

Tracy_Panda

回复 2# sunzhiguolu 50G

Tracy_Panda

回复 3# sunzhiguolu

額，不好意思，結果不對，顯示的是Can't find string terminator "'" anywhere before EOF at -e line 1.
我輸入的是perl -nle 'next if(/\A\s*#/);print if(/\bdp=(\d+)/i and $1>3)' f test.vcf > test2.vcf。

不知道哪出了一些問題。。

sunzhiguolu

回复 5# Tracy_Panda
50GB 的话, 用一般的工具软件的确打开就 over 了。

sunzhiguolu

回复 5# Tracy_Panda
perl -nle 'next if(/\A\s*#/);print if(/\bdp=(\d+)/i and $1>3)' test.vcf > test2.vcf

Tracy_Panda

回复 7# sunzhiguolu

D:\复旦\KC>perl -nle 'next if(/\A\s*#/);print if(/\bdp=(\d+)/i and $1>3)' test.vcf > test2.vcf
Can't find string terminator "'" anywhere before EOF at -e line 1.

不好意思，結果還是如此。。。因爲樣本特殊，所以很多沒用的信息會得到，但是也要過濾掉大部分數據，所以我也在想起其他辦法能不能實現。。。

哈哈，如果可以的話可以接着幫忙想一下哪出錯啦

sunzhiguolu

如果是在 Windows 端运行的话,

1. perl -nle "next if(/\A\s*#/);print if(/\bdp=(\d+)/i and $1>3)" test.vcf > test2.vcf

复制代码

sunzhiguolu

回复 8# Tracy_Panda

不好意思，結果還是如此。。。因爲樣本特殊，所以很多沒用的信息會得到，但是也要過濾掉大部分數據，所以我也在想起其他辦法能不能實現。。。

文件格式是否固定?  处理条件如何?  能否给出预期的结果?