7
$\begingroup$

I am trying to use the seqkit replace command to replace chromosome names in the format chr_I, chr_II, ... to I, II, .... I am using the following command:

seqkit replace -p "(.)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

My keyValues.txt file contains the following:

chr_I   I
chr_II  II
...

The two columns are separated by a tabulation.

I get the following output:

[INFO] read key-value file: keyValues.txt
[INFO] 6 pairs of key-value loaded
[ERRO] pattern "(.)" matches multiple targets in "chr_I", this will cause chaos

I don't know what I am doing wrong. I chose the (.) pattern to match the whole header but it seems to be wrong. Any help would be appreciated.

Update

All the headers are shown below:

grep '>' mySequence.fasta
>chr_I
>chr_II
>chr_III
>chr_IV
>chr_V
>chr_X
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Show some of the headers, the pattern argument is probably not suitable. I am guessing that you'll get headers with duplicate names now (which will cause chaos). $\endgroup$
    – benn
    Jun 21, 2018 at 10:38
  • $\begingroup$ I have written all the header names in my question. $\endgroup$
    – Biomagician
    Jun 21, 2018 at 11:25
  • $\begingroup$ There are 4 matches to your pattern "chr_I" $\endgroup$
    – benn
    Jun 21, 2018 at 11:38
  • $\begingroup$ Oh, so the key is a regular expression? I tried putting a dollar at the end of the key names, but I get the same error. $\endgroup$
    – Biomagician
    Jun 21, 2018 at 12:10

4 Answers 4

Reset to default
10
$\begingroup$

Why not use sed instead?

sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

Or even simpler:

sed 's/chr_//' mySequence.fasta > mySeq.fasta
Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Worked. I just don't know sed. I need to learn it. $\endgroup$
    – Biomagician
    Jun 21, 2018 at 12:12
  • 1
    $\begingroup$ Great! It works a bit like grep (which you know how to use), but then with replacement function. $\endgroup$
    – benn
    Jun 21, 2018 at 12:16
  • 1
    $\begingroup$ I would suggest using sed 's/^chr_//' instead. Anchoring the pattern to the beginning of the line should speed things up for large files. Otherwise, sed will have to scan the entire length of each line in case there's a match. $\endgroup$
    – terdon ♦
    Jun 21, 2018 at 16:00
  • 1
    $\begingroup$ @terdon, good point. But does it work with fasta? Remember the > at the beginning of the header. $\endgroup$
    – benn
    Jun 21, 2018 at 17:21
  • 2
    $\begingroup$ @b.nota Duh! I should learn to read :) Indeed, sed 's/^>chr_/>/ would be better. $\endgroup$
    – terdon ♦
    Jun 21, 2018 at 22:09
2
$\begingroup$

I like using awk when there's any kind of lookup involved: something like

awk 'FNR==NR {
    hash[">"$1]=$2;
    next;
}

if ($1 ~ /^>/) {
    print ">"hash[$1];
} else {
    print $1;
}' keyValues.txt mySequence.fasta

Basically the FNR==NR check tells awk to work only on the first file (useful primer here), keyValues.txt and create an association for each key (eg. chr_I) with its value (eg. I). The rest of the code after the next works only on mySequence.fasta, printing out the lookup value only if the line is a fasta header, as checked by the $1 ~ /^>/ condition.

Improve this answer
$\endgroup$
1
$\begingroup$

Here is an alternative to sed: cat test | tr -d "chr_"

The pipe is to output the file, while the tr command is to translate and change the file

Improve this answer
$\endgroup$
1
  • $\begingroup$ The use of the cat command is completely redundant here (and imposes a runtime overhead). Just do tr -d 'chr_' < test. $\endgroup$
    – Konrad Rudolph
    Jun 25, 2018 at 10:20
1
$\begingroup$

Your issue is that the pattern . matches a single character. To match the whole header, use .*:

seqkit replace -p "(.*)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

Alternatively as pointed out by others you may not really need keyValues.txt here as you're just removing the chr_ prefix:

seqkit replace -p "chr_(.*)" --replacement '$1' mySequence.fasta

(Don't forget the single quotes for $1 else the shell will expand it)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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