0

我的目标是从基因组测序 Fastq 文件中提取数据片段并绘制它们。我想获得每个测序读数的识别信息，然后是关于读数的两条信息。

下面我粘贴了两个从 Fastq 文件读取的内容以供参考，如果有帮助的话。

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  12_S12_L001
chr1    115227813       .       C       G       2120.73 .       AB=0.725;ABP=73.366;AC=1;AF=0.5;AN=2;AO=116;CIGAR=1X;DP=160;DPB=160;DPRA=0;EPP=254.901;EPPR=87.6977;GTI=0;LEN=1;MEANALT=3;MQM=60;MQMR=60;NS=
1;NUMALT=1;ODDS=152.168;PAIRED=0.991379;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=3761;QR=1366;RO=39;RPP=254.901;RPPR=87.6977;RUN=1;SAF=116;SAP=254.901;SAR=0;SRF=39;SRP=87.6977;SRR=0;TYPE=snp  GT:DP:RO:QR:AO:Q
A:GL    0/1:160:39:1366:116:3761:-10,0,-10
chr1    115227814       .       G       A,C,T   8.27007e-12     .       AB=0,0,0;ABP=0,0,0;AC=0,0,0;AF=0,0,0;AN=2;AO=120,11,35;CIGAR=1X,1X,1X;DP=84826;DPB=84826;DPRA=0,0,0;EPP=263.587,26.8965,79.0118;EPPR
=183840;GTI=0;LEN=1,1,1;MEANALT=3,3,3;MQM=60,60,60;MQMR=59.9996;NS=1;NUMALT=3;ODDS=115105;PAIRED=1,1,1;PAIREDR=0.990917;PAO=0,0,0;PQA=0,0,0;PQR=0;PRO=0;QA=4206,292,1061;QR=2822527;RO=84660;RPP=263.587,26.
8965,79.0118;RPPR=183840;RUN=1,1,1;SAF=120,11,35;SAP=263.587,26.8965,79.0118;SAR=0,0,0;SRF=84660;SRP=183840;SRR=0;TYPE=snp,snp,snp      GT:DP:RO:QR:AO:QA:GL    0/0:84826:84660:2822527:120,11,35:4206,292,1
061:0,-10,-10,-10,-10,-10,-10,-10,-10,-10

上面，您可以看到每个读取都以进行读取的染色体编号开始，以及读取在第 1 列和第 2 列中该染色体上的位置。在第 4 列中有参考碱基对，第 5 列包含变体读。然后在第 8 列中还有一堆关于读取的其他信息，其中每个部分用分号分隔。

我在这里关心的两个数字是：RO=和AO=。

我想创建一个仅包含 1、2、4、5 列信息的输出文件，然后将 AO/RO 的分数放入最后一列。

作为从第一行开始的输出示例，我想要以下输出：

chr1    115227813    C    G    0.74838
chr1    115227814    G    A,C,T    0.00142

其中 0.74838 由 RO=39 和 AO=116 计算得出，因此 116/(39+116)=0.74838。并且由 RO=84660 和 AO=120 计算，因此 120/(84660+120)=0.00142

希望这可以澄清我正在寻找的输出。

4

1 回答 1

2

这需要一些研究来找出如何在awk. 通过谷歌群组中的线程发现它很有趣！

这个想法是用来gensub()获取variable=value给定行中的，然后将其打印回来，删除该行的其余内容。所以如果我们有hello hello;AO=23;bla bla bla，我们就得到23。

awk 'v {
         ro=gensub(/^.*;RO=([0-9]*).*$/, "\\1", "1"); 
         printf "%s %f\n", f, (ao/(ao + ro)); v=0
     }
     /^chr/ {ao=gensub(/^.*;AO=([0-9]*).*$/,"\\1", "1");
             v=1;
             f=$1 FS $2 FS $4 FS $5
            }' file

基本上，我们寻找以 . 开头的行chr。在这些中，我们捕获了第 1、第 2、第 4 和第 5 个值。然后，我们捕获旁边的任何内容AO=（仅数字）。由于RO=出现在下一行，我们设置了一个标志，以便在读取下一行时搜索它。然后，我们获取该值并打印完整的数据集。最后我们取消设置标志，所以我们再次开始循环。

测试

$ awk 'v {ro=gensub(/^.*;RO=([0-9]*).*$/, "\\1", "1"); printf "%s %f\n", f, (ao/(ao + ro)); v=0} /^chr/ {ao=gensub(/^.*;AO=([0-9]*).*$/,"\\1", "1"); v=1; f=$1 FS $2 FS $4 FS $5}' a
chr1 115227813 C G 0.748387
chr1 115227814 G A,C,T 0.001415
于 2015-04-04T23:05:11.593 回答