2017-02-17 27 views
2

У меня есть некоторые файлы fastaq, которые мне нужно проанализировать. Основная проблема заключается в том, что инструмент анализа, с которым я сейчас работаю, принимает только ACTG как нуклеотиды, а не остальные номенклатуры в коде IUPAC (R, W и т. Д.).Замещение определенных нуклеотидов в файлах FastaQ в Linux

Я сделал этот код, чтобы изменить конкретные нуклеотиды:

awk '{ 
    split($2,a,"") ; 
    str="" ; 
    for (n in a) {nucleotide=a[n]} ; 
    if (nucleotide~/[ACTG]/) {str=str""nucleotide} 
    else { 
     if (nucleotide~/[RWMV]/) {str=str""A} 
     else { 
      if (nucleotide~/[YD]/) {str=str""C} 
      else { 
       if (nucleotide~/[SKN]/) {str=str""G} 
       else {str=str""T} 
      } 
     } 
    } 
}' | head 

Это работает, но это супер медленно. Знаете ли вы более эффективный способ сделать это?

Большое вам спасибо!

+0

'для (п в) {нуклеотида = a [n]}; 'Не работает хорошо –

+0

Каков ваш ожидаемый результат? и пример ввода? –

+0

Вы ничего не делаете с переменной 'str' в конце –

ответ

3

Для этого при условии, что у вас есть формат fastq, я рекомендую использовать специализированную библиотеку, biopython или bioperl - хорошие варианты.

кошка example.fastq

 
@ID 
AGTCGTACTGGACTGYGCSAACTG 
+ 
IIIIIIIIIIIIIIIIIIIIIIII 
@ID2 
RWMVYDSKNAAAAAAAAAAAAAAA 
+ 
IIIIIIIIIIIIIIIIIIIIIIII 

Однако решение с использованием awk

awk 'NR%4==2{gsub(/[RWMV]/,"A"); gsub(/[YD]/,"C"); gsub(/[SKN]/,"G")}1' example.fastq 

вы получите,

 
@ID 
AGTCGTACTGGACTGCGCGAACTG 
+ 
IIIIIIIIIIIIIIIIIIIIIIII 
@ID2 
AAAACCGGGAAAAAAAAAAAAAAA 
+ 
IIIIIIIIIIIIIIIIIIIIIIII 
+0

Первоначально я думал, что использовал 'sed' ....' sed '/^@/{ n; y/RWMVYDSKN/AAAACCGGG /;} 'example.fastq '.... –

+0

@Inian' sub() 'только изменяет первое вхождение, оно не работает в этом случае –

+0

Ahh! Теперь я помню! – Inian