1

出力ファイルには、3 つのサンプルに対して 800,000 行と 8 つのフィールドがあります。ここで2行を抽出するだけです。chr、位置、SNP-ID、品質、DP、QD、遺伝子型 (./.、0/0、/0/1、または 1/1) など、各行の特定の情報のみを抽出したい。これらの情報を抽出して新しいファイルを作成するスクリプトが必要です: アドバイスをお願いします。ありがとう

#chr pos SNP-ID Qual Info geno(sample1) geno(sample2) geno(sample3)
chrM 152 rs117135796 7427.14 AC=2;AF=0.333;AN=6;BaseQRankSum=-20.485;DB;DP=702;DS;Dels=0.00;FS=167.659;HaplotypeScore=2.6106;MLEAC=2;MLEAF=0.333;MQ=50.00;MQ0=0;MQRankSum=-1.507;QD=36.77;ReadPosRankSum=12.041 0/0:250,0:237:99:0,701,10320 0/0:250,0:238:99:0,713,10507 1/1:0,202:192:99:7465,572,0
chr10 5874 rs118203891 33.13 AC=1;AF=0.167;AN=6;BaseQRankSum=1.454;DB;DP=657;DS;Dels=0.00;FS=124.424;HaplotypeScore=5.1214;MLEAC=1;MLEAF=0.167;MQ=45.31;MQ0=0;MQRankSum=2.462;QD=0.15;ReadPosRankSum=-8.096 0/1:204,24:206:64:64,0,6345 0/0:203,0:193:99:0,473,6944 0/0:226,0:215:99:0,524,6448
4

2 に答える 2

2

試す：

awk -f ext.awk data.txt > summary.txt

data.txt入力データファイルはどこにあり、次のext.awkとおりです。

NR>1 {
    match($5,/(DP=[^;]+);/,a)
    DP=a[1]
    match($5,/(QD=[^;]+);/,a)
    QD=a[1]
    match($6,/^([^:]+\/[^:]+):/,a)
    gt1=a[1]
    match($7,/^([^:]+\/[^:]+):/,a)
    gt2=a[1]
    match($8,/^([^:]+\/[^:]+):/,a)
    gt3=a[1]
    print $1,$2,$3,$4,DP,QD,gt1,gt2,gt3
}

アップデート

遺伝子型が各フィールドの最初の 3 文字 ($6 から $NF) によって与えられると仮定すると、以下を試すことができます。

NR>1 {
    match($5,/(DP=[^;]+);/,a)
    DP=a[1]
    match($5,/(MQ=[^;]+);/,a)
    MQ=a[1]
    printf "%s %s %s %s %s %s ", $1,$2,$3,$4,DP,MQ
    for (i=6; i<=NF; i++) {
        printf "%s", substr($i,1,3)
        if (i<NF) printf " "
        else printf "\n"
    }
}

アップデート

あなたがしたい場合は：

  • DP<10 または MQ<50 の場合は、その行を削除します。
  • 次のように遺伝子型を変換します: (NA, 0, 1, 2)
  • 変換 。/。「な」に、
  • 0/0 を「0」に変換
  • 0/1 を「1」に変換
  • 1/1 を "2" に変換する

次に、試すことができます：

BEGIN {
    geno["./."]="NA"
    geno["0/0"]="0"
    geno["0/1"]="1"
    geno["1/1"]="2"
}
NR>1 {
    match($5,/(DP=[^;]+);/,a)
    DP=a[1]
    match(DP,/=(.*)$/,a)
    dpv=a[1]
    match($5,/(MQ=[^;]+);/,a)
    MQ=a[1]
    match(MQ,/=(.*)$/,a)
    mqv=a[1]
    if (dpv<10 || mqv<50) next
    else {
        printf "%s %s %s %s %s %s ", $1,$2,$3,$4,DP,MQ
        for (i=6; i<=NF; i++) {
            type=substr($i,1,3)
            printf "%s", geno[type]
            if (i<NF) printf " "
            else printf "\n"
        }
    }
}
于 2013-10-11T22:23:42.400 に答える
1

Perl はすてきな簡潔なプログラムを提供します:

perl -ane '
    BEGIN {$, = " "}
    @fields = @F[0..3];
    push @fields, $1, $2 if $F[4] =~ /(DP=.+?);.*(QD=.+?);/;
    push @fields, (split /:/)[0] for @F[5,6,7];
    print @fields, "\n";
' <<END
chrM 152 rs117135796 7427.14 AC=2;AF=0.333;AN=6;BaseQRankSum=-20.485;DB;DP=702;DS;Dels=0.00;FS=167.659;HaplotypeScore=2.6106;MLEAC=2;MLEAF=0.333;MQ=50.00;MQ0=0;MQRankSum=-1.507;QD=36.77;ReadPosRankSum=12.041 0/0:250,0:237:99:0,701,10320 0/0:250,0:238:99:0,713,10507 1/1:0,202:192:99:7465,572,0
chr10 5874 rs118203891 33.13 AC=1;AF=0.167;AN=6;BaseQRankSum=1.454;DB;DP=657;DS;Dels=0.00;FS=124.424;HaplotypeScore=5.1214;MLEAC=1;MLEAF=0.167;MQ=45.31;MQ0=0;MQRankSum=2.462;QD=0.15;ReadPosRankSum=-8.096 0/1:204,24:206:64:64,0,6345 0/0:203,0:193:99:0,473,6944 0/0:226,0:215:99:0,524,6448
END

chrM 152 rs117135796 7427.14 DP=702 QD=36.77 0/0 0/0 1/1 
chr10 5874 rs118203891 33.13 DP=657 QD=0.15 0/1 0/0 0/0
于 2013-10-11T22:45:21.700 に答える