0

I have a file containing a list of DNA sequence names and another containing DNA sequences. They look like this:

$ cat list.txt
seq1
seq3

$ cat sequences.txt
>seq1
AAAAA
AAAAA
>seq2
CCCCC
CCCCC
CCCCC
>seq3
TTTTT

I want to retrieve only seq1 and seq2 (listed on list.txt) and redirect them to individual files. As you can see, each sequence has different number of lines hence I cannot just say to 'sed' to pick up N number of lines after each match. I want my output like this:

$ ls
seq1.txt
seq2.txt

$ cat seq1.txt
>seq1
AAAAA
AAAAA
$ cat seq2.txt
>seq3
TTTTT

I am using this:

while read list
do
names=$(echo $list)
        sed '/$list/,/>/{/>/q}' "$PWD/sequences.txt" > "$names".dna
done < list.txt

However, the output is:

$ ls
seq1.txt
seq2.txt

$ cat seq1.txt
>seq1
AAAAA
AAAAA
>seq3
TTTTT

$ cat seq2.txt
>seq1
AAAAA
AAAAA
>seq3
TTTTT

The script is creating individual files but all contain all the matches, not individuals as I need.

Thanks in advance.

  • Shell variables don't expand inside single quotes - so it's trying to match $list as a regex (which fails resulting in the default printing of the whole file)? – steeldriver Jun 8 '18 at 18:50
  • See fastaexplode from Guy Slater's exonerate suite of bioinformatics tools. Also, the Bioinformatics StackExchange site. – Kusalananda Jun 8 '18 at 19:02
0

Can't do this with sed alone. But with awk:

awk '
    # remember the wanted sequences
    NR == FNR {seqs[$1]; next}

    $1 ~ /^>/ {
        # get the sequence name
        seq = $1
        sub(/^>/, "", seq)
        p = 0
        # if it is in the list, set up the file to print to
        if (seq in seqs) {
            f = seq ".txt"
            p = 1
        }
    }
    p {print > f}
' list.txt sequences.txt 
0
command

awk '$1 ~ /seq1/{f=1}$1~/seq2/{f=0;exit}f' sequences.txt  >seq1.txt

sed -n '/seq3/,$p' sequences.txt > seq2.txt

output

cat seq1.txt
seq1
AAAAA
AAAAA

cat seq2.txt
seq3
TTTTT

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.