Извлечение конкретной информации из файла Fastq для анализа последовательности

Моя цель — извлечь фрагменты данных из файлов 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.

Надеюсь, это прояснит вывод, который я ищу.


bash fastq awk pattern-matching sequencing
person The Nightman    schedule 04.04.2015    source источник
comment
желаемый результат обычно является лучшим объяснением   -  person fedorqui 'SO stop harming'    schedule 05.04.2015
comment
Нет проблем, я могу сделать это немного яснее.   -  person The Nightman    schedule 05.04.2015
comment
пожалуйста, покажите вывод из вашей 2-й строки данных, вам нужны все A, T, C из поля 5? Удачи.   -  person shellter    schedule 05.04.2015
comment
нет проблем, и да, я хочу все из поля 5.   -  person The Nightman    schedule 05.04.2015

Ответы (1)


arrow_upward
2
arrow_downward

Это потребовало некоторых исследований, чтобы выяснить, как сделать своего рода просмотр в awk. Было интересно узнать об этом через тред в группах google< /а>!

Идея состоит в том, чтобы использовать 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
person fedorqui 'SO stop harming'    schedule 04.04.2015
comment
Красиво, это именно то, что я ищу. Федорки большое спасибо. - person The Nightman; 05.04.2015