从 Fastq 文件中提取特定信息以进行测序分析

Question

我的目标是从基因组测序 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 列中还有一堆关于 read 的其他信息，其中每个部分由分号分隔。

我在这里关心的两个数字是：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

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

Answer 1

这需要一些研究来找出如何在 awk 中进行某种回顾。通过a thread in google groups!

发现很有意思

想法是使用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= FS  FS  FS 
            }' 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= FS  FS  FS }' a
chr1 115227813 C G 0.748387
chr1 115227814 G A,C,T 0.001415

从 Fastq 文件中提取特定信息以进行测序分析

Extracting specific information from a Fastq file for Sequencing Analysis

bash

awk

pattern-matching

sequencing

fastq

测试