Question: Hello everyone,

I got a few flanking sequences around a bunch of SNPs. Since they were export from biomart, the variant allele was like a _:

>2|28804|28804|rs544393009|dbSNP|C/T
TGCAT_TCAGC
>2|29358|29358|rs943169482|dbSNP|G/GAAAA
CCCAA_AAGCC
>2|28396|28396|rs781239227|dbSNP|T/C
AATGC_CTTGG

Can some experts help me to replace the _ in each fasta sequence with its header description? That is each fasta input will be turned in 2 or 3 flanking sequence:

>2|28804|28804|rs544393009|dbSNP|C
TGCATCTCAGC
>2|28804|28804|rs544393009|dbSNP|T
TGCATTTCAGC
And that is what I am desperately seeking for!!!

Thanks in advance!

Answer: That's not that hard by matching the SNPs and flanking sequences using a regular expression.

To make it easy to process, let's linearize the FASTA sequences,
i.e. converting to tabular format. For example:

$ seqkit fx2tab seqs.fa
2|28804|28804|rs544393009|dbSNP|C/T TGCAT_TCAGC

Hereby, we can follow the original answer with little modification.
Note that the 1 or 2 before allels are used to keeping order, which will be
removed in the end.

# FASTA -> tabular
# replace 'C/T\tTGCAT_TCAGC' with '1C\tTGCATCTCAGC'
# tabular -> FASTA
$ seqkit fx2tab seqs.fa | \
perl -pe 's/(\w+)\/(\w+)\t(\w+)_(\w+)/1$1\t$3$1$4/' | \
seqkit tab2fx > seqs.allel1.fa

# FASTA -> tabular
# replace 'C/T\tTGCAT_TCAGC' with '2T TGCATTTCAGC'
# tabular -> FASTA
$ seqkit fx2tab seqs.fa | \
perl -pe 's/(\w+)\/(\w+)\t(\w+)_(\w+)/2$2\t$3$2$4/' | \
seqkit tab2fx > seqs.allel2.fa

Then sort records in the two files by sequence ID in alphabetical order

# sort by ID
# remove the 1/2 before alleles
$ seqkit sort seqs.allel*.fa |\
perl -pe 's/\|\d(\w+)$/\|$1/'
>2|28396|28396|rs781239227|dbSNP|T
AATGCTCTTGG
>2|28396|28396|rs781239227|dbSNP|C
AATGCCCTTGG
>2|28804|28804|rs544393009|dbSNP|C
TGCATCTCAGC
>2|28804|28804|rs544393009|dbSNP|T
TGCATTTCAGC
>2|29358|29358|rs943169482|dbSNP|G
CCCAAGAAGCC
>2|29358|29358|rs943169482|dbSNP|GAAAA
CCCAAGAAAAAAGCC

#### Original answer:

That's not that hard by matching the SNPs and flanking sequences using a regular expression.

# replace 'C/T TGCAT_TCAGC' with 'C TGCATCTCAGC'
$ perl -pe 's/(\w+)\/(\w+) (\w+)_(\w+)/$1 $3$1$4/' seqs.fa > seqs.allel1.fa
# replace 'C/T TGCAT_TCAGC' with 'T TGCATTTCAGC'
$ perl -pe 's/(\w+)\/(\w+) (\w+)_(\w+)/$2 $3$2$4/' seqs.fa > seqs.allel2.fa

If the order between two allels is not important, just sort them by sequence ID in alphabetical order:

$ seqkit sort seqs.allel*.fa
>2|28804|28804|rs544393009|dbSNP|C TGCATCTCAGC
ACTGN
>2|28804|28804|rs544393009|dbSNP|T TGCATTTCAGC
ACTGN
>2|29358|29358|rs943169482|dbSNP|G CCCAAGAAGCC
actgn
>2|29358|29358|rs943169482|dbSNP|GAAAA CCCAAGAAAAAAGCC
actgn

If the order is important, it can be done with a trick.

# replace 'C/T TGCAT_TCAGC' with '1C TGCATCTCAGC'
$ perl -pe 's/(\w+)\/(\w+) (\w+)_(\w+)/1$1 $3$1$4/' seqs.fa > seqs.allel1.fa
# replace 'C/T TGCAT_TCAGC' with '2T TGCATTTCAGC'
$ perl -pe 's/(\w+)\/(\w+) (\w+)_(\w+)/2$2 $3$2$4/' seqs.fa > seqs.allel2.fa

# remove the 1/2 before allels after sort
$ seqkit sort seqs.allel*.fa | perl -pe 's/\|\d(\w+) /\|$1 /'
>2|28804|28804|rs544393009|dbSNP|C TGCATCTCAGC
ACTGN
>2|28804|28804|rs544393009|dbSNP|T TGCATTTCAGC
ACTGN
>2|29358|29358|rs943169482|dbSNP|G CCCAAGAAGCC
actgn
>2|29358|29358|rs943169482|dbSNP|GAAAA CCCAAGAAAAAAGCC
actgn