fastq 文件中的修剪序列和质量答案

【问题标题】：trim sequences and quality in fastq filefastq 文件中的修剪序列和质量
【发布时间】：2015-03-20 20:53:45
【问题描述】：

我在目录中有一堆 fastq 文件，我想将序列修剪 2 个核苷酸和质量（如果读取有 51 个碱基对并且以 CTG 或 TTG 结尾）。

这是我编写的 shell 脚本，但我遇到了一些错误，需要帮助，因为我是 shell 脚本的新手

输入：

@HWI-ST1072:187:C35YUACXX:7:1101:1609:1983 1:N:0:ACAGTG
NGGAGAAAGAGAGTGTGTTTTTAGGGGGAGATTTTTAAAATGGTTGTTTTG
+
#0<BFFFFFFFFF<BFFFIIFFFFFIIIBFFFFFIIFIIIIIFFBFFFFFF
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTATTCGGGAGGTTGAGCTG
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFIIIFIIFIIFFFFIIFF
@HWI-ST1072:187:C35YUACXX:7:1101:9351:2210 1:N:0:ACAGTG
CGGTTTTGTTTTATTTTGTATGATTAGGAGGGTTTTGGAGGTTTAGTTACC
+
BBBFFFFFFFFFFIIIIIFFIIFIIIIIIIIIFFIIFIFIIFFIIIFIIII
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTAT
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFI

输出：

@HWI-ST1072:187:C35YUACXX:7:1101:1609:1983 1:N:0:ACAGTG
NGGAGAAAGAGAGTGTGTTTTTAGGGGGAGATTTTTAAAATGGTTGTTT
+
#0<BFFFFFFFFF<BFFFIIFFFFFIIIBFFFFFIIFIIIIIFFBFFFF
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTATTCGGGAGGTTGAGC
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFIIIFIIFIIFFFFII
@HWI-ST1072:187:C35YUACXX:7:1101:9351:2210 1:N:0:ACAGTG
CGGTTTTGTTTTATTTTGTATGATTAGGAGGGTTTTGGAGGTTTAGTTACC
+
BBBFFFFFFFFFFIIIIIFFIIFIIIIIIIIIFFIIFIFIIFFIIIFIIII
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTAT
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFI

脚本：

for sample in *.fastq;do
    name=$(echo ${sample} | sed 's/.fastq//')
    while read line;do
        if [ ${line:0:1} == "@" ] ; then
                head="${line}"
                $echo $head
        elif [ "${head}" ] && [ "${line}" ] ; then
                length=${#line}
                if [ "${length}" = 51 -a "${line}" =~ *CTG|*TTG ] ; then
                        sequence= substr($line,0,49)
                        #echo $sequence
                fi
        elif [ ${line:0:1} == "+" ] ; then
                plus="${line}"
                #echo $plus
        elif [ "${plus}" ] && [ "${line}" ] ; then
                quality= substr($line,0,49)
                #echo $quality
        fi
        print "${head}\n${sequence}\n${plus}\n${quality}" > ${name}_new.fq
   done < $sample
done

【问题讨论】：

我在创建 substr 时出错！有没有办法可以拆分行并保存在变量中
shell 是一个调用工具的环境。它具有编程语言结构，可让您对这些调用进行排序。 awk 是处理文本文件的 UNIX 命令。因此，到目前为止您所做的完全是错误的方法 - 在 shell 中执行此操作的方法是编写一个 awk 脚本来解析您的文本文件，然后从 shell 调用它。

标签： bash shell unix awk fastq

【解决方案1】：

不要 100% 了解您在做什么，但解决了一些问题。试试下面

#!/bin/bash
for sample in *.fastq; do
  name="${sample/.fastq/}"
  while read -r line; do
    if [[ $line == '@'* ]]; then
      head="$line" && echo "$head" >> "${name}_new.fq"
    elif [[ -n $head && ${#line} == 51 && $line =~ (CTG|TTG)$ ]]; then
      sequence="${line:0:49}" && echo "$sequence" >> "${name}_new.fq"
    elif [[ $line == '+'* ]]; then
      plus="$line" && echo "$line" >> "${name}_new.fq"
    else
      quality="$line" && echo "$quality" >> "${name}_new.fq"
    fi
  done < "$sample"
done

示例输出

> cat sample_new.fq

> cat sample.fastq
@HWI-ST1072:187:C35YUACXX:7:1101:1609:1983 1:N:0:ACAGTG
NGGAGAAAGAGAGTGTGTTTTTAGGGGGAGATTTTTAAAATGGTTGTTTTG
+
#0<BFFFFFFFFF<BFFFIIFFFFFIIIBFFFFFIIFIIIIIFFBFFFFFF
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTATTCGGGAGGTTGAGCTG
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFIIIFIIFIIFFFFIIFF
@HWI-ST1072:187:C35YUACXX:7:1101:9351:2210 1:N:0:ACAGTG
CGGTTTTGTTTTATTTTGTATGATTAGGAGGGTTTTGGAGGTTTAGTTACC
+
BBBFFFFFFFFFFIIIIIFFIIFIIIIIIIIIFFIIFIFIIFFIIIFIIII
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTAT
+
#0<BFFFFFFFFFFIIBFFIIIIIIFIIIFFIIFI

> ./abovescript

> cat sample_new.fq
@HWI-ST1072:187:C35YUACXX:7:1101:1609:1983 1:N:0:ACAGTG
NGGAGAAAGAGAGTGTGTTTTTAGGGGGAGATTTTTAAAATGGTTGTTT
+
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTATTCGGGAGGTTGAGC
+
@HWI-ST1072:187:C35YUACXX:7:1101:9351:2210 1:N:0:ACAGTG
CGGTTTTGTTTTATTTTGTATGATTAGGAGGGTTTTGGAGGTTTAGTTACC
+
BBBFFFFFFFFFFIIIIIFFIIFIIIIIIIIIFFIIFIFIIFFIIIFIIII
@HWI-ST1072:187:C35YUACXX:7:1101:1747:1995 1:N:0:ACAGTG
NGGTTGTGGTGGTGGGTATTTGTAGTTTTATTTAT
+

【讨论】：

它删除序列的每 2 个核苷酸！！！但我只想在序列以 CTG 或 TTG @BroSlow 结尾时删除 2 个核苷酸
@user2243831 我想我不太明白。如果一行以# 开头，你想做什么？例如查看更新不是 51 个字符且与其他参数匹配的行（例如，刚刚打印以 # 开头的行）。
如果序列有 2 个条件（长度应为 51 并且末尾有 CTG 或 TTG），我只想将序列修剪 2 个核苷酸。可能还有一些其他序列是 51，但如果它们没有有 CTG 或 TTG 我不修剪它们。# 行也应该根据条件删除@BroSlow
@user2243831 再试一次。如果不是这样，您需要更新预期输出。
我只需要根据条件修改序列的第二行，我可以打印从 0 到 49 的第 4 行！我们需要 substr 函数@BroSlow