vcf文件需要用perl處理取出DP>3的點(diǎn)

Tracy_Panda

全基因組數(shù)據(jù)，經(jīng)過SNP calling產(chǎn)生了很大的vcf文件，但是很多點(diǎn)實(shí)際上質(zhì)量很差沒有計(jì)算的意義，所以想把DP>1的數(shù)據(jù)篩出來進(jìn)行後續(xù)計(jì)算，請(qǐng)教大家，不知道如何處理比較好，因爲(wèi)是全基因組數(shù)據(jù)，文件過大利用excel根本打不開，所以excel的方法暫時(shí)不考慮了。文件格式如下：
#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列下的第一個(gè)參數(shù)；
因爲(wèi)文件很大，所以無法上傳，（壓縮之後還是不行），所以只節(jié)選了文件中的一小部分，希望大家能幫個(gè)忙，謝謝了先�。。。。�

sunzhiguolu

文件有多大?

sunzhiguolu

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

復(fù)制代碼

Tracy_Panda

回復(fù) 2# sunzhiguolu 50G

Tracy_Panda

回復(fù) 3# sunzhiguolu

額，不好意思，結(jié)果不對(duì)，顯示的是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

回復(fù) 5# Tracy_Panda
50GB 的話, 用一般的工具軟件的確打開就 over 了。

sunzhiguolu

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

Tracy_Panda

回復(fù) 7# sunzhiguolu

D:\復(fù)旦\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.

不好意思，結(jié)果還是如此。。。因爲(wèi)樣本特殊，所以很多沒用的信息會(huì)得到，但是也要過濾掉大部分?jǐn)?shù)據(jù)，所以我也在想起其他辦法能不能實(shí)現(xiàn)。。。

哈哈，如果可以的話可以接著幫忙想一下哪出錯(cuò)啦

sunzhiguolu

如果是在 Windows 端運(yùn)行的話,

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

復(fù)制代碼

sunzhiguolu

回復(fù) 8# Tracy_Panda

不好意思，結(jié)果還是如此。。。因爲(wèi)樣本特殊，所以很多沒用的信息會(huì)得到，但是也要過濾掉大部分?jǐn)?shù)據(jù)，所以我也在想起其他辦法能不能實(shí)現(xiàn)。。。

文件格式是否固定?  處理?xiàng)l件如何?  能否給出預(yù)期的結(jié)果?