5
$\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
$\endgroup$
  • 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 '18 at 10:38
  • $\begingroup$ I have written all the header names in my question. $\endgroup$ – charlesdarwin Jun 21 '18 at 11:25
  • $\begingroup$ There are 4 matches to your pattern "chr_I" $\endgroup$ – benn Jun 21 '18 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$ – charlesdarwin Jun 21 '18 at 12:10
9
$\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
$\endgroup$
  • 1
    $\begingroup$ Worked. I just don't know sed. I need to learn it. $\endgroup$ – charlesdarwin Jun 21 '18 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 '18 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 '18 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 '18 at 17:21
  • 2
    $\begingroup$ @b.nota Duh! I should learn to read :) Indeed, sed 's/^>chr_/>/ would be better. $\endgroup$ – terdon Jun 21 '18 at 22:09
1
$\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.

$\endgroup$
0
$\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

$\endgroup$
  • $\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 '18 at 10:20

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.