2017-02-17 65 views
2

我有一些我需要分析的fastaq文件。主要问题是我目前使用的分析工具只接受ACTG作为核苷酸,而不接受IUPAC代码(R,W等)中的其余术语。替换Linux中的FastaQ文件中的特定核苷酸

我做了这个代码更改的特定核苷酸:

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

'为(N a)中的{核苷酸= a [n]};'工作不好 –

+0

您的预期产出是?和示例输入? –

+0

你对末尾的变量'str'没有任何作用 –

回答

3

对于这个假设你有fastq格式,我建议使用专门的库,biopythonbioperl是不错的选择。

猫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

啊!现在我想起了! – Inian