обрезать последовательности и качество в файле fastq

У меня есть куча файлов fastq в каталоге, и я хочу обрезать последовательность на 2 нуклеотида и качество (если чтение имеет 51 пару оснований и заканчивается CTG или TTG).

вот что я написал как сценарий оболочки, но я получаю некоторые ошибки, нужна помощь, так как я новичок в сценариях оболочки

Вход:

@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

bash shell fastq unix awk
person user2243831    schedule 06.02.2014    source источник
comment
У меня ошибка при создании подстроки! Есть ли способ разделить строки и сохранить в переменной   -  person user2243831    schedule 06.02.2014
comment
оболочка - это среда, из которой можно вызывать инструменты. Он имеет конструкции языка программирования, позволяющие упорядочивать эти вызовы. awk — это команда UNIX для обработки текстовых файлов. Поэтому то, что вы сделали до сих пор, является совершенно неправильным подходом - способ сделать это в оболочке - написать awk-скрипт для анализа вашего текстового файла, а затем вызвать его из оболочки.   -  person Ed Morton    schedule 06.02.2014

Ответы (1)


arrow_upward
1
arrow_downward

Не на 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
+
person Reinstate Monica Please    schedule 06.02.2014
comment
удаление каждых 2 нуклеотидов последовательности!!! но я хочу удалить 2 нуклеотида, только если последовательность заканчивается CTG или TTG @BroSlow - person user2243831; 06.02.2014
comment
@user2243831 user2243831 Наверное, я не очень понимаю. Что вы хотите сделать, если строка начинается с #? например см. обновление, где строка, которая не состоит из 51 символа и соответствует другим параметрам (например, строка, начинающаяся с #, только что напечатана). - person Reinstate Monica Please; 06.02.2014
comment
я просто хочу обрезать последовательности на 2 нуклеотида, если они имеют 2 условия (длина должна быть 51 и иметь CTG или TTG в конце). могут быть некоторые другие последовательности, которые равны 51, но если они не имеют CTG или TTG, я не обрезаю их Строка .# также должна быть удалена в соответствии с условиями @BroSlow - person user2243831; 06.02.2014
comment
@ user2243831 Попробуйте еще раз. Если это не так, вам нужно обновить ожидаемый результат. - person Reinstate Monica Please; 06.02.2014
comment
мне нужна только вторая строка последовательности для изменения в зависимости от условий, 4 строки, которые я могу напечатать от 0 до 49! нам нужна функция substr @BroSlow - person user2243831; 06.02.2014
comment
я изменю ввод и вывод @BroSlow - person user2243831; 06.02.2014
comment
@ user2243831 No {line:0:49} уже получает 49-символьную подстроку из 0 с учетом ваших условий. - person Reinstate Monica Please; 06.02.2014
comment
теперь, если вы видите, что первые 2 последовательности имеют TTG и CTG в конце, и они равны 51, мне нужно обрезать их, как вы видите в выводе, а другие 2 последовательности разные, одна из них равна 51, но у них нет конечных последовательностей, мне нужно сохранить его и последняя последовательность не удовлетворяет обоих, я оставляю ее @BroSlow - person user2243831; 06.02.2014
comment
@user2243831 user2243831 Правильно ... приведенный выше скрипт уже делает это и получает точно такой же вывод, за исключением того, что с удаленной строкой #, что, как я думал, вы хотели? - person Reinstate Monica Please; 06.02.2014
comment
Я вижу, ваш скрипт обрезает каждую 51 строку последовательности в файле, даже если у них нет концов TTG или CTG @BroSlow - person user2243831; 06.02.2014
comment
@user2243831 user2243831 Нет, если в вашей среде не происходит что-то странное. Покажите вывод, который вы получаете из приведенного выше скрипта. - person Reinstate Monica Please; 06.02.2014
comment
Это то, что я получаю для последовательности трех строк CGGTTTTGTTTTATTTTGTATGATTAGGAGGGTTTTGGAGGTTTAGTTA - person user2243831; 06.02.2014
comment
@user2243831 user2243831 Вы уверены, что у вас есть именно тот сценарий, который указан выше? - person Reinstate Monica Please; 06.02.2014
comment
Мне не хватает # строк в выходных файлах, они мне нужны - person user2243831; 06.02.2014
comment
@HWI-ST1072:187:C35YUACXX:7:1101:1609:1983 1:N:0:ACAGTG NGGAGAAAGAGAGTGTGTTTTTAGGGGGAGATTTTAAATGGTTGTTT + #0‹BFFFFFFFFF‹BFFFIIFFFFFIIIBFFFFFIIFIIIIIFFBFFFFFF я не вижу изменения 4 строки на 49 - person user2243831; 06.02.2014