linux - 在 Linux 中替换 FastQ 文件中的特定核苷酸

Question

我有一些 fastaq 文件需要分析。主要问题是我目前使用的分析工具只接受 ACTG 作为核苷酸，而不接受 IUPAC 代码中的其余命名法(R、W 等)。

我编写了这段代码来改变特定的核苷酸:

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

它正在运行，但速度非常慢。你知道更有效的方法吗？

非常感谢!

Answer 1

假设你有fastq 格式，我建议使用专门的库，biopython 或bioperl 是不错的选择。

cat example.fastq

@IDAGTCGTACTGGACTGYGCSAACTG+IIIIIIIIIIIIIIIIIIIIIIII@ID2RWMVYDSKNAAAAAAAAAAAAAAA+IIIIIIIIIIIIIIIIIIIIIIII

However, solution using awk

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

你明白了，

@IDAGTCGTACTGGACTGCGCGAACTG+IIIIIIIIIIIIIIIIIIIIIIII@ID2AAAACCGGGAAAAAAAAAAAAAAA+IIIIIIIIIIIIIIIIIIIIIIII